;IDL code to calculate the Mie/Rayleigh backscattering cross-sections and ratio
;Carey, 10/05
pro mie_ray_plot, lam, x1, x2
; lam: wavelength (cm) - user provided
; Refractive index (x1,x2) (m=n-ik where n=x1, k=x2) - user provided 
; e.g., see Battan 1973
; open text outputfiles
  openw,96,'mie_ray_error.dat'
  openw,95,'norm_mie_ray.dat'
  openw,94,'mie_ray.dat'
; do not declare s1, s2 in the main code
  nang = 2
; declare arrays for plotting
  mie=FLTARR(1000)
  ray=FLTARR(1000)
  mor=FLTARR(1000)
  perror=FLTARR(1000)
  ind=0
; loop through x: size parameter
  for x=0.01,10.0,0.01 do begin 
; calculate diameter, D, diameter from size parameter, x = 2*pi*a/lam = pi*d/lam
; first get pi
  pi = 4.D0*atan(1.D0)
  d = (x*lam)/pi
; refractive index needs to be double complex
  refrel = dcomplex(x1,x2)
; call bohren and huffman mie program in IDL (bhmie.pro)
  bhmie, x,refrel,nang,s1,s2,qext,qsca,qback,gsca
;  print,x, refrel, qext,qsca,qback,gsca
; Multiply Qback by factor 4*pi to be consistent with 
; radar backscatter efficiency(Qback) and multiply Qback by pi*r^2 = pi*(d/2)^2
; to get radar backscatter cross-section for Mie scattering.  
; Net result is to multiply Qback by 4*pi for normalized radar backscatter
; cross-section and to multiply Qback by (pi*d)^2 to get 
; radar backscatter cross-section.
;
; get normalized Mie backscattering cross-section
  nsigmie=qback*(4*pi)
; get Mie backscatter cross-section (Eqn. 5.63 Stephens 1994)
  sigmamie = qback*(pi*d)^2
; calculate Rayleigh back-scatter cross-section (Eqn. 5.65 Stephens 1994)
  K2=abs((refrel^2-1)/(refrel^2+2))^2
; Rayleigh back-scattering cross-section
  sigmaray = pi^5/lam^4*K2*d^6 
; normalized Rayleigh back-scattering cross-section
  nsigray = sigmaray/(pi*(d/2)^2)
; print out results
  printf, 94, d, sigmamie, sigmaray
  printf, 95, x, nsigmie, nsigray, nsigmie/nsigray 
  printf, 96, x, abs(sigmamie-sigmaray)/sigmamie*100. 
; copy to arrays for plotting
  mie(ind)=nsigmie
  ray(ind)=nsigray
  mor(ind)=sigmamie/sigmaray
  perror(ind)=abs(sigmamie-sigmaray)/sigmamie*100.
  ind=ind+1
  endfor

  !P.Font=1
  Device, SET_FONT = "Helvetica Bold", /TT_FONT
  SET_PLOT, 'PS'
  Device, filename='nsigma_mie_ray.ps', /COLOR, /Helvetica

  plot, mie, /XLOG, /YLOG, Xticks=3, XtickName=[0.01,0.1,1,10], xr=[1,1000], $
      yr=[0.0000001,100000], XTitle='Size Parameter', $ 
      YTitle='Normalized BackScattering Cross-Section',color=0, background=255

  oplot, ray, color=0, linestyle=2
  Device, /CLOSE

  !P.Font=1
  Device, SET_FONT = "Helvetica Bold", /TT_FONT
  SET_PLOT, 'PS'
  Device, filename='mie_ovr_ray.ps', /COLOR, /Helvetica

  plot, mor, /XLOG, /YLOG, Xticks=3, XtickName=[0.01,0.1,1,10], xr=[1,1000], $
      yr=[0.00001,10], XTitle='Size Parameter', $ 
      YTitle='Sigma(Mie)/Sigma(Rayleigh)',color=0, background=255

  Device, /CLOSE

  !P.Font=1
  Device, SET_FONT = "Helvetica Bold", /TT_FONT
  SET_PLOT, 'PS'
  Device, filename='mie_ray_error.ps', /COLOR, /Helvetica

  plot, perror, /XLOG, /YLOG, Xticks=3, XtickName=[0.01,0.1,1,10], xr=[1,1000], $
      yr=[0.01,100], XTitle='Size Parameter', $ 
      YTitle='Absolute Error (%)',color=0, background=255

  Device, /CLOSE


stop
end