1、 物理系 2007 级本科生王 峰 0710130948辛瑞平 0710130944指导老师：赵虎廖红波【摘要】本文将介绍使用 matlab 软件将核磁共振（NMR）实验做可视化，旨在指导以后的同学们更好的经行近代物理实验的学习和操作【关键字】核磁共振 弛豫时间的测量 matlab1939 年美国哥伦比亚的 II拉比（Rabi）等人在原子束实验中首次观察到核磁共振（Nuclear Magnetic Resonance,简称 NMR）现象。1945 年 12 月，美国哈佛大学的珀塞尔（Puccell 生于 1912 年）等人，报道了他们在石蜡样品中观察到质子的核磁共振吸收信号；1946 年 1
2、月，美国斯坦福大学布洛赫（Bloch 生于 1905 年）等人，也报道了他们在水样品中观察到质子的核感应信号。两个研究小组用了稍微不同的方法，几乎同时在凝聚物质中发现了核磁共振。因此，布洛赫和珀塞尔荣获了 1952 年的诺贝尔物理学奖。一 理论分析从经典理论观点来讨论核磁共振问题。把经典理论核矢量模型用于微观粒子是不严格的，但是它对某些问题可以做一定的解释。数值上不一定正确，但可以给出一个清晰的物理图象，帮助我们了解问题的实质。所以我们模拟的行为是明显而清晰地经典观点。现在将我们地工作一一介绍如下：1 单个核的拉摩进动如果陀螺不旋转，当它的轴线偏离竖直方向时，在重力作用下，它就会倒下来。但是如
3、果陀螺本身做自转运动，它就不会倒下而绕着重力方向做进动，如图 1 所示。图 1 陀螺进动由于原子核具有自旋和磁矩，所以它在外磁场中的行为同陀螺在重力场中的行为是完全一样的。设核的角动量为 ，磁矩为 ，外磁场为 ，由经典理论可知PB（1）dt由于， ，所以有（2）Bdt)()(xyxz zzyyzyxBdtdt若设稳恒磁场为 ，且 轴沿 方向，即 ， ，则上式将变为000yx0Bz（3）00dtBdtzxyyx这样我们可以利用解微分方程组的办法将磁矩的三个分量解出，同时以分量作为末了点的坐标，从而可以画出单个核的拉摩进动过程。下面我们只依照解析解的情况对于这一结果做一定性的描述。以期与运行结果对
4、比。从上面的方程看出，磁矩分量 是一个常数（ ） ，即磁矩 在 方向上的z0dtz0B投影将保持不变。以 代入，有00B（4）常 数AtyxLyx20)(sinco由此可知，核磁矩 在稳恒磁场中的运动特点是：（1）它围绕外磁场 做进动，进动的角频率为 ，和 与 的夹角 无0B00B0关；（2）它在 平面上的投影 是常数；xyL（3）它在外磁场 方向上的投影 为常数。0z其运动图像如图 2 所示。图 2图 3如果在与 垂直的方向上加一个旋转磁场 ，且 ， 的角频率和转动方0B1B01B向与磁矩 的进动角频率和进动方向都相同，如图 3 所示。这时，和核磁矩 除了受到 的作用之外，还要受到旋转磁场
5、的影响。也就是说 除了要围绕 进动之外，还要0 10绕 进动。1B这样上面的（2）方程组就将有如下变化： ，)cos(),sin(11 tBtByx ，0z（5） )sin()cos()(in)cos()( 11010tBtBdt tdt yxxyxz xzzzy zyyzyx 所以 与 之间的夹角 将发生变化。由核磁矩的势能0（6） cos0E可知， 的变化意味着核的能量状态变化。当 值增加时，核要从旋转磁场 中吸收能量。1B这就是核磁共振。产生共振的条件为（7）00B如果旋转磁场 的转动角频率 与核磁矩 的进动角频率 不相等，即 ，100则角度 的变化不显著。平均说来， 角的变化为零。原子
6、核没有吸收磁场的能量，因此就观察不到核磁共振信号。2 弛豫过程与 Bloch 方程2.1 体磁化强度上面讨论的是单个核的核磁共振。但我们在实验中研究的样品不是单个核磁矩，而是由这些磁矩构成的体磁化强度矢量 ；另外，我们研究的系统并不是孤立的，而是与周围M物质有一定的相互作用。只有全面考虑了这些问题，才能建立起核磁共振的理论。因为体磁化强度矢量 是单位体积内核磁矩 的矢量和， 即：（8）NiM1由于（2）式，可得：（9）)(Bdt它表明磁化强度矢量 围绕着外磁场 做进动，进动的角频率 。而处在0B 00B热平衡状态时，体磁化强度 沿 Z 方向，记为 。M2.2 射频脉冲磁场 瞬态作用1B引入一个
7、旋转坐标系 ， 方向与 方向重合，坐标旋转角频率 ,则),(zyx0 0M在新坐标系中静止。在 轴上加恒定磁场 ,作用时间为脉宽 pt，在射频脉冲作用前1处在热平衡状态，方向与 轴（ 轴）重合（即 ） ，施加射频脉冲作用，则 将以z0MM频率 绕 轴进动。1Bx图 4M转过的角度 （如图 4（a） ）称为倾倒角,如果脉冲宽度恰好使ptB1或 ，称这种脉冲为 09或 8脉冲。 09脉冲作用下 M将倒在 上，2/ y018脉冲作用下将倒向 方向。由 可知，只要射频场足够强，则 值均可zptB1pt以做到足够小而满足 （ 为下面将要介绍的弛豫时间）,这意味着射频脉冲21,Ttp作用期间弛豫作用可以忽
8、略不计。这一过程我们称之为射频脉冲磁场 的瞬态作用。1B2.3 弛豫过程弛豫和射频诱导激发是两个相反的过程，当两者的作用达到动态平衡时，实验上可以观测到稳定的共振讯号。弛豫因涉及到体磁化强度的纵向分量和横向分量变化，故分为纵向弛豫和横向弛豫。2.3.1 纵向弛豫（ ）01:MT纵向弛豫又称为自旋晶格弛豫。宏观样品是由大量小磁矩的自旋系统和它们所依附的晶格系统组成。系统间不断发生相互作用和能量变换，纵向弛豫是指自旋系统把从射频磁场中吸收的能量交给周围环境，转变为晶格的热能。自旋核由高能态无辐射地返回低能态，能态粒子数差 按下式规律变化n（10）)/exp(10Tt式中， 为时间 时的能态粒子差，
9、 为粒子数的差异与体磁化强度 M的纵向0t1分量 的变化一致，粒子数差增加 也相应增加，故 称为纵向弛豫时间。ZMZM1T是自旋体系与环境相互作用时的速度量度， 的大小主要依赖于样品核的类型和样1T品状态，所以对 的测定可知样品核的信息。12.3.2 横向弛豫（ ）0:2MT横向弛豫又称为自旋自旋弛豫。自旋系统内部也就是说核自旋与相邻核自旋之间进行能量交换，不与外界进行能量交换，故此过程体系总能量不变。自旋自旋弛豫过程，由非平衡进动相位产生时的体磁化强度的横向分量 0 恢复到平衡态时相位无关=0 表征,所需的特征时间记为 。由于 与体磁化强度的横向分量 的弛豫时间有M2T2 M关，故 也称横向
10、弛豫时间。自旋自旋相互作用也是一种磁相互作用，进动相位相关主2T要来自于核自旋产生的局部磁场。射频场 ，外磁场空间分布不均匀都可看成是局部磁场。1B2.4 脉冲作用后体磁化强度 M的行为自由感应衰减（FID）信号介绍完上面的知识，我们来分析一个简单的过程，这将作为下面我们演示的基础。设 时刻加上射频场 ,到 时绕 旋转 而倾倒在 轴上，这时射频0t1Bpt1B09y场 消失，核磁矩系统将由弛豫过程回复到热平衡状态。由上面所述其中 的变1B 0zM化速度取决于 , 和 的衰减速度取决于 ,在旋转坐标系看来，没有1TxM0y 2T进动，恢复到平衡位置的过程如图 5（a）所示。在实验室坐标系看来，
11、M绕 轴旋进按z螺旋形式回到平衡位置，如图 5（b）所示。图 5在这个弛豫过程中，若在垂直于 轴方向上置一个接收线圈，便可感应出一个射频信z号，其频率与进动频率 相同，其幅值按照指数规律衰减，称为自由感应衰减信号，也写0作 FID 信号。经检波并滤去射频以后，观察到的 FID 信号是指数衰减的包络线，如图4（a）所示。FID 信号与 M在 平面上横向分量的大小有关，所以 脉冲的 FID 信号xy 09幅值最大， 脉冲的幅值为零。018图 6 自由感应衰减信号 实验中由于恒定磁场 不可能绝对均匀，样品中不同位置的核磁矩所处的外场大小有0B所不同，其进动频率各有差异，实际观测到的 FID 信号是各
12、个不同进动频率的指数衰减信号的叠 加，如图 6（b）所示，设 为磁场不均匀所等效的横向弛豫时间，则总的 FID 信2T号的衰减速度由 和 两者决定，可以用一个称为表观横向弛豫时间 来等效：2T 2T（11）2*21若磁场域不均匀，则 越小，从而 也越小，FID 信号衰减也越快。T22.5 Bloch 方程前面分别分析了外磁场和弛豫过程对核磁化强度矢量 的作用。当上述两种作用同时M存在时，描述核磁共振现象的基本运动方程为（12）kTjiMTBdtMzyx 102)(1)( 该方程称为布洛赫方程。式中 ， ， 分别是 ， ， 方向上的单位矢量。将其展开ijkxz为：（13）102)()(TMBdt
13、MTdt zxyxz yzxzy xZzyx 3 自旋回波法测量横向弛豫时间 2T（ 00189脉冲序列方式）自旋回波是一种用双脉冲或多个脉冲来观察核磁共振信号的方法，它特别适用于测量横向弛豫时间 2T，谱线的自然线宽是由自旋自旋相互作用决定的，但在许多情况下，由于外磁场不够均匀，谱线就变宽了，与这个宽度相对应的横向弛豫时间是前面讨论过的表观横向弛豫时间 ，而不是 了，但用自旋回波法仍可以测出横向弛豫时间 。22 2T图 7实际应用中，常用两个或多个射频脉冲组成脉冲序列，周期性的作用于核磁矩系统。比如在 09射频脉冲作用后，经过 时间再施加一个 射频脉冲，便组成一个018脉冲序列，这些脉冲序列
14、的脉宽 和脉距 应满足下列条件：18 pt,21Ttp（14）*2,脉冲序列的作用结果如图 5 所示，在 射频脉冲后即观察到 FID 信0089 09号；在 射频脉冲后面对应于初始时刻的 处可以观察到一个“回波”信号。这种回12波信号是在脉冲序列作用下核自旋系统的运动引起的，所以称为自旋回波。 ()()() +() ()图 8 自旋回波矢量图解0019以下用图 8 来说明自旋回波的产生过程。图 8（a）表示体磁化强度 在 射频脉冲作用下绕 轴转到 轴上；0Mxy图 8（b）表示脉冲消失后核磁矩自由进动受到 不均匀的影响，样品中部分磁矩的0B进动频率不同，引起磁矩的进动频率不同，使磁矩相位分散并
15、呈扇形展开。为此可把 看M成是许多分量 之和。从旋转坐标系看来，进动频率等于 的分量相对静止，大于i 0的分量（图中以 代表）向前转动，小于 的分量（图中以 为代表）向后转动；01M02图 8（c）表示 射频脉冲的作用使磁化强度各分量绕 轴翻转 ，并继续它们0 x018原来的转动方向运动；图 8（d）表示 时刻各磁化强度分量刚好汇聚到 轴上；2t y图 8（e）表示 以后，用于磁化强度各矢量继续转动而又呈扇形展开。因此，在2t处得到如图 7 所示的自旋回波信号。2t由此可知，自旋回波与 FID 信号密切相关，如果不存在横向弛豫，则自旋回波幅值应与初始的 FID 信号一样，但在 时间内横向弛豫作
16、用不能忽略，体磁化强度各横向分量相应减小，使得自旋回波信号幅值小于 FID 信号的初始幅值，而且脉距 越大则自旋回波幅值越小，并且回波幅值 与脉距 存在以下关系：U（15）2/0Tte式（15）中 ， 是 射频脉冲刚结束时 FID 信号的初始幅值，实验中只要t09改变脉距 ，则回波的峰值就相应的改变，若依次增大 测出若干个相应的回波峰值，便得到指数衰减的包络线。对(15)式两边取对数，可以得到直线方程（16）20/lnTU上式中 作为自变量，则直线斜率的倒数便是 。22这就是我们测 的方法。2T二 程序设计基于上面所述的三个过程我们进行了编程可视化现逐一介绍如下：1 单个核的拉摩进动运行结果：
17、有两条黄色箭头代表 y 与 z 方向的磁场。红色箭头为总磁场，而绿色的场线束代表旋转着的磁矩。右侧的三个下拉菜单分别控制：面板上设有三个参数 0B， 1， B可以进行调解控制动画的变化。拖动条控制视角。设有停止与退出按钮。点击 plot 按钮后可以观察到在实际情况中单条磁矩的拉莫尔进动。另外值得一提的是此面板的设计是完全靠代码实现的，没有依靠 GUI 工具箱，相对比较简洁。2 弛豫过程（自由感应衰减（FID）信号）-1 -0.50 0.51 -1-0.500.5 100.511.52中中中中中中中中-1 -0.5 0 0.5 100.20.40.60.811.21.41.61.82 中中中中中
18、中中中中中中中0 1 2 3 4 5 6 7 8 9 10-1-0.500.51中中中中中中中此面板依旧是采用程序代码设计。分为三个绘图区，分别为：磁矩尖端的实际运动轨迹（左上） ；磁矩大小随时间的衰减（右上） ；与所测到的自感衰减信号随时间变化的情况（左下） 。点击“绘图”按钮后，会随着时间三幅图同时运作，展示出对应的变化关系。 3 自旋回波法测量横向弛豫时间 2T（ 00189脉冲序列方式）界面简洁，功能全面。左侧为 运动模拟区 ：描绘若干条磁矩在旋转系中的运行行为。运动模拟区下方为 值改变条 ，通过拉动其来改变实验的 值。 右侧上方为 信号绘图区 ，为在动系 y 轴的方向设置探测器来接收
19、辐射的电测波信号，横轴为时间；纵轴为所测得的信号幅值。信号绘图区下方为 数据记录区 ，分为：序号 （记录改变 值，并运行次数。 ） ；脉冲长度 （读取通过拉动条改变实验的 值） ；回波强度最大值 （即公式（15）中的 值，这样只要测出不同 值下的 值，利用UU公式（16）就能够通过线性拟合，求出 的值。完全按照实验的操作方式。另外附加了一2T个功能，即遇到相同的 值时，不记录实验数据）在数据记录区下面是 状态区 ，其作用是提示用户当前程序的运行状态，指导操作。在面板的最下部设计四个 控制按钮 ：运行： 取拉动条所取的 值进行模拟；数据分析： 打开数据分析面板，对所得的实验数据，经行数据拟合等功
20、能；清零： 将所有实验数据清零；退出： 关闭此面板。上面是数据处理面板。试将前面所的到的实验数据经行拟合，从而得到 的值。左侧为2T绘图区；右侧为结果显示区，即将所得到的拟合后公式显示出来，并展示所求得的 的值。2下面为具体的操作实例：三 总结与感想本实验模拟能够基本上满足实验教学的需要，且图像清晰，形象直观，操作简便。可以作为近代物理实验教学时的课件，或者给同学们提供实验预习之用。但由于需要将方程完全解出，才能一次性将图画出，所以解方程的时间较长是本可视化的一大缺憾可以在以后的学习中继续体会，试图找出能够缩短运行时间的办法。这次做这个实验的可视化可以说是十分艰难，我们从开始对近代物理实验一无
21、所知到题目的完成，很大程度上是老师耐心指导的结果。要做好一个物理实验可视化，首先要对这个实验有相当程度上的了解，而且在图像和原理上不能含糊，在这一点上尤其感谢廖红波老师的耐心指导和指正，让我们的物理概念得到了澄清，物理图像更加清晰。感谢赵虎老师给我们的这么有意思的题目。使我们在做的过程中学到了太多课程上所不能学到的东西，希望在未来的日子里还能有机会做类似的可视化，在提高我们才能的同时，为我们系的实验教学做自己的一点点贡献。 参考文献：1FD-PNMR- 型脉冲核磁共振实验仪实验指导书2 FD-PNMR-II 脉冲核磁共振仪说明书3近代物理实验.熊俊.北京师范大学出版社.2007 年 8 月第一
22、版4精通 MATLAB GUI 设计.陈垚光 毛涛涛 王正林 王玲编著. 北京，电子工业出版社附录：具体的 M 文件1. 单个核的拉摩进动%*%核磁共振主函数%*function HCGZ_MAINclcclear global a b m %共享变量hf=figure(Position,100,80,800,550,Name,核磁共振,.Color,.2 .1 .1,NumberTitle,off);ha=axes(Position,0.065 0.095 0.65 0.85,Color,.2 .1 .15,.%设置画图区域XColor,w,YColor,w,ZColor,w,Box,off
23、);hviewb=uicontrol(hf,Style,slider,Position,10,53,10,468,.%纵向视角控制句柄Value,-30,Min,-90,Max,0,.Callback,b=get(hviewb,Value););hviewa=uicontrol(hf,Style,slider,Position,52,15,520,10,.%横向视角控制句柄Value,-40,Min,-90,Max,0,.Callback,a=get(hviewa,Value););uicontrol(hf,Style,text,Position,650 450 70 15,String,定磁
24、场 B0); %文本框hB0=uicontrol(hf,Style,popupmenu,Position,650,420,100,30,.%选择定磁场下拉菜单String,1|1.5|2, Value,3);uicontrol(hf,Style,text,Position,650 390 70 15,String,旋转磁场 B1); hB1=uicontrol(hf,Style,popupmenu,Position,650,360,100,30,.%选择旋转磁场下拉菜单String,0.2|0.4|0.6|0.8|1|0.001, Value,3);uicontrol(hf,Style,tex
25、t, Position,650 330 100 15,String,旋转磁场角速度 wB1); hwB1=uicontrol(hf,Style,popupmenu,Position,650,300,100,30,.%选择角速度下拉菜单String,1|2|3|4|5, Value,3); hexit=uicontrol(hf,Style,pushbutton,. %退出按钮Position,700,10,80,30, String,EXIT,.CallBack,m=1;,close(hf);hstop=uicontrol(hf,Style,pushbutton,. %停止按钮Position,
26、650,130,100,30,String,STOP,.Callback,m=1;);hplot=uicontrol(hf,. %作图按钮Style,pushbutton,.Position,650,200,100,30,String,PLOT,.Callback,m=1;.B0UD=1 1.5 2;.B0=B0UD(get(hB0,Value);.B1UD=0.2 0.4 0.6 0.8 1 0.001;. B1=B1UD(get(hB1,Value);.wB1UD=1 2 3 4 5;.wB1=wB1UD(get(hwB1,Value);.m=0;,.hcgz(B0,B1,wB1););%
27、调用函数%*%核磁共振被调用函数%*function hcgz(B0,B1,wB1)global a b mw=wB1;r=15;F=(t,y)r*(B1*sin(w*t)*y(3)-B0*y(2); %微分方程r*(B0*y(1)-B1*cos(w*t)*y(3);r*(B1*cos(w*t)*y(2)-B1*sin(w*t)*y(1);ux=(3(1/2)*(B0/3 + (2(1/2)*B1)/3)/(B02 + B12)(1/2);uz=sqrt(1-ux2); %磁矩的初始值T,Y=ode45(F,0:0.005:5,ux,0,uz); %求解方程Bx=B1*cos(w*T);By=
28、B1*sin(w*T);Bz=B0*ones(length(T),1);%计算旋转磁场随时间的变化B=Bx,By,Bz;B1=Bx,By,0*Bz;grid onhold ona=-30;b=-40;a1=-a;b1=-b;view(a1,b1); %设置视角axis(-1,1,-1,1,0,1.5);N=50;qqq=zeros(N,1);h0=quiver3(0,0,0,0,0,B0,y); %画定磁场 B0h1=quiver3(qqq,qqq,qqq,Y(1:N,1),Y(1:N,2),Y(1:N,3); %画磁矩的变化h2=quiver3(0,0,0,B(1,1),B(1,2),B(1
29、,3); %画旋转磁场h3=quiver3(0,0,0,B1(1,1),B1(1,2),B1(1,3); %画旋转磁场在水平方向的投影h4=plot3(0.9*B(1,1) 0.9*B1(1,1),0.9*B(1,2) 0.9*B1(1,2),0.9*B(1,3) 0.9*B1(1,3),w:);h5=plot3(0.9*B(1,1) 0,0.9*B(1,2) 0,0.9*B(1,3) 0.9*B0,w:);h6=plot3(0.9*Y(1:N,1),0.9*Y(1:N,2),0.9*Y(1:N,3),c);set(h6,EraseMode,xor);while 1 %用擦图的方法画动画if
30、m=1delete(h0,h1,h2,h3,h4,h5,h6); break; %设置停止开关endfor i=N+1:length(T)if m=1break; %设置停止开关enda1=-a;b1=-b;view(a1,b1);delete(h1,h2,h3,h4,h5)h1=quiver3(qqq,qqq,qqq,Y(i-N:i-1,1),Y(i-N:i-1,2),Y(i-N:i-1,3),g);h2=quiver3(0,0,0,B(i,1),B(i,2),B(i,3),r);h3=quiver3(0,0,0,B1(i,1),B1(i,2),B1(i,3),y);h4=plot3(0.9
31、*B(i,1) 0.9*B1(i,1),0.9*B(i,2) 0.9*B1(i,2),0.9*B(i,3) 0.9*B1(i,3),c:);h5=plot3(0.9*B(i,1) 0,0.9*B(i,2) 0,0.9*B(i,3) 0.9*B0,c:);set(h6,Xdata,0.9*Y(i-N:i,1),Ydata,0.9*Y(i-N:i,2),Zdata,0.9*Y(i-N:i,3);drawnow ;endendhold off2 弛豫过程（自由感应衰减（FID）信号）function CYGC_MAINclcclear global ha hb hchf=figure(Positio
32、n,80,80,850,600,.Name,核磁共振弛豫过程,.Color,.8 0.9 .9,.NumberTitle,off);ha=axes(Position,0.04 0.38 0.45 0.55,. %设置画图区 1Color,.8 0.9 .9);hb=axes(Position,0.53 0.38 0.45 0.55,. %设置画图区 2Color,.8 0.9 .9);hc=axes(Position,0.1 0.05 0.55 0.25,. %设置画图区 3Color,.8 0.9 .9);hdo=uicontrol(hf,. %设置画图按钮Style,pushbutton,
33、.Position,680,80,100,50,String,绘 图,.Callback,CYGC); %调用函数function CYGC(cbo,eventdata,h)global ha hb hcB1=0 ;B0=5 ;w=0;r=8;T1=15;T2=2;F=(t,y)r*(B1*sin(w*t)*y(3)-B0*y(2)-y(1)/T2; r*(B0*y(1)-B1*cos(w*t)*y(3)-y(2)/T2;r*(B1*cos(w*t)*y(2)-B1*sin(w*t)*y(1)-y(3)/T1+0.2;T,Y=ode45(F,0:0.005:10,0,1,0); %求解微分方程
34、得弛豫过程磁矩的变化数据subplot(ha) %在 ha 区域画静止系中的图像hold ongrid onhha=plot3(Y(1,1),Y(2,1),Y(1,2),Y(2,2),Y(1,3),Y(2,3);%画磁矩hha2=plot3(0,Y(1,1),0,Y(1,2),0,Y(1,3),r); %画磁矩的 Z 分量axis(-1.2,1.2,-1.2,1.2,0,2);view(30,40);title(磁矩随时间的变化);hold offsubplot(hb) %在 hb 区域画静止系中的图像a=Y(1,2);b=Y(1,3);hold ongrid onhhb=plot(Y(1,2
35、),Y(1,3);%画磁矩axis(-1.2,1.2,0,2);title(旋转系中磁矩随时间的变化);hold offsubplot(hc)%在 hc 区域画自感应衰减信号hold ongrid onhhc=plot(T(1),Y(1,1);%画自感应衰减信号hhc1=plot(T(1),exp(-T(1)/T2),r);%画上包络线hhc2=plot(T(1),-exp(-T(1)/T2),r);%画下包络线axis(T(1),T(end),-1.2,1.2);title(自感应衰减信号);hold offset(hha,EraseMode,none);set(hha2,EraseMode
36、,xor);set(hhb,EraseMode,none);set(hhc,EraseMode,none);set(hhc1,EraseMode,none);set(hhc2,EraseMode,none);for i=1:length(T)-1 %用擦图的方法画动画图set(hha,Xdata,Y(i,1),Y(i+1,1),Ydata,Y(i,2),Y(i+1,2),Zdata,Y(i,3),Y(i+1,3);set(hha2,Xdata,0,Y(i,1),Ydata,0,Y(i,2),Zdata,0,Y(i,3); if Y(i,1)0 %判断是否有残留图像，若有则删除delete(ch
37、ildrenhandles);end;if mark=0 %mark=0 表示脉冲长度有变化，重新计算数据 tao=get(handles.slider1,value);%将滑动条的值赋值给变量 taokk=kk+1;num_taozhi(kk)=tao; %记录数据作为文本框要显示的内容num_xuhao(kk)=kk+0.0000001;numtaozhi=num2str(num_taozhi);%将数值转化为文本numxuhao=num2str(num_xuhao);set(handles.taozhi,string,numtaozhi);%设置文本框显示文本set(handles.xu
38、hao,string,numxuhao);ZXHB(tao); %调用函数计算数值endset(handles.zhuangtai,string,状态：正在画图.) %显示状态栏k=zeros(1,N); for i=1:length(T)-2 %for 循环实现动画显示if abs(T(i)-2*tao)0.01 %找出 t=tao 时刻的数据num_fenzhi(kk)=abs(sum(YYy(i,:)numfenzhi=num2str(num_fenzhi);set(handles.fenzhi,string,numfenzhi);endendfor i=1:length(T)-2axe
39、s(handles.axes3) %在 hb 中画旋转系中磁矩随时间的变化hh1=quiver3(k,k,k,YYx(i,:),YYy(i,:),Yz(i,:),k);axis(-0.8,0.8,-0.8,0.8,0,0.8);axes(handles.axes2) %在 hc 中画自旋回波信号hold onhh2=plot(T(i),T(i+1),sum(Yy(i,:),sum(Yy(i+1,:);axis(T(1),T(end),-11,11);drawnow;endset(handles.zhuangtai,string,状态：运行结束)function ZXHB(tao)global
40、YYx YYy Yx Yy Yz N T markr=90;T1=40; T2=3;N=11; %画 N 条磁矩B0=linspace(5.25,5.35,N); %不均匀磁场B=15,0,0;Yx=;Yy=;Yz=;T=;%加 90 度脉冲%*F1=(t,y)r*(B(2)*y(3)-B(3)*y(2);r*(B(3)*y(1)-B(1)*y(3);r*(B(1)*y(2)-B(2)*y(1);opts=odeset(events,f); %打开事件函数t0,y0,t90final=ode45(F1,0,Inf,0,0,1,opts);for i=1:Nyx=y0(:,1);yy=y0(:,
41、2);yz=y0(:,3);tt=t0;%90 度脉冲后%*F=(t,y)-r*B0(i)*y(2)-y(1)/T2;r*B0(i)*y(1)-y(2)/T2;-y(3)/T1+N/T1;t,y=ode23(F,tt(end):0.02:tt(end)+tao,yx(end,1),yy(end,1),yz(end,1);yx=yx;y(:,1);yy=yy;y(:,2);yz=yz;y(:,3);tt=tt;t;%加 180 度脉冲%*B=15*cos(r*5.3*tt(end),15*sin(r*5.3*tt(end),0; %动系中的磁场方向F1=(t,y)r*(B(2)*y(3)-B(3
42、)*y(2);r*(B(3)*y(1)-B(1)*y(3);r*(B(1)*y(2)-B(2)*y(1);t,y=ode23(F1,tt(end):0.0001:tt(end)+2*t90final,yx(end,1),yy(end,1),yz(end,1);yx=yx;y(:,1);yy=yy;y(:,2);yz=yz;y(:,3);tt=tt;t;%180 度脉冲后%*t,y=ode23(F,tt(end):0.02:tt(end)+2*tao,yx(end,1),yy(end,1),yz(end,1);yx=yx;y(:,1);yy=yy;y(:,2);yz=yz;y(:,3);tt=t
43、t;t;Yx=Yx,yx;Yy=Yy,yy;Yz=Yz,yz; %将脉冲过程中的数据统一到一起T=T,tt;endYYx=Yx.*cos(5.3*r*(T-t90final)+Yy.*sin(5.3*r*(T-t90final); %将静止系中的数据转换到旋转系YYy=Yy.*cos(5.3*r*(T-t90final)-Yx.*sin(5.3*r*(T-t90final);mark=1;function Stop,isterminal,direction=f(t,y) %事件函数Stop=y(3);isterminal=1;direction=;%*%数据拟合面板%*function var
44、argout = sjfx(varargin)gui_Singleton = 1;gui_State = struct(gui_Name, mfilename, .gui_Singleton, gui_Singleton, .gui_OpeningFcn, sjfx_OpeningFcn, .gui_OutputFcn, sjfx_OutputFcn, .gui_LayoutFcn, , .gui_Callback, );if nargin endif nargoutvarargout1:nargout = gui_mainfcn(gui_State, varargin:);elsegui_mainfcn(gui_State, varargin:);endfunction sjfx_OpeningFcn(hObject, eventdata, handles, varargin)handles.output = hObject;gui
