核磁共振的可视化.doc

物理系2007级本科生王 峰 0710130948辛瑞平 0710130944指导老师：赵虎 廖红波 【摘要】本文将介绍使用matlab软件将核磁共振（NMR）实验做可视化，旨在指导以后的同学们更好的经行近代物理实验的学习和操作【关键字】核磁共振 弛豫时间的测量 matlab1939年美国哥伦比亚的II拉比（Rabi）等人在原子束实验中首次观察到核磁共振（Nuclear Magnetic Resonance,简称NMR）现象。1945年12月，美国哈佛大学的珀塞尔（Puccell生于1912年）等人，报道了他们在石蜡样品中观察到质子的核磁共振吸收信号；1946年1月，美国斯坦福大学布洛赫（Bloch生于1905年）等人，也报道了他们在水样品中观察到质子的核感应信号。两个研究小组用了稍微不同的方法，几乎同时在凝聚物质中发现了核磁共振。因此，布洛赫和珀塞尔荣获了1952年的诺贝尔物理学奖。一 理论分析从经典理论观点来讨论核磁共振问题。把经典理论核矢量模型用于微观粒子是不严格的，但是它对某些问题可以做一定的解释。数值上不一定正确，但可以给出一个清晰的物理图象，帮助我们了解问题的实质。所以我们模拟的行为是明显而清晰地经典观点。现在将我们地工作一一介绍如下：1单个核的拉摩进动如果陀螺不旋转，当它的轴线偏离竖直方向时，在重力作用下，它就会倒下来。但是如果陀螺本身做自转运动，它就不会倒下而绕着重力方向做进动，如图1所示。图 1陀螺进动由于原子核具有自旋和磁矩，所以它在外磁场中的行为同陀螺在重力场中的行为是完全一样的。设核的角动量为，磁矩为，外磁场为，由经典理论可知 （1）由于，所以有 （2）若设稳恒磁场为，且轴沿方向，即，则上式将变为 （3）这样我们可以利用解微分方程组的办法将磁矩的三个分量解出，同时以分量作为末了点的坐标，从而可以画出单个核的拉摩进动过程。 下面我们只依照解析解的情况对于这一结果做一定性的描述。以期与运行结果对比。从上面的方程看出，磁矩分量是一个常数（），即磁矩在方向上的投影将保持不变。以代入，有 （4）由此可知，核磁矩在稳恒磁场中的运动特点是：（1）它围绕外磁场 做进动，进动的角频率为，和与的夹角无关；（2）它在平面上的投影是常数；（3）它在外磁场方向上的投影为常数。其运动图像如图 2所示。图 2图 3如果在与垂直的方向上加一个旋转磁场，且， 的角频率和转动方向与磁矩的进动角频率和进动方向都相同，如图3所示。这时，和核磁矩除了受到的作用之外，还要受到旋转磁场的影响。也就是说除了要围绕进动之外，还要绕进动。这样上面的（2）方程组就将有如下变化：， （5）所以与之间的夹角将发生变化。由核磁矩的势能 （6） 可知，的变化意味着核的能量状态变化。当值增加时，核要从旋转磁场中吸收能量。这就是核磁共振。产生共振的条件为 （7）如果旋转磁场的转动角频率与核磁矩的进动角频率不相等，即，则角度的变化不显著。平均说来，角的变化为零。原子核没有吸收磁场的能量，因此就观察不到核磁共振信号。2弛豫过程与Bloch方程2.1体磁化强度上面讨论的是单个核的核磁共振。但我们在实验中研究的样品不是单个核磁矩，而是由这些磁矩构成的体磁化强度矢量；另外，我们研究的系统并不是孤立的，而是与周围物质有一定的相互作用。只有全面考虑了这些问题，才能建立起核磁共振的理论。因为体磁化强度矢量是单位体积内核磁矩的矢量和， 即： （8） 由于（2）式，可得： （9）它表明磁化强度矢量围绕着外磁场做进动，进动的角频率。而处在热平衡状态时，体磁化强度沿Z方向，记为。2.2射频脉冲磁场瞬态作用引入一个旋转坐标系，方向与方向重合，坐标旋转角频率,则在新坐标系中静止。在轴上加恒定磁场,作用时间为脉宽，在射频脉冲作用前处在热平衡状态，方向与轴（轴）重合（即），施加射频脉冲作用，则将以频率绕轴进动。图 4 转过的角度（如图4（a）称为倾倒角,如果脉冲宽度恰好使或，称这种脉冲为或脉冲。脉冲作用下将倒在上，脉冲作用下将倒向方向。由可知，只要射频场足够强，则值均可以做到足够小而满足（为下面将要介绍的弛豫时间）,这意味着射频脉冲作用期间弛豫作用可以忽略不计。这一过程我们称之为射频脉冲磁场的瞬态作用。2.3弛豫过程弛豫和射频诱导激发是两个相反的过程，当两者的作用达到动态平衡时，实验上可以观测到稳定的共振讯号。弛豫因涉及到体磁化强度的纵向分量和横向分量变化，故分为纵向弛豫和横向弛豫。2.3.1纵向弛豫（）纵向弛豫又称为自旋晶格弛豫。宏观样品是由大量小磁矩的自旋系统和它们所依附的晶格系统组成。系统间不断发生相互作用和能量变换，纵向弛豫是指自旋系统把从射频磁场中吸收的能量交给周围环境，转变为晶格的热能。自旋核由高能态无辐射地返回低能态，能态粒子数差按下式规律变化 （10） 式中，为时间时的能态粒子差，为粒子数的差异与体磁化强度的纵向分量的变化一致，粒子数差增加也相应增加，故称为纵向弛豫时间。是自旋体系与环境相互作用时的速度量度，的大小主要依赖于样品核的类型和样品状态，所以对的测定可知样品核的信息。2.3.2横向弛豫（）横向弛豫又称为自旋自旋弛豫。自旋系统内部也就是说核自旋与相邻核自旋之间进行能量交换，不与外界进行能量交换，故此过程体系总能量不变。自旋自旋弛豫过程，由非平衡进动相位产生时的体磁化强度的横向分量0恢复到平衡态时相位无关=0表征,所需的特征时间记为。由于 与体磁化强度的横向分量的弛豫时间有关，故也称横向弛豫时间。自旋自旋相互作用也是一种磁相互作用，进动相位相关主要来自于核自旋产生的局部磁场。射频场，外磁场空间分布不均匀都可看成是局部磁场。2.4脉冲作用后体磁化强度的行为自由感应衰减（FID）信号介绍完上面的知识，我们来分析一个简单的过程，这将作为下面我们演示的基础。设时刻加上射频场,到时绕旋转而倾倒在轴上，这时射频场消失，核磁矩系统将由弛豫过程回复到热平衡状态。由上面所述其中的变化速度取决于,和的衰减速度取决于,在旋转坐标系看来，没有进动，恢复到平衡位置的过程如图5（a）所示。在实验室坐标系看来，绕轴旋进按螺旋形式回到平衡位置，如图5（b）所示。图 5在这个弛豫过程中，若在垂直于轴方向上置一个接收线圈，便可感应出一个射频信号，其频率与进动频率相同，其幅值按照指数规律衰减，称为自由感应衰减信号，也写作FID信号。经检波并滤去射频以后，观察到的FID信号是指数衰减的包络线，如图4（a）所示。FID信号与在平面上横向分量的大小有关，所以脉冲的FID信号幅值最大，脉冲的幅值为零。图 6自由感应衰减信号 实验中由于恒定磁场不可能绝对均匀，样品中不同位置的核磁矩所处的外场大小有所不同，其进动频率各有差异，实际观测到的FID信号是各个不同进动频率的指数衰减信号的叠加，如图6（b）所示，设为磁场不均匀所等效的横向弛豫时间，则总的FID信号的衰减速度由和两者决定，可以用一个称为表观横向弛豫时间来等效： （11）若磁场域不均匀，则越小，从而也越小，FID信号衰减也越快。2.5 Bloch方程前面分别分析了外磁场和弛豫过程对核磁化强度矢量的作用。当上述两种作用同时存在时，描述核磁共振现象的基本运动方程为 （12）该方程称为布洛赫方程。式中，分别是，方向上的单位矢量。将其展开为： （13）3 自旋回波法测量横向弛豫时间（脉冲序列方式）自旋回波是一种用双脉冲或多个脉冲来观察核磁共振信号的方法，它特别适用于测量横向弛豫时间，谱线的自然线宽是由自旋自旋相互作用决定的，但在许多情况下，由于外磁场不够均匀，谱线就变宽了，与这个宽度相对应的横向弛豫时间是前面讨论过的表观横向弛豫时间，而不是了，但用自旋回波法仍可以测出横向弛豫时间。图 7 实际应用中，常用两个或多个射频脉冲组成脉冲序列，周期性的作用于核磁矩系统。比如在射频脉冲作用后，经过时间再施加一个射频脉冲，便组成一个脉冲序列，这些脉冲序列的脉宽和脉距应满足下列条件： （14）脉冲序列的作用结果如图5所示，在射频脉冲后即观察到FID信号；在射频脉冲后面对应于初始时刻的处可以观察到一个“回波”信号。这种回波信号是在脉冲序列作用下核自旋系统的运动引起的，所以称为自旋回波。图 8自旋回波矢量图解以下用图8来说明自旋回波的产生过程。图8（a）表示体磁化强度在射频脉冲作用下绕轴转到轴上；图8（b）表示脉冲消失后核磁矩自由进动受到不均匀的影响，样品中部分磁矩的进动频率不同，引起磁矩的进动频率不同，使磁矩相位分散并呈扇形展开。为此可把看成是许多分量之和。从旋转坐标系看来，进动频率等于的分量相对静止，大于的分量（图中以代表）向前转动，小于的分量（图中以为代表）向后转动；图8（c）表示射频脉冲的作用使磁化强度各分量绕轴翻转，并继续它们原来的转动方向运动；图8（d）表示时刻各磁化强度分量刚好汇聚到轴上；图8（e）表示以后，用于磁化强度各矢量继续转动而又呈扇形展开。因此，在处得到如图7所示的自旋回波信号。由此可知，自旋回波与FID信号密切相关，如果不存在横向弛豫，则自旋回波幅值应与初始的FID信号一样，但在时间内横向弛豫作用不能忽略，体磁化强度各横向分量相应减小，使得自旋回波信号幅值小于FID信号的初始幅值，而且脉距越大则自旋回波幅值越小，并且回波幅值与脉距存在以下关系： （15）式（15）中，是射频脉冲刚结束时FID信号的初始幅值，实验中只要改变脉距，则回波的峰值就相应的改变，若依次增大测出若干个相应的回波峰值，便得到指数衰减的包络线。对(15)式两边取对数，可以得到直线方程 （16）上式中作为自变量，则直线斜率的倒数便是。这就是我们测的方法。二 程序设计基于上面所述的三个过程我们进行了编程可视化现逐一介绍如下：1单个核的拉摩进动运行结果： 有两条黄色箭头代表y与z方向的磁场。红色箭头为总磁场，而绿色的场线束代表旋转着的磁矩。右侧的三个下拉菜单分别控制：面板上设有三个参数，可以进行调解控制动画的变化。拖动条控制视角。设有停止与退出按钮。点击plot按钮后可以观察到在实际情况中单条磁矩的拉莫尔进动。另外值得一提的是此面板的设计是完全靠代码实现的，没有依靠GUI工具箱，相对比较简洁。2弛豫过程（自由感应衰减（FID）信号）此面板依旧是采用程序代码设计。分为三个绘图区，分别为：磁矩尖端的实际运动轨迹（左上）；磁矩大小随时间的衰减（右上）；与所测到的自感衰减信号随时间变化的情况（左下）。 点击“绘图”按钮后，会随着时间三幅图同时运作，展示出对应的变化关系。 3自旋回波法测量横向弛豫时间（脉冲序列方式）界面简洁，功能全面。左侧为运动模拟区：描绘若干条磁矩在旋转系中的运行行为。运动模拟区下方为值改变条，通过拉动其来改变实验的值。右侧上方为信号绘图区，为在动系y轴的方向设置探测器来接收辐射的电测波信号，横轴为时间；纵轴为所测得的信号幅值。信号绘图区下方为数据记录区，分为：序号（记录改变值，并运行次数。）；脉冲长度（读取通过拉动条改变实验的值）；回波强度最大值（即公式（15）中的值，这样只要测出不同值下的值，利用公式（16）就能够通过线性拟合，求出的值。完全按照实验的操作方式。另外附加了一个功能，即遇到相同的值时，不记录实验数据）在数据记录区下面是状态区，其作用是提示用户当前程序的运行状态，指导操作。在面板的最下部设计四个控制按钮：运行：取拉动条所取的值进行模拟；数据分析：打开数据分析面板，对所得的实验数据，经行数据拟合等功能；清零：将所有实验数据清零；退出：关闭此面板。上面是数据处理面板。试将前面所的到的实验数据经行拟合，从而得到的值。左侧为绘图区；右侧为结果显示区，即将所得到的拟合后公式显示出来，并展示所求得的的值。下面为具体的操作实例：三 总结与感想 本实验模拟能够基本上满足实验教学的需要，且图像清晰，形象直观，操作简便。可以作为近代物理实验教学时的课件，或者给同学们提供实验预习之用。但由于需要将方程完全解出，才能一次性将图画出，所以解方程的时间较长是本可视化的一大缺憾可以在以后的学习中继续体会，试图找出能够缩短运行时间的办法。这次做这个实验的可视化可以说是十分艰难，我们从开始对近代物理实验一无所知到题目的完成，很大程度上是老师耐心指导的结果。要做好一个物理实验可视化，首先要对这个实验有相当程度上的了解，而且在图像和原理上不能含糊，在这一点上尤其感谢廖红波老师的耐心指导和指正，让我们的物理概念得到了澄清，物理图像更加清晰。感谢赵虎老师给我们的这么有意思的题目。使我们在做的过程中学到了太多课程上所不能学到的东西，希望在未来的日子里还能有机会做类似的可视化，在提高我们才能的同时，为我们系的实验教学做自己的一点点贡献。 参考文献：1FD-PNMR-型脉冲核磁共振实验仪实验指导书2 FD-PNMR-II脉冲核磁共振仪说明书3近代物理实验.熊俊.北京师范大学出版社.2007年8月第一版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);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,定磁场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,text, 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,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););%调用函数%*%核磁共振被调用函数%*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=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,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 m=1 delete(h0,h1,h2,h3,h4,h5,h6); break; %设置停止开关 endfor i=N+1:length(T) if m=1 break; %设置停止开关 end a1=-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*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(Position,80,80,850,600,. Name,核磁共振弛豫过程,. Color,.8 0.9 .9,. NumberTitle,off);ha=axes(Position,0.04 0.38 0.45 0.55,. %设置画图区1 Color,.8 0.9 .9);hb=axes(Position,0.53 0.38 0.45 0.55,. %设置画图区2 Color,.8 0.9 .9);hc=axes(Position,0.1 0.05 0.55 0.25,. %设置画图区3 Color,.8 0.9 .9);hdo=uicontrol(hf,. %设置画图按钮 Style,pushbutton,. 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); %求解微分方程得弛豫过程磁矩的变化数据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),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 off set(hha,EraseMode,none); set(hha2,EraseMode,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&Y(i+1,1)0; %判断是否有残留图像，若有则删除 delete(childrenhandles);end;if mark=0 %mark=0表示脉冲长度有变化，重新计算数据 tao=get(handles.slider1,value);%将滑动条的值赋值给变量tao kk=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.xuhao,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)-2axes(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 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(:,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)*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=tt;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,ist