aUCBLogo Demos and Tests / blochequations3
be blochequations3
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}
]
Mges=FloatArray 3
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.05*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=0*us
t_1=t_0+1*us
t_2=t_1+5*us
t_3=t_2+t_1/2
t_4=t_3+2*us
t_5=t_4+t_1
ti=(list t_0 t_1 t_2 t_3 t_4 t_5)
t=0
dt=0.1/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
]
compile [compute]
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
[ ifelse and not empty? ti [ti.1 < t and2 t < ti.2]
[ B1=B1c
][ if ti.2 < t and2 not empty? ti
[ B1=0
ti=butFirst butFirst ti
]
]
compute
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 compute
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
; M.j=M.j/Norm M.j
]
t=t+dt
]
Mges=FloatArray 3
repeat n
[ j=repcount
if (Norm M.j) > M0 [M.j=M.j/Norm M.j]
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
end
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