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='gyro_frequency.ps'

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


   xr=[1e4,1e0]
   yr=[1e-3,1e2]
;   xtitle='Distance [Rs]'
   xtitle='Magnetic Field Strength [nT]'
   ytitle='Gyrofrequency [s!U-1!N]'
   ytickname=['0.001','0.01','0.1','1.0','10','100']

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

   oplot,[10000,1],[100,100]

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

;   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.)  ;
      
      
      speed=sqrt(2*ke*1.6e-6/1.67e-24)/3.0e10
      
      m=1.67e-24
      e=4.8e-10                 ; statcoulomb
      
      r=((1+findgen(1200))/800)
      b_nt=2.36*r^(-2.06)
      b=b_nt*1e-9*1e5           ; gauss
      b_t=b_nt*1e-9             ; tesra
      au=1.5e13
      larmor3=3.3*(gamma*938)*(beta)/b_t/au
      
      gfreq=(speed*3.0e10)/(2*!pi*larmor3*au)


      ss=min(where(r ge 1.0))
      ss=max(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),110, 'L1', al=0.2
      if k eq 2 then xyouts, b_nt(ss), gfreq(ss), string(gfreq(ss), format='(f6.4)'), al=-0.05, color=blue
      if k eq 3 then xyouts, b_nt(ss), gfreq(ss), '!C'+string(gfreq(ss), format='(f6.4)'), al=1.05, color=violet


      ss=min(where(r ge 1.0/215*100))
      ss=max(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),110, '~100 Rs', al=0.8
      if k eq 2 then xyouts, b_nt(ss), gfreq(ss), string(gfreq(ss), format='(f5.3)'), al=-0.05, color=blue
      if k eq 3 then xyouts, b_nt(ss), gfreq(ss), '!C'+string(gfreq(ss), format='(f5.3)'), al=1.05, color=violet

      ss=max(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),110, 'GSO', al=0.5
      if k eq 2 then xyouts, b_nt(ss), gfreq(ss), string(gfreq(ss), format='(f4.2)'), al=-0.05, color=blue
      if k eq 3 then xyouts, b_nt(ss), gfreq(ss), '!C'+string(gfreq(ss), format='(f4.2)'), al=1.05, color=violet


      ss=min(where(r ge 1.0/215*10))
      ss=max(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),110, '~10 Rs', al=0.5
      if k eq 2 then xyouts, b_nt(ss), gfreq(ss), string(gfreq(ss), format='(f3.1)'), al=-0.05, color=blue
      if k eq 3 then xyouts, b_nt(ss), gfreq(ss), '!C'+string(gfreq(ss), format='(f3.1)'), al=1.05, color=violet


      case k of
         0:       oplot, b_nt, gfreq, color=color, lines=0
         1:       oplot, b_nt, gfreq, color=color, lines=1
         2:       oplot, b_nt, gfreq, color=color, lines=2
         3:       oplot, b_nt, gfreq, color=color, lines=0
      endcase

   endfor;k

   pend,/display,bbox=0

endfor;eps

end