loadct,0,/s
ycolors2
red=gt_ycolors2(4,0,0)
green=gt_ycolors2(0,3,0)
blue=gt_ycolors2(0,0,4)
gray=gt_ycolors2(2,2,2)
violet=gt_ycolors2(3,1,4)

for eps=0,1 do begin

   xs=18
   ys=18
   outfil='larmor_frequency.ps'

   pset, xs=xs, ys=ys, eps=eps, outfil=outfil


   xr=[1e4,1e0]
   xr=[1e0,1e4]

   yr=[1e-2,1e3]

;   xtitle='Distance [Rs]'
   xtitle='Magnetic Field Strength [nT]'
   ytitle='Proton Larmor Radius [Re]'
   ytickname=['0.01','0.1','1.0','10','100','1000']

   plot_oo,xr, yr, yr=yr,ytitle=ytitle, xtitle=xtitle, xr=xr, ytickname=ytickname,$
            xs=9,/nodata

   oplot,[10000,1],[1000,1000]

   clip_xyouts, 0.98,0.95,'1000 MeV', color=violet, al=1
   clip_xyouts, 0.98,0.91,'100 MeV', color=blue, al=1
   clip_xyouts, 0.98,0.87,'30 MeV', color=green, al=1
   clip_xyouts, 0.98,0.83,'10 MeV', color=red, al=1

;   ke=10.^(findgen(301)/300*5-1)
   
   for k=0,3 do begin
      
      case k of
         0: begin
            ke=10.
            color=red
         end
         1: begin
            ke=30.
            color=green
         end
         2: begin
            ke=100.
            color=blue
         end
         3: begin
            ke=1000.
            color=violet
         end
      endcase
      

      m0c2=938.                 ; rest mass in MeV

      beta=sqrt(1.0-((m0c2/(m0c2+ke))^2.))
      gamma=1./sqrt(1-beta^2.)  ;

      
      m=1.67e-24
      e=4.8e-10                 ; statcoulomb
      c=3.0e10

      b_nt=10.^(dindgen(40000)/10000)
      b_t=b_nt*1e-9             ; tesra
      b_g=b_t*1e4               ; gauss

      larmor3=gamma*m*c*beta*c/(e*b_g)/6.371e8; Earth Radius

      ss=min(where(b_nt ge 5))
      if k eq 0 then oplot, b_nt(ss)*[1,1],yr, lines=2
      if k eq 0 then xyouts, b_nt(ss),1100, 'L1', al=0.8
      if k eq 0 then xyouts, b_nt(ss), larmor3(ss)*1, '!C'+string((larmor3(ss)*1.), format='(f4.1)'), al=1.20, color=red
      if k eq 2 then xyouts, b_nt(ss), larmor3(ss)*1, string(round(larmor3(ss)*1.), format='(i2)'), al=-0.05, color=blue
      if k eq 3 then xyouts, b_nt(ss), larmor3(ss)*1, string(round(larmor3(ss)*1.), format='(i3)'), al=-0.00, color=violet

      ss=min(where(b_nt ge 10))
      if k eq 0 then oplot, b_nt(ss)*[1,1],yr, lines=2
      if k eq 0 then xyouts, b_nt(ss),1100, '~100 Rs', al=0.2
      if k eq 0 then xyouts, b_nt(ss), larmor3(ss)*1, '!C'+string(larmor3(ss)*1, format='(f3.1)'), al=1.05, color=red
      if k eq 2 then xyouts, b_nt(ss), larmor3(ss)*1, string(round(larmor3(ss)*1), format='(i2)'), al=-0.20, color=blue
      if k eq 3 then xyouts, b_nt(ss), larmor3(ss)*1, string(round(larmor3(ss)*1), format='(i2)'), al=-0.20, color=violet


      ss=min(where(b_nt ge 100))
      if k eq 0 then oplot, b_nt(ss)*[1,1],yr, lines=2
      if k eq 0 then xyouts, b_nt(ss),1100, 'GSO', al=0.5
      if k eq 0 then xyouts, b_nt(ss), larmor3(ss)*1, '!C'+string(larmor3(ss)*1, format='(f4.2)'), al=1.05, color=red
      if k eq 2 then xyouts, b_nt(ss), larmor3(ss)*1, string(larmor3(ss)*1, format='(f3.1)'), al=-0.10, color=blue
      if k eq 3 then xyouts, b_nt(ss), larmor3(ss)*1, string(larmor3(ss)*1, format='(f3.1)'), al=-0.10, color=violet

      ss=min(where(b_nt ge 1000))
      if k eq 0 then oplot, b_nt(ss)*[1,1],yr, lines=2
      if k eq 0 then xyouts, b_nt(ss),1100, '~10 Rs', al=0.5
      if k eq 0 then xyouts, b_nt(ss), larmor3(ss)*1, '!C'+string(larmor3(ss)*1, format='(f5.3)'), al=1.05, color=red
      if k eq 2 then xyouts, b_nt(ss), larmor3(ss)*1, string(larmor3(ss)*1, format='(f4.2)'), al=-0.05, color=blue
      if k eq 3 then xyouts, b_nt(ss), larmor3(ss)*1, string(larmor3(ss)*1, format='(f4.2)'), al=-0.05, color=violet



      oplot, b_nt, larmor3*1, color=color

   endfor;k

   pend,/display,bbox=0

endfor;eps

end