be blochequations2 n=100 M0=1 M=Array n dM_dt=Array n repeat n [ M.repcount=FloatArray {0 1 0}*M0 dM_dt.repcount=FloatArray {0 0 0} ] B=FloatArray {0 1 0} B1w=FloatArray {0 1 0} dB=Array n dBmax=0.1 repeat n [ j=repcount dB.j=FloatArray {0 0 0} dB.j.2=(sqr sqr rnd)*dBmax if rnd < 0.5 [dB.j.2=-dB.j.2] ] Tesla=1 B0=1*Tesla B1c=0.2*Tesla B1=B1c Mo=[] cMo=[] sec=1 us=1e-6*sec MHz=1/us gamma=42*MHz/Tesla T1=1*sec T2=0.01*sec T2e=T2 t=0 dt=0.01/gamma/B0 omega=gamma*B0*60 size=350 it=0 cM=IntArray n repeat n [ j=repcount cM.j=HSB dB.j.2/dBmax*180 1 1 ] setUpdateGraph false hideTurtle WindowMode perspective enableCylinderLines enableTextureFont setLabelAlign 1 1 rotatescene red=RGB 1 0 0 green=RGB 0 1 0 black=RGB 0 0 0 video=false if video [(VideoStart "blochequations.divx 25)] graphics=true cleanit=true ch=char 255 forever [ repeat 5 [ B.1=B1*sin -omega*t B.2=B0 B.3=0 ;B1*cos -omega*t repeat n [ j=repcount Mj=M.j Bj=B+dB.j dM_dtj=dM_dt.j dM_dtj.1=(gamma*(cross Mj Bj).1- Mj.1 /T2e) dM_dtj.2=(gamma*(cross Mj Bj).2+(M0-Mj.2)/T1) dM_dtj.3=(gamma*(cross Mj Bj).3- Mj.3 /T2) d2M_dt2=dM_dtj-dM_dt.j dM_dt.j=dM_dtj M.j=Mj+dM_dtj*dt+d2M_dt2*dt/2 if M.j > M0 [M.j=M.j/Norm M.j] ] t=t+dt ] Mges=FloatArray 3 repeat n [ j=repcount Mges=Mges+M.j ] Mges=Mges/n Mo=fPut Mges*size Mo cMo=fPut HSB 360*t/dt/50000 1 1 cMo it=it+1 if graphics [ ifelse cleanit [ clearScreen axes B1w.1=B1/B1c*sin -omega*t B1w.2=0 B1w.3=0 ;B1*cos -omega*t setPC green setPenSize 8 setPosXYZ B1w*size setPenSize 0 disableCylinderLines disableLighting repeat n [ j=repcount Home setPC cM.j setPosXYZ M.j*size ] Home setPC red setPenSize 8 enableCylinderLines enableLighting setPosXYZ Mges*size setPenSize 0 disableCylinderLines disableLighting repeat it [ j=repcount setPC cMo.j setPosXYZ Mo.j ] enableLighting PenUp setXYZ 300 0 0 setHeading 90 setPC black Label Word "t= t setXYZ 300 -50 0 Label Word "M= Norm Mges setXYZ 300 -100 0 Label Word "f= omega/60/MHz ][ PenUp setPosXYZ Mo.2 PenDown setPC cMo.1 setPosXYZ Mo.1 ] updateGraph if video [VideoFrame] ] if Key? [ ch=readChar if ch==Char 27 [break] if ch=="0 [B1=0] if ch=="1 [B1=B1c] if ch=="+ [omega=omega*1.001] if ch=="- [omega=omega/1.001] if ch=="g [graphics=not graphics] if ch=="c [cleanit=not cleanit] if ch=="r [ repeat n [ M.repcount=FloatArray {0 1 0}*M0 dM_dt.repcount=FloatArray {0 0 0} ] t=0 clearMo ] if ch==Char wxk_delete [ clearMo ] ] rotatescene::once ch=0 ] if video [VideoEnd] be clearMo clearScreen axes setPenSize 0 disableCylinderLines Mo=[] cMo=[] Mo=fPut Mges*size Mo cMo=fPut HSB 360*t/dt/50000 1 1 cMo it=1 end be axes setPC black setPenSize 5 enableCylinderLines home pd setposxyz [1 0 0]*200 pu setposxyz [1 0 0]*250 seth 90 label "x home seth 90 label "o pd setposxyz [0 1 0]*200 pu setposxyz [0 1 0]*250 seth 90 label "z home pd setposxyz [0 0 1]*200 pu setposxyz [0 0 1]*250 seth 90 label "y home pd end be rotatescene dphi=1 phi=20 dtheta=1 eye=array 3 light=array 3 r=900 l=r ;200 dr=1.1 p={0 0 0} dx={10 0 0} theta=30 center={0 0 0} upvector={0 1 0} upv=upvector ddphi=dphi/3 be once while [or ch==char 255 Key?] [ eye.1=r*(cos theta)*sin phi eye.2=r* sin theta eye.3=r*(cos theta)*cos phi upv.1=(cos 90-theta)*sin phi+180 upv.2= sin 90-theta upv.3=(cos 90-theta)*cos phi+180 setEye eye+center-p center-p upv light.1=l*(cos theta)*sin phi light.2=l* sin theta light.3=l*(cos theta)*cos phi setLightPos light redraw if ch==char 255 [ ch=readCharExt ifElse (BitAnd MouseButtons 16)==16 [ dx=cross eye-center upvector dx=(dx/Norm dx)*10 dy=cross dx eye-center dy=(dy/Norm dy)*10 dz=eye-center dz=(dz/Norm dz)*10 if ch==wxk_right [p=p+dx] if ch==wxk_left [p=p-dx] if ch==wxk_up [p=p+dy] if ch==wxk_down [p=p-dy] if ch==wxk_prior [p=p+dz] if ch==wxk_next [p=p-dz] ][ if ch==wxk_right [phi=phi+dphi] if ch==wxk_left [phi=phi-dphi] if ch==wxk_up [theta=theta+dtheta] if ch==wxk_down [theta=theta-dtheta] if ch==wxk_prior [r=r/dr] if ch==wxk_next [r=r*dr] ] ] ifelse Key? [ch=readChar][ch=0] ] ch=0 end end end