pro rayleigh_canting_power_law_relations, r,am, bm, freq_ghz, temp, phi, sigma, ah, bh, av, bv, ak, bk, arho, brho

; computes the power laws that relate radar cross sections and phase delays to 
; aspect ratio and mass of prolate spheroids under the Rayleigh approximation
; here axes a>b and a is assumed to be oriented vertically
; See Jung et al 2010 Equations 3-7 or so and Rhyzkov et al., 2011 and 1998.
; r - aspect ratio of oblate spheroid of b axis to a axis
; bm - exponent of the mass-dimensional reationship
; am - prefactor of tthe mass-dimensional relationship (g/cm^bm)
; freq_ghz - radar frequency in GHz
; temp - temperature in kelvins
; phi - mean canting angle
; sigma - standard deviation of the canting angle
; output
; 
; 

common epsi, eps_spheroid, d_equi


;;; test with these parameters
;r=0.3
;am=0.001
;bm=2.0
;temp=260.
;freq_ghz=6.
;phi=0.
;sigma=40. ; degrees

sigma=sigma*(!dpi/180.d)  ; radians

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

; compute a D vector - cm
D=0.03d & while(max(D)) lt 10. do D=[D,(max(D)*1.1d)] 

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

;;;;; following Jung et al.2010 - canting angle terms
A=(1.d/8.d)*((3.d)+((4.d)*cos((2.d)*phi)*exp((-2.d)*(sigma^2)))+(cos((4.d)*phi)*exp((-8.d)*(sigma^2))))
B=(1.d/8.d)*((3.d)-((4.d)*cos((2.d)*phi)*exp((-2.d)*(sigma^2)))+(cos((4.d)*phi)*exp((-8.d)*(sigma^2))))
C=(1.d/8.d)*((1.d)-(cos((4.d)*phi)*exp((-8.d)*(sigma^2))))
Ck=cos(2.d*phi)*exp(-2.d*(sigma^2))

; compute the backscatter cross sections at these D's

; initialize variables
sigma_h=dblarr(n_elements(D))
sigma_v=dblarr(n_elements(D))
fh=complexarr(n_elements(D))
fv=complexarr(n_elements(D))
M=dblarr(n_elements(D))
re_fhmfv=dblarr(n_elements(D))
sig_diff_ov_re_fhmfv=dblarr(n_elements(D))
d_equi_vec=dblarr(n_elements(D))

zhh_bracketed_term=dblarr(n_elements(D))
zvv_bracketed_term=dblarr(n_elements(d))
kdp_term=dblarr(n_elements(d))
zhv_bracketed_term=dblarr(n_elements(d))


for j=0,n_elements(D)-1 do begin
  if r lt 1. then begin
    x=rayleigh_backscatter_cross_section_spheroid(r, am, bm, D[j], freq_ghz, temp)
  endif else begin
    x=rayleigh_backscatter_cross_section_prolate_spheroid(r, am, bm, D[j], freq_ghz, temp)    
  endelse
sigma_h[j]=real_part(x[0])
sigma_v[j]=real_part(x[1])  ; for prolates, v>h
fh[j]=x[2]
fv[j]=x[3]
M[j]=am*(D[j]^bm)
re_fhmfv[j]=real_part(fh[j]-fv[j])  ; fv>fh so term should be negative
if re_fhmfv[j] gt 0. and r gt 1. then stop
sig_diff_ov_re_fhmfv[j]=(sigma_h[j]-sigma_v[j])/real_part(fh[j]-fv[j])
d_equi_vec[j]=d_equi

zhh_bracketed_term[j]=(A*(fh[j]^2))+(B*(fv[j]^2))+(2.*C*real_part(fh[j]*conj(fv[j])))
zvv_bracketed_term[j]=(B*(fh[j]^2))+(A*(fv[j]^2))+(2.*C*real_part(fh[j]*conj(fv[j])))
if zvv_bracketed_term[j] lt zhh_bracketed_term[j] and r gt 1. then stop
kdp_term[j]=Ck*re_fhmfv[j]
zhv_bracketed_term[j]=(C*((fh[j]^2)+(fv[j]^2)))+(A*(fh[j]*conj(fv[j])))+(B*(fv[j]*conj(fh[j])))

