function rayleigh_backscatter_cross_section_spheroid, r, am, bm, D, freq_ghz, temp

common epsi, eps_spheroid, d_equi

; this function computes the rayleigh backscatter cross section (cm^2) of a prolate spheroid ice crystal with aspect ratio r and mass-dimensional realtion am*(D^bm) at a frequency of freg_ghz
; For prolate spheroid, r>1
; ; Based on paper by Rhyzkov et al., 1998 (Journal of applied meteorology, Vol 37, pg 125)
; ; note that equation 24 of Rhyzkov et al., 2011 is incorrect.  
; dielectric constant by Matrosov et al. (1996): R. F. Renking, R. A. Kropfli, and B. W. Bartram, 1996: Estimation of ice hydrometeor types and shapes from radar polarization measurements. J. Atmos. Oceanic Technol., 13, 85–96.
; Uses the Van De Hulst (1981: Light Scattering by Small Particles. Dover Publications, 470 pp.) Rayleigh approximation
; Refractive index from Matzler and Wugmuller 1987; Journal of Physics D Applied Physics · November 2000 DOI: 10.1088/0022-3727/20/12/013
; 
; Inputs:
; r=a/b - unitless ratio of b to a dimension of the prolate spheroid and should be greater than 1
; am - prefactor of the mass dimensional relatonship (units of g/(cm^bm)
; bm - exponent of the mass dimensional relationship (unitless)
; D - maximum dimension of the oblate spheroid (cm)
; freq_ghz - radar frequency in Ghz
; temperature - Kelvins
; spheroid_type - 'P' prolate.  'O' oblate
; output
; 4 element vector sigma_h, sigma_v, fh, fv
; rayleigh_backscatter_cross_section_oblate in cm^2 in a 2-element vector with [horizontal,vertical]

;r=0.6
;am=0.001 
;bm=2.0
;D=0.1   ; cm
;freq_ghz=6. 
;temp=263.
;spheroid_type='O'

;m_liq=liquid_refractive_index_mw_turner( freq_ghz, temp)
;k_sqrd_liq=abs((((m_liq^2)-1.)/((m_liq^2)+2.))^2)
;k_sqrd_ice=abs((((epsilon_solid^2)-1.)/((epsilon_solid^2)+2.))^2)
;print, k_sqrd_liq, k_sqrd_ice

cc=299792458.d*100.d ; cm/s
lambda_radar=cc/(freq_ghz*1.e9)  ; cm
lambda_radar_mm=lambda_radar*10.

; calcualte the real part of the refractive index for pure ice.   Eqn 10 of Matzler and Wugmuller
epsilon_prime=(3.1884d)+(0.00091d)*(temp-273.)

; imaginary part of the refractive index figures 5, 6, and equation 13 of Matzler and Wugmuller
if abs((temp-273.)-(-15.)) lt abs((temp-273.)-(-5.)) then begin
A=3.5d-4 & B=3.6d-5 & C=1.2d
endif else begin
  A=6.0d-4 & B=6.5d-5 & C=1.07d
endelse

epsilon_dprime=(A/freq_ghz)+B*(freq_ghz^C)

;epsilon_prime=3.18 
;epsilon_dprime=1.e-3

epsilon_solid=complex(epsilon_prime, epsilon_dprime)


D_mm=double(D)*10.d

;if spheroid_type eq 'P' then begin
;a=D_mm/2.  ; axis of rotation and maximum dimension
;b=a*r
;endif
;if spheroid_type eq 'O' then begin
b=D_mm/2.d  ; axis of rotation and maximum dimension
a=b*double(r)
;endif
;if spheroid_type ne 'O' and spheroid_type ne 'P' then stop


volume_spheroid=((4.d)/(3.d))*!dpi*(a)*(b^2) ; cubic mm
D_equi=2.d*((volume_spheroid*3.d/(4.d*!pi))^(1.d/3.d))  ; equivlent volume diameter of a sphere, mm

;mass_spheroid=am*((D_equi)^bm) ; g using D in cm and am in g/cm^bm
mass_spheroid=am*((D)^bm) ; g using D in cm and am in g/cm^bm

density_spheroid=(mass_spheroid/volume_spheroid)*1000.d  ; g per cubic cm

if density_spheroid gt 0.92d then density_spheroid=0.92d

n=(density_spheroid/0.92d)*(epsilon_solid-1.d)/(epsilon_solid+2.d)  ; Matrosov et al. (1996) equation 7

eps_spheroid=((2.d*n)+1.d)/((1.d)-n)
;print, (eps_spheroid-1.)/(eps_spheroid+2.)



;if spheroid_type eq 'P' then begin
;  e_sqrd=(1.-(b/a)^2) & e=sqrt(e_sqrd)
;  La=(((1./(2.*e))*(alog((1.+e_sqrd)/(1.-e_sqrd))))-1.)*((1.-(e_sqrd))/e_sqrd)  ; equation 6 and van de Hulst 1981 pg 71.
;endif
;if spheroid_type eq 'O' then begin
  f_sqrd=((b/a)^2)-1.d
  f=sqrt(f_sqrd)
  La=(((1.d)+f_sqrd)/f_sqrd)*((1.d)-((atan(f))/f)) ; equation 7 of Ryz et al. 1998 and van de Hulst (1981) pg 71
;endif

Lb=((1.d)-La)/2.d

zeta_a=1.d/(La+(1.d/(eps_spheroid-1.d))) ; equation 5 of Ryz et al. 1998
zeta_b=1.d/(Lb+(1.d/(eps_spheroid-1.d)))

fa=(((!pi^2)*(D_equi^3))/(6.d*lambda_radar_mm^2))*zeta_a  ; equation 4 and lambda_radar should be in mm, D_equi in mm
fb=(((!pi^2)*(D_equi^3))/(6.d*lambda_radar_mm^2))*zeta_b

cross_section_a1=4.*!pi*(abs(fa))^2
cross_section_b1=4.*!pi*(abs(fb))^2

cross_section_a2=real_part(((!Pi^5)*(D_equi^6)/(9.d*(lambda_radar_mm^4)))*(zeta_a^2))
cross_section_b2=real_part(((!Pi^5)*(D_equi^6)/(9.d*(lambda_radar_mm^4)))*((zeta_b^2)))

cross_section_h=cross_section_b2/100.d ; square cm
cross_section_v=cross_section_a2/100.d  ; square cm

fh=fb/10.d ; amplitude in cm
fv=fa/10.d

;print, d, real_part(La), real_part(f), n

d_equi=d_equi/10.

return, [cross_section_h, cross_section_v,fh,fv]
end