aUCBLogo Demos and Tests / nmr_carterspectra

setCaseIgnored false

to nmr_carterspectra
   setUpdateGraph "false
   hideTurtle
   clearScreen
   WindowMode
   setpensize 0
;   disableTextureFont
   enablePolySmooth
   ch=0
nmr
stop
   forever
   [
      clearScreen
      ;pr repcount
      catch "error [ignore runResult [nmr]]
      err=Error
      if :err != []
      [   clearText
         print :err.2
         waitMS 200
      ]
      updategraph
      updateVars
      dispatchMessages
      waitMS 100
      if ch==Char 27 [stop]
   ]
end

to nmr
   ni=1000
   x=rSeqFA 400 -400 :ni

   spectrum=readfile "../La2RuO5/56Mhz/0_700K_ausw.dat
   specx=:spectrum.1
   specy=:spectrum.2
   specy=:specy/(max :specy)
   specn=count :specy
   specy_cor=:specy-rSeqFA 0.0 0.05 specn
   specy=:specy_cor*400-200

;   B=specx
;   nni=count specx
;   BI=Array nni
;   repeat nni
;   [   i=repcount
;      BI.i=list B.i specy_cor.i
;   ]
;   foreach BI [pr ?]

   spectrum2=readfile "../La2RuO5/64Mhz/1K_2_ausw.dat
   specx2=:spectrum2.1/64*56
   specy2=:spectrum2.2
   specy2=:specy2/(max :specy2)*400-200

;   B=specx2
;   nni=count specx2
;   BI=Array nni
;   repeat nni
;   [   i=repcount
;      BI.i=list B.i spectrum2.(2).i
;   ]
;   foreach BI [pr ?]

   for [p 0 0]
   [   intensity=FloatArray :ni

      I=7/2
      MHz=1e6
      Tesla=1
      gamma=6.061*:MHz/:Tesla
      markerx=56.25*:MHz
      nu0=56.25*:MHz
      nuref=56*:MHz
      nuQ=4.8*:MHz
      eta=0.85
      expo=0
      K1=0 ;0.0045
      K2=0
      ntheta=200
      nphi=200
      I0=1
      Inu0=:nu0
      nnu=3*:Tesla*:gamma
      smoothwidth=0.225*:MHz
      compwidth=10*:MHz
      maxy=200
      maxx=400
      ticksx=5
      tick=15
      txtx=25

      xmax=(-:nu0+:nnu/2)/:gamma
      xmin=(-:nu0-:nnu/2)/:gamma
      w=:xmax-:xmin
      wq=(:xmax-:xmin)*2^(int -0.5+(ln :w/(:xmax-:xmin))/ln 2)
      markerx=:markerx/:gamma

      nmrloop
      if ch==char 27 [break]
   ]
;   B=rSeqFA -xmin -xmax ni
;   BI=Array ni
;   repeat ni
;   [   i=repcount
;      BI.i=list B.i smoothedIntensity.i
;   ]
;   foreach BI [pr ?]
end

to nmrloop
   repeat 1
   [
;      computeIntensityRatios
;stop
      intensity=NMRCarterIntensity
         :I :nu0 :nuQ :nuref
         :eta :K1 :K2
         :ntheta :nphi :Inu0 :nnu
         :intensity :expo
;      smoothedIntensity=:intensity
      smoothedIntensity=GaussianSmooth
         :nnu :smoothwidth :compwidth :intensity
      mi=(max :smoothedIntensity)
      if :mi != 0
      [   intensityScaled=:smoothedIntensity/:mi*I0*400-maxy
      ]
      hideTurtle
      clearScreen
      WindowMode
      setpensize 0
      penup
      setXY -400 :intensityScaled.1
      pendown
      setXY :x :intensityScaled

      specxsc=(:specx+:xmax)/w*maxx*2-maxx
      setPC HSB 120 1 1
      penup
      setXY :specxsc.1 :specy.1
      pendown
      setXY :specxsc :specy
      penup
      setPC 0

      specxsc2=(:specx2+:xmax)/w*maxx*2-maxx
      setPC HSB 0 1 1
      penup
      setXY :specxsc2.1 :specy2.1
      pendown
      setXY :specxsc2 :specy2
      penup
      setPC 0

      markerxsc=(:markerx+:xmax)/w*maxx*2-maxx
      setPC HSB 240 1 1
      penup
      setXY :markerxsc -:maxy
      pendown
      setXY :markerxsc :maxy
      penup
      setPC 0

      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
         seth 90
         label :rx
         setY -:maxy
      ]

      setXY 0 -250
      Label "B (T)

      setXY 0 -280
      Label :p

      updateGraph
      if Key?
      [   ch=readChar
         if :ch==Char 27 [break]
         if :ch=="s
         [   penup
            setpensize 0
         ;   boundingbox
            catch "error
            [   saveScreenVector "nmr_carterspectrum.pdf
            ]
            err=error
            if :err != []
            [   print :err.2
            ]
         ]
      ]
      (GC "true)
   ]
end

to computeIntensityRatios
   r=(FloatArray 2*:I -:I+1/2)
   j=1
   ratios={{1}@0 {4 3}@0 {9 8 5}@0 {16 15 12 7}@0 {25 24 12 16 9}@0}
   for [m -:I+1 :I]
   [   j=Int :m-1/2
      factor=:ratios.(Int :I+1/2).abs :j
      intensity=FloatArray :ni
      intensity=NMRCarterIntensity_m
         :I :nu0 :nuQ :nuref
         :eta :K1 :K2
         :ntheta :nphi :Inu0 :nnu
         :intensity :m
      smoothedIntensity=GaussianSmooth :nnu :smoothwidth :intensity
      mi=(max :smoothedIntensity)
      if :mi != 0
      [   intensityScaled=:smoothedIntensity/:mi*factor*40-200
      ]
      penup
      setpc HSB 180+180*:m/:I 1 1
      setXY -400 :intensityScaled.1
      pendown
      setXY :x :intensityScaled
      setitem :j :r 0+:intensity
   ]
;      show :r
end

to readfile f
   local [x y d]
   x=[]
   y=[]
   openRead :f
   setReader :f
   while [not eof?]
   [   d=readlist
      x=fput :d.1 :x
      y=fput :d.2 :y
   ]
   setReader []
   close :f
   output list reverse FloatArray :x reverse FloatArray :y
end