aUCBLogo Demos and Tests / blochequations6
be blochequations6
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 xy
video=false
video=true
if video
[ (VideoStart "blochequations6_vartau_200us_12.divx 20)
]
for [i 0 20 0.05]
[ bloch i
if video
[ VideoFrame
]
]
if video
[ VideoEnd
]
end
be bloch param1
MSamples=500
IMAX=100000
tscale=IMAX/MSamples
t=rSeqFA 0 tscale IMAX
Mx=FloatArray IMAX
My=FloatArray IMAX
P0=FloatArray IMAX
P1=FloatArray IMAX
signalOn=true
frequency=1 ;in MHz
B1c=0.1732 ;0.05 ;0.2 ;0.1
dBmax=0.1 ;0.01 ;0.2 ;0.01
tau=param1 ;2.4 ;5
tauabs=[20];[25 50 100 200 500 1000 2000 5000]
T1=100
T2=1 ;0.03 ;0.07 ;0.1
T2e=T2
spins=100 ;10
filterpara=300
ww=0 ;0.001
expo=1 ;1.4142 ;1.33 ;1 ;1.5 ;1.333
Mx1=FloatArray IMAX
Mx0=FloatArray IMAX
n1=0
n0=0
running=true
showAll=false
; showAll=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.2 1.2 4 4 "lines
setScreenRange -900 -300 900 300
][ g=XYGraph frm 400 300 0 tscale -1.2 1.2 4 4 "lines
setScreenRange -400 -300 400 300
]
g'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
]
bs=BoxSizer wxVertical
BoxSizerAdd bs g'g 100 wxExpand+wxALIGN_CENTER 0
if ShowAll
[ BoxSizerAdd bs g2'g 100 wxExpand+wxALIGN_CENTER 0
BoxSizerAdd bs g3'g 150 wxExpand+wxALIGN_CENTER 0
]
FrameSetSizer frm bs
ifelse showAll
[ FrameSetClientSize frm 640 950
][ FrameSetClientSize frm 640 480
]
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="blochequations5.html
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
(NMRBlochSimSpectrum
Mx My P1
signalOn frequency MSamples B1c dBmax
tau tauab spins
T1 T2 T2e ww expo)
Mx1=FloatArray Mx
GraphSetCurrent g'g
clearScreen
(g'plot t Mx HSBA 240 1 1 1)
(g'plot t My HSBA 320 1 1 1)
(g'plot t P1 HSBA 0 0 0.5 1)
Uqre=lowPassFilter2 Mx1*(sin t*360*frequency)/4 filterpara
Uqim=lowPassFilter2 Mx1*(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 [tau=] tau pu
updateGraph
setSaveSize list g'width*2 g'height*2
saveScreen (word "blochequations_on_ tauab "us.png)
if key? or2 not running [break]
if showAll
[
signalOn=false
(NMRBlochSimSpectrum
Mx My P0
signalOn frequency MSamples B1c dBmax
tau tauab spins
T1 T2 T2e ww expo)
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