aUCBLogo Demos and Tests / nmr_carterspectra_f_HgCrS

setCaseIgnored false

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

to zzz
   ni=1000
   x=rSeqFA -400 400 :ni
;comment [
   spectrum=readfile "../HgCrS/46Mhz/2_700K_149.fsp
   specx=:spectrum.1
   specx=(-specx+46.66)*1000000
   specy=:spectrum.2
   specy=:specy/(max :specy)
   specn=count specy
   specy=:specy*400-200

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

   for [p 0 0]
   [
      intensity=FloatArray :ni
      I=3/2
      kHz=1e3
      MHz=1e6
      markerx=46*:MHz
      nu0=46.68*MHz
      nuref=nu0
      nuQ=900*:kHz
      eta=0.07
      expo=0.5
      K1=0.0005
      K2=0.0005
      ntheta=200
      nphi=200
      I0=1
      Inu0=:nu0
      nnu=4*MHz
      smoothwidth=0.025*:MHz
      compwidth=1*:MHz
      nmr
      dnu=0.11*MHz
      I0=1      Inu0=:nu0-dnu      nmr
      I0=0.3   Inu0=:nu0-2*dnu   nmr
      I0=0.3   Inu0=:nu0+2*dnu   nmr
      I0=0.05   Inu0=:nu0+3*dnu   nmr
      I0=0.2   Inu0=:nu0+4*dnu   nmr
      I0=0.05   Inu0=:nu0+7*dnu   nmr

      I0=1
      maxy=200
      maxx=400
      ticksx=8
      tick=15
      txtx=25
      filename="nmr_carterspectrum_HgCrS.pdf

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

      if ch==char 27 [break]
   ]
end

to nmr
   intens=FloatArray :ni
;comment[
   intens=NMRCarterIntensity
      :I :nu0 :nuQ :nuref
      :eta :K1 :K2
      :ntheta :nphi :Inu0 :nnu
      :intens :expo
   intensity=intensity+intens*I0
;]
comment [
   for [theta 0 90 2]
   [   intens=NMRCarterIntensity_oneTheta
      :I :nu0 :nuQ :nuref
      :eta :K1 :K2
      :theta :nphi :Inu0 :nnu
      :intens :expo
      factor=abs (sin theta)
      if factor != 0
      [   factor=factor/factor^expo
      ]
      intensity=intensity+intens*I0*factor
   ]
]
end

to draw
   repeat 1
   [
;      computeIntensityRatios
;stop
      smoothedIntensity=:intensity
;comment[
      smoothedIntensity=GaussianSmooth
         :nnu :smoothwidth :compwidth :intensity

      mi=(max :smoothedIntensity)
      ifelse :mi != 0
      [   intensityScaled=(:smoothedIntensity/:mi)*I0*400-maxy
      ][   intensityScaled=:intensity
      ]
      hideTurtle
      clearScreen
      WindowMode
      setpensize 0
      penup
      setXY -400 :intensityScaled.1
      pendown
      setXY :x :intensityScaled
;comment [
;      specxsc=(:specx*:gamma+:xmax)/w*maxx*(-2)+maxx
      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
;]
comment [
      specxsc2=(:specx2*:gamma+:xmax)/w*maxx*(-2)+maxx
      setPC HSB 0 1 1
      penup
      setXY :specxsc2.1 :specy2.1
      pendown
      setXY :specxsc2 :specy2
      penup
      setPC 0
]
      markerxsc=maxx-(:markerx+:xmax)/w*maxx*2
      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/1e6
         setY -:maxy
      ]

      setXY 0 -250
      Label "f (MHz)

      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 filename
            ]
            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
   firstline=readList
   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