aUCBLogo Demos and Tests / nmr_qm_sim
to nmr_qm_sim
xy
end
be XYGraph frm width height ~
xmin xmax ymin ymax ticksx ticksy mode
local [g mx my mw mh x0 y0 w0 h0 zx zy zx0 zy0 p0 p pz
mouseleft mouseright inmotion c
xa ya xxa yya
xmin0 xmax0 ymin0 ymax0
maxx maxy]
g=(Graph frm
wxdefault_frame_style+wxfull_repaint_on_resize ;+wxStay_on_Top
[0 0] List width height [Graph])
;g=[]
norefresh
setUpdateGraph false
hideTurtle
mx=0
my=0
x0=0
y0=0
mw=0
mh=0
w0=0
h0=0
zx=1
zy=1
zx0=1
zy0=1
c=0
p0=[0 0 0]
p=[0 0 0]
pz=[0 0 0]
mouseleft=false
mouseright=false
inmotion=false
yminh=ymax
ymax=ymin
ymin=yminh
xminh=xmax
xmax=xmin
xmin=xminh
xmin0=xmin
xmax0=xmax
ymin0=ymin
ymax0=ymax
maxx=width-70
maxy=height-50
be init
WindowOnLeftDown g
[
if not mouseleft
[ GraphSetCurrent g
p0=MousePos+p
x0=p0.1+mx
y0=p0.2+my
mouseleft=true
]
]
WindowOnLeftUp g
[
mouseleft=false
GraphSetCurrent g
p0=p0-MousePos+p
mx=p.1
my=p.2
]
WindowOnRightDown g
[
if not mouseright
[ GraphSetCurrent g
p0=MousePos+pz
w0=p0.1+mw
h0=p0.2+mh
mouseright=true
]
]
WindowOnRightUp g
[
mouseright=false
GraphSetCurrent g
p0=p0-MousePos+pz
mw=p.1
mh=p.2
zx0=zx0*(1+mw/800)
zy0=zy0*(1+mh/600)
]
;comment [
WindowOnMotion g
[
if mouseleft and2 not inmotion
[ GraphSetCurrent g
p=p0-MousePos
mx=p.1
my=p.2
xmin=(xmin0+xmax0)/2-(xmax0-xmin0)/2*zx-mx/maxx*(xmax0-xmin0)
xmax=(xmin0+xmax0)/2+(xmax0-xmin0)/2*zx-mx/maxx*(xmax0-xmin0)
ymin=(ymin0+ymax0)/2-(ymax0-ymin0)/2*zy-my/maxy*(ymax0-ymin0)
ymax=(ymin0+ymax0)/2+(ymax0-ymin0)/2*zy-my/maxy*(ymax0-ymin0)
clearScreen
drawplot
updateGraph
]
if mouseright and2 not inmotion
[ GraphSetCurrent g
pz=p0-MousePos
mw=pz.1
mh=pz.2
zx=zx0*(1+mw/800)
zy=zy0*(1+mh/600)
xmin=(xmin0+xmax0)/2-(xmax0-xmin0)/2*zx-mx/maxx*(xmax0-xmin0)
xmax=(xmin0+xmax0)/2+(xmax0-xmin0)/2*zx-mx/maxx*(xmax0-xmin0)
ymin=(ymin0+ymax0)/2-(ymax0-ymin0)/2*zy-my/maxy*(ymax0-ymin0)
ymax=(ymin0+ymax0)/2+(ymax0-ymin0)/2*zy-my/maxy*(ymax0-ymin0)
clearScreen
drawplot
updateGraph
]
]
;]
end
be plot xa_ ya_ c_
xa=xa_
ya=ya_
c=c_
drawplot
end
be drawplot
local [x y y2 i h w hq wq hx hy rx ry]
WindowMode
setPenSize [.5 .5]
setPC 0
PenUp
setXY -maxx -maxy
PenDown
setY maxy
setX maxx
setY -maxy
setX -maxx
tick=10
txtx=25
setH 90
setLabelSize [25 25]
setLabelAlign 0 0
setPrintPrecision 4
setLabelFont "Times
h=ymax-ymin
hq=(ymax0-ymin0)*2^(round -0.5+(ln abs h/(abs ymax0-ymin0))/ln 2)
w=xmax-xmin
wq=(xmax0-xmin0)*2^(int -0.5+(ln abs w/(abs xmax0-xmin0))/ln 2)
for [y -ticksy ticksy]
[
hy=-hq*(y/ticksy-0.5)
hy=hy+(remainder -(ymax+ymin)/2 hq/ticksy)
ry=hy+(ymax+ymin)/2
hy=-hy*maxy*2/h
if hy > maxy or2 hy < -maxy [continueLoop]
;(pr maxy h hy)
;if this == xy::g3 [pr hy]
setY hy
PenDown
setX -maxx+tick
setX -maxx
PenUp
setX maxx-tick
PenDown
setX maxx
PenUp
setX -maxx-txtx
Label ry
;pr ry
setX -maxx
]
setY -maxy
for [x -ticksx 2*ticksx]
[ hx=-wq*(x/ticksx-0.5)
hx=hx+(remainder -(xmax+xmin)/2 wq/ticksx)
rx=hx+(xmax+xmin)/2
hx=-hx*maxx*2/w
if hx > maxx or2 hx < -maxx [continueLoop]
setX hx
PenDown
setY -maxy+tick
setY -maxy
PenUp
setY maxy-tick
PenDown
setY maxy
PenUp
setY -maxy-txtx
label rx
setY -maxy
]
setPenSize [0 0]
x=(xmax+xmin)/2
y=(ymax+ymin)/2
fx=2/(xmax-xmin)*maxx
fy=2/(ymax-ymin)*maxy
xxa=saturateBelow -maxx saturateAbove maxx (-xa+x)*fx
yya=saturateBelow -maxy saturateAbove maxy (-ya+y)*fy
ifelse mode == "pixels
[ setPixelXY xxa yya c
][
setPC c
PenUp
; x=maxx
; repeat count xxa
; [ j=repcount
; if xxa.j < x [x=xxa.j y=yya.j]
; ]
; setXY x y
setXY xxa.1 yya.1
PenDown
setXY xxa yya
PenUp
]
updateGraph
;pause
end
be clean
GraphSetCurrent g
clearScreen
end
end
be delta i j
output ifelse i==j [1][0]
end
to xy
intensity=FloatArray 1
for [p 0.01 1 0.01]
[ sim p
]
end
to sim param1
MSamples=500
IMAX=1000
tscale=10e-6
MSamples=IMAX/tscale
t=rSeqFA 0 tscale IMAX
signalOn=true
frequency=4.2575e6 ;in MHz
gamma=42.575e6 ;MHz/Tesla
B0=param1 ;0.11744 ;1.1744
B1=0.2 ;0.05 ;0.2 ;0.1
dBmax=0.1 ;0.01 ;0.2 ;0.01
tau=2.5 ;10 ;2.4 ;5
tau2=tau
tau3=tau
tauabs=[80] ;[25 50 100 200 500 1000 2000 5000]
taubc=50
T1=100 ;3
T2=1 ;0.03 ;0.07 ;0.1
T2e=T2
spins=1000 ;10
filterpara=300
ww=0 ;0.001
expo=1 ;1.4142 ;1.33 ;1 ;1.5 ;1.333
B=Array 3
omega=360*frequency
t0=1e-6
ta=t0+1e-6
pulses=(t>t0)*(t<ta)
B.1=(sin t*omega)*B1*pulses
B.2=;rSeqFA 0 0 IMAX
(cos t*omega)*B1*pulses
B.3=rSeqFA B0 B0 IMAX
I=1/2
I_=Array 3
m=(FloatArray 2*I+1 0)
Iplus=(Array 2*I+1 0)
Iminus=(Array 2*I+1 0)
for [j 0 2*I]
[ m.j=-I+j
Iplus.j=(ComplexArray 2*I+1 0)
Iminus.j=(ComplexArray 2*I+1 0)
]
for [j 0 2*I]
[ for [k 0 2*I]
[ Iplus.j.k=(sqrt I*(I+1)-m.j*(m.j+1))*(delta 2*(m.j+1) 2*m.k)
Iminus.j.k=(sqrt I*(I+1)-m.j*(m.j-1))*(delta 2*(m.j-1) 2*m.k)
]
]
I_.1=(Iplus+Iminus)/1
I_.2=-(Iplus-Iminus)/1i
I_.3=-((MatrixXMatrix I_.1 I_.2)-(MatrixXMatrix I_.2 I_.1))/2i
dt=tscale/IMAX
h_bar=6.62620e-34/(2*pi)
qcc=2 ;in MHz
theta=0
phi=0
eta=0
S=I_*h_bar
Hz=S*dt/(1i*h_bar)*gamma
ifelse I > 1/2
[ Hq=((MatrixXMatrix I_.3 I_.3)*3-I+(I+1))
*dt/(1i*4*I*(2*i-1))*qcc/2
*(3*(sqr cos theta)-1-eta*(sqr sin theta)*(cos 2*phi))
][ Hq=(Array 2 0)
Hq.0=(ComplexArray 2 0)
Hq.1=(ComplexArray 2 0)
; Hq.(1).0=0.001
; Hq.(1).1=-0.001
]
chi0=(ComplexArray 2*I+1 0)
chi0.0=1
;stop
showAll=false
running=true
if not Name? "frm
[
frm=Frame [][MyFrame]
wxresize_border+wxcaption+wxsystem_menu+wxclose_box
+wxfull_repaint_on_resize
[0 0][400 300]
ifelse showAll
[ g=XYGraph frm 900 300 0 tscale -1 1 4 4 "lines
setScreenRange -900 -300 900 300
][ g=XYGraph frm 400 300 0 tscale -1.5 1.5 4 4 "lines
setScreenRange -400 -300 400 300
g2=XYGraph frm 400 300 0 tscale -1 1 4 4 "lines
setScreenRange -400 -300 400 300
g3=XYGraph frm 400 300 0 100 -1 1 4 4 "lines
setScreenRange -400 -300 400 300
]
g'init
g2'init
g3'init
if showAll
[ g2=XYGraph frm 900 300 0 tscale -1.2 1.2 4 4 "lines
setScreenRange -900 -300 900 300
g2'init
g3=XYGraph frm 900 300 0 tscale -2 2 4 4 "lines
setScreenRange -900 -300 900 300
g3'init
]
; txtRe1=StaticText frm [ContrastRe]
bs=BoxSizer wxvertical
; bsh=BoxSizer wxHORIZONTAL
; bsh2=BoxSizer wxHORIZONTAL
; BoxSizerAdd bs bsh 100 wxExpand+wxALIGN_CENTER 0
; BoxSizerAdd bs bsh2 100 wxExpand+wxALIGN_CENTER 0
; bsv=BoxSizer wxVertical
; bsv2=BoxSizer wxVertical
; BoxSizerAdd bsh bsv 100 wxExpand+wxALIGN_CENTER 0
; BoxSizerAdd bsh bsv2 100 wxExpand+wxALIGN_CENTER 0
BoxSizerAdd bs g'g 100 wxExpand+wxALIGN_CENTER 0
BoxSizerAdd bs g2'g 100 wxExpand+wxALIGN_CENTER 0
BoxSizerAdd bs g3'g 100 wxExpand+wxALIGN_CENTER 0
if showAll
[ BoxSizerAdd bs g3'g 150 wxExpand+wxALIGN_CENTER 0
]
; bst=BoxSizer wxHORIZONTAL
; BoxSizerAdd bs bst 20 wxExpand+wxALIGN_CENTER 0
; bsRe =BoxSizer wxVertical
; bsAbs=BoxSizer wxVertical
; BoxSizerAdd bst bsRe 20 wxExpand+wxALIGN_CENTER 0
; BoxSizerAdd bst bsAbs 20 wxExpand+wxALIGN_CENTER 0
; BoxSizerAdd bsRe txtRe1 100 wxExpand+wxALIGN_CENTER 0
; BoxSizerAdd bsRe txtRe2 100 wxExpand+wxALIGN_CENTER 0
; BoxSizerAdd bsRe txtRe3 100 wxExpand+wxALIGN_CENTER 0
; BoxSizerAdd bsAbs txtAbs1 100 wxExpand+wxALIGN_CENTER 0
; BoxSizerAdd bsAbs txtAbs2 100 wxExpand+wxALIGN_CENTER 0
; BoxSizerAdd bsAbs txtAbs3 100 wxExpand+wxALIGN_CENTER 0
FrameSetSizer frm bs
ifelse showAll
[ FrameSetClientSize frm 1000 950
][ FrameSetClientSize frm 640 960
]
running=true
stopping=false
FrameOnClose frm [running=false]
FrameOnChar frm [running=false]
GraphOnChar g'g [running=false]
if showAll
[ GraphOnChar g2'g [running=false]
GraphOnChar g3'g [running=false]
]
]
;]
setPrintPrecision 3
fn="nmr_qm_sim.html
closeall
openWrite fn
setWriter fn
type "<html>
<body>
while [not key? and2 running]
[ if empty? tauabs [break]
tauab=first tauabs
tauabs=butfirst tauabs
j=repcount
signalOn=true
Mx=(FloatArray IMAX 0)
My=(FloatArray IMAX 0)
Mz=(FloatArray IMAX 0)
chi=(NMR_QMsim Hz Hq chi0 B S Mx My Mz)
ii=int 2*I+1
chir=Array ii
chir.1=real chi.1*conjugate chi.1
chir.ii=real chi.ii*conjugate chi.ii
comment [
chit=transposeMatrix chi
repeat IMAX
[ i=repcount-1
Mx.i=(real 0.0+chit.i*(MatrixXVector S.1 chit.i))/h_bar
My.i=(real 0.0+chit.i*(MatrixXVector S.2 chit.i))/h_bar
Mz.i=(real 0.0+chit.i*(MatrixXVector S.3 chit.i))/h_bar
]
]
i=(0+sqr items int ta*IMAX/tscale IMAX-1 Mx)
intensity=lput i intensity
intensityn=intensity/(max intensity)
xI=rseqfa 0 100 (count intensity)
GraphSetCurrent g'g
clearScreen
(g'plot t chir.ii HSBA 320 1 1 1)
(g'plot t chir.1 HSBA 240 1 1 1)
setpc 0 pu setxy 0 275 pd label word [Up and Down levels, p=] param1 pu
updateGraph
setSaveSize list g'width*2 g'height*2
; saveScreen (word "nmr_qm_sim_chi.png)
GraphSetCurrent g2'g
clearScreen
(g2'plot t My HSBA 320 1 1 1)
(g2'plot t Mx HSBA 240 1 1 1)
; (g2'plot t Mz HSBA 0 0 0.5 1)
setpc 0 pu setxy 0 275 pd label [Mx, My, Mz] pu
updateGraph
setSaveSize list g2'width*2 g2'height*2
; saveScreen (word "nmr_qm_sim_M.png)
GraphSetCurrent g3'g
clearScreen
(g3'plot xI intensityn HSBA 320 1 1 1)
setpc 0 pu setxy 0 275 pd label [Spectrum] pu
updateGraph
setSaveSize list g3'width*2 g3'height*2
; saveScreen (word "nmr_qm_sim_M.png)
if key? or2 not running [break]
if showAll
[ signalOn=false
(NMR_QMsim Hz Hq chi0 B)
Mx0=FloatArray Mx
GraphSetCurrent g2'g
clearScreen
(g2'plot t Mx HSBA 240 1 1 1)
(g2'plot t My HSBA 320 1 1 1)
(g2'plot t P0 HSBA 0 0 0.5 1)
Uqre=lowPassFilter2 Mx0*(sin t*360*frequency)/4 filterpara
Uqim=lowPassFilter2 Mx0*(cos t*360*frequency)/4 filterpara
Uq =sqrt (sqr Uqre )+(sqr Uqim )
(g'plot t Uq HSBA 0 1 1 0.1)
(g'plot t Uqre HSBA 240 1 1 0.1)
(g'plot t Uqim HSBA 120 1 1 0.1)
setpc 0 pu setxy 0 275 pd label word [A off, tauAB=] tauab pu
updateGraph
setSaveSize list g2'width*2 g2'height*2
saveScreen (word "blochequations_off_ tauab "us.png)
if key? or2 not running [break]
GraphSetCurrent g3'g
clearScreen
Mxd=Mx1-Mx0
(g3'plot t Mxd HSBA 120 1 1 1)
(g3'plot t P1 HSBA 0 0 0.5 1)
Uqre=lowPassFilter2 Mxd*(sin t*360*frequency)/6 filterpara
Uqim=lowPassFilter2 Mxd*(cos t*360*frequency)/6 filterpara
Uq =sqrt (sqr Uqre )+(sqr Uqim )
(g'plot t Uq HSBA 0 1 1 0.1)
(g'plot t Uqre HSBA 240 1 1 0.1)
(g'plot t Uqim HSBA 120 1 1 0.1)
setpc 0 pu setxy 0 275 pd label word [on-off, tauAB=] tauab pu
updateGraph
setSaveSize list g3'width*2 g3'height*2
saveScreen (word "blochequations_diff_ tauab "us.png)
]
GC
(type "
<img src="blochequations_on_ tauab "us.png" alt="on tauab "us" width="100%"></br>
<img src="blochequations_off_ tauab "us.png" alt="off tauab "us" width="100%"></br>
<img src="blochequations_diff_ tauab "us.png" alt="diff tauab "us" width="100%"></br>)
]
type "
</body>
</html>
setWriter []
close fn
end