function zhv_forward_model_rayleigh_oblate_canting, r, am, bm, iwc, Dm, alpha, phi, sigma_in, freq_ghz, temp

; computes horizontal and vertical polarization Z assuming Rayleigh scattering from oblates oriented in the horizontal plane (r<1) or prolates in the vertical plane (r>1)
;  and low incidence angle given...
; r - aspect ratio of oblate spheroids if r<1 and prolate spheroids r>1.  For r>1, Zv>Zh and kdp<0
; Mass-D relationship am*D^bm in cgs units
; iwc - g/cm3
; Dm - mass mean maximum dimension, cm
; alpha - shape ffactor of assumed modified gamma distribuiton off snow
; freq_ghz - radar frequency in GHz
; temp - temperature in Kelvins
; phi is the mean canting angle
; sigma_in is the standard deviatio of the canting angle in degrees
; results returned in mm^6/m^3
;
; Adpated from  Jung et al., 2010.

common pwr_law, ah, bh, av, bv

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

sigma=sigma_in
;;; calc K^2 for liquid at this freq

m_liq=liquid_refractive_index_mw_turner( freq_ghz, temp)
freq_ghz=double(freq_ghz)
k_sqrd_liq=abs((((m_liq^2)-1.)/((m_liq^2)+2.d))^2)

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

; get the power law relations - equation 16 of Ryzkhov et al. 1998 and page 39 of Notebook (3/4/2021-)
;rhyzkov_power_law_relations, r,am, bm, freq_ghz, temp, ah, bh, av, bv, a2, b2, a3, b3

; get power law relationship for backscatter cross section functions.  These are the terms that include the complex amplitudes functions in equations 3-10 of Jung et al. 2010
rayleigh_canting_power_law_relations, r,am, bm, freq_ghz, temp, phi, sigma, ah, bh, av, bv, ak, bk, arho, brho
;rayleigh_canting_power_law_relations_ryz, r,am, bm, freq_ghz, temp, phi, sigma, ah, bh, av, bv, ak, bk, arho, brho

D0=Dm/(alpha+bm+1.d)
N0=iwc/(am*(D0^(bm+1.d))*gamma(alpha+bm+1.d))

;Zh=(1.e12)*(ah/9.*k_sqrd_liq)*N0*D0^((bh)+1.)*gamma(alpha+(bh)+1.)
;Zv=(1.e12)*(av/9.*k_sqrd_liq)*N0*D0^((bv)+1.)*gamma(alpha+(bv)+1.)

Zh=(ah*4.d*(lambda_radar^4)/((!pi^4)*k_sqrd_liq))*N0*(D0^((bh)+1.))*gamma(alpha+(bh)+1.d)
Zv=(av*4.d*(lambda_radar^4)/((!pi^4)*k_sqrd_liq))*N0*(D0^((bv)+1.))*gamma(alpha+(bv)+1.d)

;kdp=((180.*!pi*ak)/(6.*lambda_radar))*N0*(D0^((bk)+1.))*gamma(alpha+(bk)+1.)

kdp=(180.d*lambda_radar*ak/!pi)*N0*(D0^((bk)+1.d))*gamma(alpha+(bk)+1.)

rho_hv=(((4.d*(arho*lambda_radar^4))/((!pi^4)*k_sqrd_liq))*N0*(D0^((brho)+1.d))*gamma(alpha+(brho)+1.))/sqrt(Zh*Zv)

Zh=Zh*(1.d12)
Zv=Zv*(1.d12)

kdp=kdp*1.e5  ; degrees per cm to degrees per km
kdp=kdp*10.   ; test for match to observations.  Unsure  of this.

return, [Zh, Zv, kdp, rho_hv, ah, bh]
end