endfor

result=linfit(alog10(D),alog10(zhh_bracketed_term),yfit=fit)
ah=10.d^result[0]
bh=result[1]

result=linfit(alog10(D),alog10(zvv_bracketed_term),yfit=fit)
av=10.d^result[0]
bv=result[1]

result=linfit(alog10(D),alog10(abs(kdp_term)),yfit=fit)
ak=10.d^result[0]
if r gt 1. then ak=(-1.d)*ak
bk=result[1]

result=linfit(alog10(D),alog10(zhv_bracketed_term),yfit=fit)
arho=10.d^result[0]
brho=result[1]


; plot results to test validity of the power laws and equations
do_plot=-1
if do_plot eq 1 then begin
  set_plot, 'Z'
  erase
  device, set_resolution=[2600,2400]
  loadct, 5
  !p.background=255   ;
  ;!p.charsize=2.0
  !p.multi=[0,2,2]
  xl=0.08 & xr=0.90 & yb=0.08 & yt=0.90 & sx=0.12 & sy=0.08
  numplots_x=2 & numplots_y=2
  position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos

  pnum=0

  plot,D,zhh_bracketed_term,color=0, $
    psym=2,symsize=1.,position=pos[pnum,*],$
    ytitle='radar cross section cm^2',/xlog, /ylog,$
    title='Cross section h -  Power law fits', xtitle='Maximum D (cm)', $
    charthick=5, charsize=3., thick=10, xthick=5, ythick=5

  oplot, D, (ah*(D^bh)), psym=-1, color=0, thick=3 ; 3.5 instead of 9 gives the best fit.    ;/(9.*(lambda_radar^4)/(!pi^5))

oplot, D, zvv_bracketed_term, color=170, psym=2
oplot, D, (av*(D^bv)), psym=-1, color=170, thick=3
;
  pnum=1

  plot,D,zvv_bracketed_term,color=0, $
    psym=2,symsize=1.,position=pos[pnum,*],$
    ytitle='radar cross section cm^2',/xlog, /ylog,$
    title='Cross section v -  Power law fits', xtitle='Maximum D (cm)', $
    charthick=5, charsize=3., thick=10, xthick=5, ythick=5

  oplot, D, (av*(D^bv)), psym=-1, color=150, thick=3 ; 3.5 instead of 9 gives the best fit. ;/(9.*(lambda_radar^4)/(!pi^5))

  pnum=2

  plot,D,kdp_term,color=0, $
    psym=2,symsize=1.,position=pos[pnum,*],$
    ytitle='re_fhmfv',/xlog, /ylog,$
    title='re_fhmfv - Power law fits', xtitle='Maximum D (cm)', $
    charthick=5, charsize=3., thick=10, xthick=5, ythick=5

  oplot, D, (ak*(D^bk)), psym=-1, color=150, thick=3 ; 3.5 instead of 9 gives the best fit. ;*(!pi^2)/(6.*lambda_radar^2)

  pnum=3

  plot,D,zhv_bracketed_term,color=0, $
    psym=2,symsize=1.,position=pos[pnum,*],$
    ytitle='zhv_bracketed_term',/xlog, /ylog,$
    title='zhv_bracketed_term - Power law fits', xtitle='Maximum D (cm)', $
    charthick=5, charsize=3., thick=10, xthick=5, ythick=5

  oplot, D, (arho*(D^brho)), psym=-1, color=150, thick=3 ; 3.5 instead of 9 gives the best fit. ; *2.*(!pi^2)/(3.*lambda_radar^3)



  write_gif, '/Users/u0029340/Documents/idl_code/hid_gmi/power_law_test_canting.gif', tvrd() 
  
  
endif

return
end