aUCBLogo Demos and Tests / blochequations7


to blochequations7
   
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 0List 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 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 [-ticksy ticksy]
      
[   
         
hy=-hq*(y/ticksy-0.5)
         
hy=hy+(remainder -(ymax+ymin)/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 [-ticksx 2*ticksx]
      
[   hx=-wq*(x/ticksx-0.5)
         
hx=hx+(remainder -(xmax+xmin)/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
   
MSamples=500
   
IMAX=300000
   
tscale=IMAX/MSamples
   
t=rSeqFA tscale IMAX 
   
Mx=FloatArray IMAX
   
My=FloatArray IMAX
   
P0=FloatArray IMAX
   
P1=FloatArray IMAX
   
signalOn=true
   
frequency=;in MHz

   
B1c=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=;0.03 ;0.07 ;0.1
   
T2e=T2
   
spins=1000 ;10
   
filterpara=300
   
ww=;0.001
   
expo=;1.4142 ;1.33 ;1 ;1.5 ;1.333

   
Mx1=FloatArray IMAX
   
Mx0=FloatArray IMAX
   
n1=0
   
n0=0
   
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.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
      
]
;      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'100 wxExpand+wxALIGN_CENTER 0
      
if showAll
      
[   BoxSizerAdd bs g2'100 wxExpand+wxALIGN_CENTER 0
         
BoxSizerAdd bs g3'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 480
      
]
      
running=true
      
stopping=false
      
FrameOnClose frm [running=false]
      
FrameOnChar frm [running=false]
      
GraphOnChar g'[running=false]
      
if showAll
      
[   GraphOnChar g2'[running=false]
         
GraphOnChar g3'[running=false]
      
]
   
]
;]
   
setPrintPrecision 3
   
fn="blochequations4.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
      
(NMRBlochSim2 Mx My P1 
         
signalOn frequency MSamples B1c dBmax 
         
tau tau2 tau3
         
tauab taubc 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)/filterpara
      
Uqim=lowPassFilter2 Mx1*(cos t*360*frequency)/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 pu setxy 0 275 pd label word [A ontauAB=] tauab pu
      
updateGraph
      
setSaveSize list g'width*g'height*2
      
saveScreen (word "blochequations_on_ tauab "us.png)
      
      
if key? or2 not running [break]
      
if showAll
      
[   signalOn=false
         
(NMRBlochSim2 Mx My P0 
            
signalOn frequency MSamples B1c dBmax 
            
tau tau2 tau3
            
tauab taubc 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)/filterpara
         
Uqim=lowPassFilter2 Mx0*(cos t*360*frequency)/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 pu setxy 0 275 pd label word [A offtauAB=] tauab pu
         
updateGraph
         
setSaveSize list g2'width*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)/filterpara
         
Uqim=lowPassFilter2 Mxd*(cos t*360*frequency)/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 pu setxy 0 275 pd label word [on-offtauAB=] tauab pu
         
updateGraph
         
setSaveSize list g3'width*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