pro invert_zh_zdr_kdp_for_iwc_dm_r_bm_canting_test, zh, zdr, kdp, ret_vector, unc_vector, first_guess_vector, $
zh_unc, zdr_unc, kdp_unc, alpha_guess, temperature,fx_vector,fx_vector_uncertainty,am_uncert,alpha_uncert,canting_sdv,freq_ghz,$
info_content,info_vector,num_ret_param,num_ind_obs


common observables_snow_response, dBZe_h_resp,zdr_resp,kdp_resp,rhohv_resp,precip_rate_resp,iwc_resp,alpha_resp,bm_resp,aspect_ratio_resp,Dm_resp,temperature_resp,sigma_resp

common in_situ_snow_properties, Lx, Nx, alpha, mass_agg, dm_agg, dbz_cband_agg,mass_rimed, dm_rimed, dbz_cband_rimed,mass_graup, dm_graup, dbz_cband_graup

lookup_flag=0 ; 1 means use the lookup tables in obsrevables_snow_response to calucalte Kx and Kb.  0 means use finite diffferences

alpha_assumed=alpha_guess
canting_angle_sdv_assumed=canting_sdv
canting_angle_assumed=0.

; code implements a modifed verion Ryhzkov et al 2010 or Jung et al. 2011.  See notes 3/4/2021- pages 38-47.

;first guess vector and ret_vector structure and units
;iwc (g/cm3)
;Dm - mass mean size (cm)
;r - aspect ratio of hydrometeors
;bm - M-D exponent

iwc_guess=(first_guess_vector[0])  ; g/cm3
dm_guess=(first_guess_vector[1])
r_guess=(first_guess_vector[2])
bm_guess=first_guess_vector[3]

;lambda=55.  ; mm C Band
;freq_ghz=6.

; assuming that am in the m=am*D^bm relationship is correlated to am. This is an assumption that should be tested and quantified for uncertainty.
;bm=3.7038+0.75*alog10(am) ; fit for convective clouds from heysmfield et al. (2010) from Xu and Mace (2017)
am_regress=10.^((bm_guess-4.0)/0.75)

zhzv=zhv_forward_model_rayleigh_oblate_canting(r_guess, am_regress, bm_guess, iwc_guess, Dm_guess, alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)
fx_vector=[10.*alog10(zhzv[0]), (10.*alog10(zhzv[0]))-(10.*alog10(zhzv[1])) , zhzv[2]] ; simulated measurement - what we are calculating from x.


;kdp_first=kdp_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess, iwc_guess, Dm_guess, alpha_assumed, freq_ghz, temperature)
;zhzv=zhv_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess, iwc_guess, Dm_guess, alpha_assumed, freq_ghz, temperature)
;z_iwc_forward_model_vivek, N0_guess, alpha_guess,zh_guess, iwc_guess, rho_guess ; dbz
;zdr_forward_model, zdr_guess, rho_guess, r_guess, zdr ; db
;kdp_forward_model, kdp_guess, lambda, r_guess, rho_guess, iwc_guess ; degrees per km

  y_vector=[zh, zdr, kdp] ; units of degress per km for kdp.  zh and zdr input in db
  x_vector=first_guess_vector  ; x is what we are trying to retrieve 
  x_vector[0]=alog(x_vector[0])
  x_vector[1]=alog(x_vector[1])
  
  xa_vector=x_vector  ; this needs to be replaced with in situ data statistics.  see below where sa is set


;;;;;; use the uncertainties provided.  Need to know covariances.  
; Need to add in forward model uncertainty due to assumptions.  alpha, bm-am correlation, etc.
;kb=fltarr(ny,nass)  ; sensitivity of measurements to assumptions
;sb=fltarr(nass,nass)  ; covariance of assumptions - here am and bm
;am_scaled_sdv=0.75
;bm_sdv=0.2
;sb[0,0]=(am_scaled_sdv)^2 ; assuming a fractional variance of a factor of 2.
;sb[1,1]=(bm_sdv)^2  ; estimated
;sb[0,1]=0.8*(sqrt(sb[0,0]))*(sqrt(sb[1,1]))
;s_eps_forward_model=((transpose(kb)##transpose(sb))##(kb))

; calculate  ddbz_dam, dkdp_dam, dzdr_dam because this is not in the lookup table
am_plus=am_regress+(0.25*am_regress)
zhzv_plus=zhv_forward_model_rayleigh_oblate_canting(r_guess, am_plus, bm_guess, iwc_guess, Dm_guess, alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)
am_minus=am_regress-(0.25*am_regress)
zhzv_minus=zhv_forward_model_rayleigh_oblate_canting(r_guess, am_minus, bm_guess, iwc_guess, Dm_guess, alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)
ddbz_dam=(((10.*alog10(zhzv_plus[0]))-(10.*alog10(zhzv_minus[0])))/(((am_plus))-(am_minus)))*am_regress
zdr_plus=10.*alog10(zhzv_plus[0])-10.*alog10(zhzv_plus[1])
zdr_minus=10.*alog10(zhzv_minus[0])-10.*alog10(zhzv_minus[1])
dzdr_dam=((zdr_plus-zdr_minus)/((am_plus)-(am_minus)))*am_regress
dkdp_dam=((zhzv_plus[2]-zhzv_minus[2])/((am_plus)-(am_minus)))*am_regress

if lookup_flag eq 0 then begin
  zhzv_plus=zhv_forward_model_rayleigh_oblate_canting(r_guess, am_regress, bm_guess, iwc_guess, Dm_guess, alpha_assumed+(alpha_assumed/4.), canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)
  zhzc_minus=zhv_forward_model_rayleigh_oblate_canting(r_guess, am_regress, bm_guess, iwc_guess, Dm_guess, alpha_assumed-(alpha_assumed/4.), canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)
  ddbz_dalpha=(((10.*alog10(zhzv_plus[0]))-(10.*alog10(zhzv_minus[0])))/(alpha_assumed/4.))
  zdr_plus=10.*alog10(zhzv_plus[0])-10.*alog10(zhzv_plus[1])
  zdr_minus=10.*alog10(zhzv_minus[0])-10.*alog10(zhzv_minus[1])
  dzdr_dalpha=((zdr_plus-zdr_minus)/(alpha_assumed/4.))
  dkdp_dalpha=((zhzv_plus[2]-zhzv_minus[2])/(alpha_assumed/4.))
  
  zhzv_plus=zhv_forward_model_rayleigh_oblate_canting(r_guess, am_regress, bm_guess, iwc_guess, Dm_guess, alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed+(canting_angle_sdv_assumed/4.), freq_ghz, temperature)
  zhzc_minus=zhv_forward_model_rayleigh_oblate_canting(r_guess, am_regress, bm_guess, iwc_guess, Dm_guess, alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed-(canting_angle_sdv_assumed/4.), freq_ghz, temperature)
  ddbz_dsigma=(((10.*alog10(zhzv_plus[0]))-(10.*alog10(zhzv_minus[0])))/(canting_angle_sdv_assumed/4.))
  zdr_plus=10.*alog10(zhzv_plus[0])-10.*alog10(zhzv_plus[1])
  zdr_minus=10.*alog10(zhzv_minus[0])-10.*alog10(zhzv_minus[1])
  dzdr_dsigma=((zdr_plus-zdr_minus)/(canting_angle_sdv_assumed/4.))
  dkdp_dsigma=((zhzv_plus[2]-zhzv_minus[2])/(canting_angle_sdv_assumed/4.))
  
endif else begin

; calcualte d()_dalpha.  get from lookup table
;common observables_snow_response, dBZe_h_resp,zdr_resp,kdp_resp,rhohv_resp,precip_rate_resp,iwc_resp,alpha_resp,bm_resp,aspect_ratio_resp,Dm_resp,temperature_resp,sigma_resp
iwc_i=where(abs(iwc_resp-iwc_guess*1.e6) eq min(abs(iwc_resp-iwc_guess*1.e6))) & iwc_i=iwc_i[0]
alpha_i=where(abs(alpha_resp-alpha_assumed) eq min(abs(alpha_resp-alpha_assumed))) & alpha_i=alpha_i[0]
bm_i=where(abs(bm_resp-bm_guess) eq min(abs(bm_resp-bm_guess))) & bm_i=bm_i[0]
r_i=where(abs(aspect_ratio_resp-r_guess) eq min(abs(aspect_ratio_resp-r_guess))) & r_i=r_i[0]
dm_i=where(abs(dm_resp-dm_guess) eq min(abs(dm_resp-dm_guess))) & dm_i=dm_i[0]
sigma_i=where(abs(sigma_resp-canting_angle_sdv_assumed) eq min(abs(sigma_resp-canting_angle_sdv_assumed))) & sigma_i=sigma_i[0]
t_i=where(abs(temperature_resp-temperature) eq min(abs(temperature_resp-temperature))) & t_i=t_i[0]

if alpha_i gt 0 and alpha_i lt n_elements(alpha_resp)-1 then begin
  alpha_i_plus=alpha_i+1 & alpha_i_minus=alpha_i-1
endif else begin
  if alpha_i eq 0 then begin
    alpha_i_plus=alpha_i+1 & alpha_i_minus=alpha_i
  endif else begin
    alpha_i_plus=alpha_i & alpha_i_minus=alpha_i-1
  endelse
endelse

ddbz_dalpha=(dbze_h_resp[iwc_i,alpha_i_plus,bm_i,r_i, dm_i,t_i, sigma_i]-dbze_h_resp[iwc_i,alpha_i_minus,bm_i,r_i, dm_i,t_i, sigma_i])/(alpha_resp[alpha_i_plus]-alpha_resp[alpha_i_minus])
dkdp_dalpha=(kdp_resp[iwc_i,alpha_i_plus,bm_i,r_i, dm_i,t_i, sigma_i]-kdp_resp[iwc_i,alpha_i_minus,bm_i,r_i, dm_i,t_i, sigma_i])/(alpha_resp[alpha_i_plus]-alpha_resp[alpha_i_minus])
dzdr_dalpha=(zdr_resp[iwc_i,alpha_i_plus,bm_i,r_i, dm_i,t_i, sigma_i]-zdr_resp[iwc_i,alpha_i_minus,bm_i,r_i, dm_i,t_i, sigma_i])/(alpha_resp[alpha_i_plus]-alpha_resp[alpha_i_minus])

if sigma_i gt 0 and sigma_i lt n_elements(sigma_resp)-1 then begin
  sigma_i_plus=alpha_i+1 & sigma_i_minus=sigma_i-1
endif else begin
  if sigma_i eq 0 then begin
    sigma_i_plus=sigma_i+1 & sigma_i_minus=sigma_i
  endif else begin
    sigma_i_plus=sigma_i & sigma_i_minus=sigma_i-1
  endelse
endelse

ddbz_dsigma=(dbze_h_resp[iwc_i,alpha_i,bm_i,r_i, dm_i,t_i, sigma_i_plus]-dbze_h_resp[iwc_i,alpha_i,bm_i,r_i, dm_i,t_i, sigma_i_minus])/(sigma_resp[sigma_i_plus]-sigma_resp[sigma_i_minus])
dkdp_dsigma=(kdp_resp[iwc_i,alpha_i,bm_i,r_i, dm_i,t_i, sigma_i_plus]-kdp_resp[iwc_i,alpha_i,bm_i,r_i, dm_i,t_i, sigma_i_minus])/(sigma_resp[sigma_i_plus]-sigma_resp[sigma_i_minus])
dzdr_dsigma=(zdr_resp[iwc_i,alpha_i,bm_i,r_i, dm_i,t_i, sigma_i_plus]-zdr_resp[iwc_i,alpha_i,bm_i,r_i, dm_i,t_i, sigma_i_minus])/(sigma_resp[sigma_i_plus]-sigma_resp[sigma_i_minus])

endelse   ; lookup flag


kb=fltarr(3,3)  ;(ny,nass)
;kb[0,0]=ddbz_dam & kb[1,0]=dzdr_dam & kb[2,0]=dkdp_dam
;kb[0,1]=ddbz_dalpha & kb[1,1]=dzdr_dalpha & kb[2,1]=dkdp_dalpha
;kb[0,2]=ddbz_dsigma & kb[1,2]=dzdr_dsigma & kb[2,2]=dkdp_dsigma

kb[0,0]=ddbz_dam & kb[0,1]=dzdr_dam & kb[0,2]=dkdp_dam
kb[1,0]=ddbz_dalpha & kb[1,1]=dzdr_dalpha & kb[1,2]=dkdp_dalpha
kb[2,0]=ddbz_dsigma & kb[2,1]=dzdr_dsigma & kb[2,2]=dkdp_dsigma

sb=fltarr(3,3)
sb[0,0]=(am_uncert)^2  ; 20% uncertainty in am from correlation with bm
sb[1,1]=(alpha_uncert)^2  ; uncertainty in alpha
sb[2,2]=(canting_sdv)^2 ; degress uncertainty in standard deviation of canting angle.  I have no clue wtf this should be.

S_forward_model=((transpose(kb)##transpose(sb))##(kb))

  s_eps=fltarr(3,3) & s_eps[0,0]=zh_unc^2 & s_eps[1,1]=zdr_unc^2 & s_eps[2,2]=kdp_unc^2
  S_y=s_eps;+s_forward_model
  
  s_y[0,1]=0.
  s_y[0,2]=0.
  s_y[1,2]=0.
  if min(finite(s_y)) eq 0 then return

;;;;;; these are guesses of the covariances of the retrieved quantities.  Need to come up with these from in situ data and determine covariances
;common in_situ_snow_properties, Lx, Nx, alpha, mass_agg, dm_agg, dbz_cband_agg,mass_rimed, dm_rimed, dbz_cband_rimed,mass_graup, dm_graup, dbz_cband_graup
i=where(abs(dbz_cband_rimed-zh) lt 5.)
if n_elements(i) lt 10 then i=where(dbz_cband_rimed gt max(dbz_cband_rimed)/10.)
mass_var=(stddev(mass_rimed[i])/mean(mass_rimed[i]))^2
dm_var=((stddev(dm_rimed[i])/mean(dm_rimed[i])))^2

    Sa=fltarr(4,4)
    Sa[0,0]=mass_var/4.  ;1.^2  ; frational variance of iwc - factor of 2.  
    Sa[1,1]=dm_var/4.  ;1.^2
    Sa[2,2]=0.15^2
    Sa[3,3]=0.15^2
    if min(finite(sa)) eq 0 then stop

;    xa_vector[0]=alog(mean(mass_rimed[i]))      ; =x_vector  ; this needs to be replaced with in situ data statistics.  see below where sa is set
;    xa_vector[1]=alog(mean(dm_rimed[i]))

;;;;;; use finite diferences to get the first derivatives of the observables in terms of the retrieved quantities.  The Jakobian - Kx matrix
      
; dzhdiwc
; dzdrdiwc
; dkdpdiwc

;zhzv=zhv_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess, iwc_guess+(iwc_guess/4.), Dm_guess, alpha_guess, freq_ghz, temperature)
;kdp_plus=kdp_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess, iwc_guess+(iwc_guess/4.), Dm_guess, alpha_guess, freq_ghz, temperature)
;zh_plus=zhzv[0]
;zv_plus=zhzv[1]
;zhzv=zhv_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess, iwc_guess-(iwc_guess/4.), Dm_guess, alpha_guess, freq_ghz, temperature)
;kdp_minus=kdp_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess, iwc_guess-(iwc_guess/4.), Dm_guess, alpha_guess, freq_ghz, temperature)
;zh_minus=zhzv[0]
;zv_minus=zhzv[1]
;
;iwc_diff=((iwc_guess+(iwc_guess/4.))-(iwc_guess-(iwc_guess/4.)))/iwc_guess ; fractional change
;
;dzhdiwc=(10.*alog10(zh_plus)-10.*alog10(zh_minus))/iwc_diff  ; db per frational change in iwc
;dzdrdiwc=(10.*alog10(zh_plus/zv_plus)-10.*alog10(zh_minus/zv_minus))/iwc_diff
;dkdpdiwc=(kdp_plus-kdp_minus)/iwc_diff
if lookup_flag eq 0 then begin
  
  yy=zhv_forward_model_rayleigh_oblate_canting(r_guess, am_regress, bm_guess, iwc_guess+(iwc_guess/4.), Dm_guess+(dm_guess/4.), alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)
  zh_plus=yy[0]
  zv_plus=yy[1]
  kdp_plus=yy[2]
  
  yy=zhv_forward_model_rayleigh_oblate_canting(r_guess, am_regress, bm_guess, iwc_guess-(iwc_guess/4.), Dm_guess+(dm_guess/4.), alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)
  zh_minus=yy[0]
  zv_minus=yy[1]
  kdp_minus=yy[2]
  
  iwc_diff=((iwc_guess+(iwc_guess/4.))-(iwc_guess-(iwc_guess/4.)))/iwc_guess ; fractional change
  
  dzhdiwc=(10.*alog10(zh_plus)-10.*alog10(zh_minus))/iwc_diff  ; db per frational change in iwc
  dzdrdiwc=(10.*alog10(zh_plus/zv_plus)-10.*alog10(zh_minus/zv_minus))/iwc_diff
  dkdpdiwc=(kdp_plus-kdp_minus)/iwc_diff
  
endif else begin

if iwc_i gt 0 and iwc_i lt n_elements(iwc_resp)-1 then begin
  iwc_i_plus=iwc_i+1 & iwc_i_minus=iwc_i-1
endif else begin
  if iwc_i eq 0 then begin
    iwc_i_plus=iwc_i+1 & iwc_i_minus=iwc_i
  endif else begin
    iwc_i_plus=iwc_i & iwc_i_minus=iwc_i-1
  endelse
endelse
denom=(iwc_resp[iwc_i_plus]-iwc_resp[iwc_i_minus])/(iwc_guess*1.e6)

dzhdiwc=(dbze_h_resp[iwc_i_plus,alpha_i,bm_i,r_i, dm_i,t_i, sigma_i]-dbze_h_resp[iwc_i_minus,alpha_i,bm_i,r_i, dm_i,t_i, sigma_i])/denom
dkdpdiwc=(kdp_resp[iwc_i_plus,alpha_i,bm_i,r_i, dm_i,t_i, sigma_i]-kdp_resp[iwc_i_minus,alpha_i,bm_i,r_i, dm_i,t_i, sigma_i])/denom
dzdrdiwc=(zdr_resp[iwc_i_plus,alpha_i,bm_i,r_i, dm_i,t_i, sigma_i]-zdr_resp[iwc_i_minus,alpha_i,bm_i,r_i, dm_i,t_i, sigma_i])/denom
endelse
;;;;;;;;;;;;;;;;;;;;;;;;;;;; mass mean D - dm
if lookup_flag eq 0 then begin
;zhzv=zhv_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess, iwc_guess, Dm_guess+(dm_guess/4.), alpha_guess, freq_ghz, temperature)
;kdp_plus=kdp_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess, iwc_guess, Dm_guess+(dm_guess/4.), alpha_guess, freq_ghz, temperature)

yy=zhv_forward_model_rayleigh_oblate_canting(r_guess, am_regress, bm_guess, iwc_guess, Dm_guess+(dm_guess/4.), alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)
zh_plus=yy[0]
zv_plus=yy[1]
kdp_plus=yy[2]
;zhzv=zhv_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess, iwc_guess, Dm_guess-(dm_guess/4.), alpha_guess, freq_ghz, temperature)
;kdp_minus=kdp_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess, iwc_guess, Dm_guess-(dm_guess/4.), alpha_guess, freq_ghz, temperature)
yy=zhv_forward_model_rayleigh_oblate_canting(r_guess, am_regress, bm_guess, iwc_guess, Dm_guess-(dm_guess/4.), alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)

zh_minus=yy[0]
zv_minus=yy[1]
kdp_minus=yy[2]
dm_diff=((dm_guess+(dm_guess/4.))-(dm_guess-(dm_guess/4.)))/dm_guess ; fractional change

dzhddm=(10.*alog10(zh_plus)-10.*alog10(zh_minus))/dm_diff  ; db per change in dm
dzdrddm=(10.*alog10(zh_plus/zv_plus)-10.*alog10(zh_minus/zv_minus))/dm_diff
dkdpddm=(kdp_plus-kdp_minus)/dm_diff
endif else begin
if dm_i gt 0 and dm_i lt n_elements(dm_resp)-1 then begin
  dm_i_plus=dm_i+1 & dm_i_minus=dm_i-1
endif else begin
  if dm_i eq 0 then begin
    dm_i_plus=dm_i+1 & dm_i_minus=dm_i
  endif else begin
    dm_i_plus=dm_i & dm_i_minus=dm_i-1
  endelse
endelse
denom=(dm_resp[dm_i_plus]-dm_resp[dm_i_minus])/dm_guess

dzhddm=(dbze_h_resp[iwc_i,alpha_i,bm_i,r_i, dm_i_plus,t_i, sigma_i]-dbze_h_resp[iwc_i,alpha_i,bm_i,r_i, dm_i_minus,t_i, sigma_i])/denom
dkdpddm=(kdp_resp[iwc_i,alpha_i,bm_i,r_i, dm_i_plus,t_i, sigma_i]-kdp_resp[iwc_i,alpha_i,bm_i,r_i, dm_i_minus,t_i, sigma_i])/denom
dzdrddm=(zdr_resp[iwc_i,alpha_i,bm_i,r_i, dm_i_plus,t_i, sigma_i]-zdr_resp[iwc_i,alpha_i,bm_i,r_i, dm_i_minus,t_i, sigma_i])/denom
endelse
;;;;;;;;;;;;;;;;;;;;;;;;;;;; aspect ratio - r
if lookup_flag eq 0 then begin

;zhzv=zhv_forward_model_rayleigh_oblate(r_guess+(r_guess/15.), am_regress, bm_guess, iwc_guess, Dm_guess, alpha_guess, freq_ghz, temperature)
;kdp_plus=kdp_forward_model_rayleigh_oblate(r_guess+(r_guess/15.), am_regress, bm_guess, iwc_guess, Dm_guess, alpha_guess, freq_ghz, temperature)
yy=zhv_forward_model_rayleigh_oblate_canting(r_guess+(r_guess/15.), am_regress, bm_guess, iwc_guess, Dm_guess, alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)
zh_plus=yy[0]
zv_plus=yy[1]
kdp_plus=yy[2]
;zhzv=zhv_forward_model_rayleigh_oblate(r_guess-(r_guess/15.), am_regress, bm_guess, iwc_guess, Dm_guess, alpha_guess, freq_ghz, temperature)
;kdp_minus=kdp_forward_model_rayleigh_oblate(r_guess-(r_guess/15.), am_regress, bm_guess, iwc_guess, Dm_guess, alpha_guess, freq_ghz, temperature)
yy=zhv_forward_model_rayleigh_oblate_canting(r_guess-(r_guess/15.), am_regress, bm_guess, iwc_guess, Dm_guess, alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)

zh_minus=yy[0]
zv_minus=yy[1]
kdp_minus=yy[2]
r_diff=(r_guess+(r_guess/15.))-(r_guess-(r_guess/15.))

dzhdr=(10.*alog10(zh_plus)-10.*alog10(zh_minus))/r_diff  ; db per change in r
dzdrdr=(10.*alog10(zh_plus/zv_plus)-10.*alog10(zh_minus/zv_minus))/r_diff
dkdpdr=(kdp_plus-kdp_minus)/r_diff
endif else begin
if r_i gt 0 and r_i lt n_elements(aspect_ratio_resp)-2 then begin
  r_i_plus=r_i+1 & r_i_minus=r_i-1
endif else begin
  if r_i eq 0 then begin
    r_i_plus=r_i+1 & r_i_minus=r_i
  endif else begin
    r_i_plus=r_i & r_i_minus=r_i-1
  endelse
endelse
denom=(aspect_ratio_resp[r_i_plus]-aspect_ratio_resp[r_i_minus])

dzhdr=(dbze_h_resp[iwc_i,alpha_i,bm_i,r_i_plus, dm_i,t_i, sigma_i]-dbze_h_resp[iwc_i,alpha_i,bm_i,r_i_minus, dm_i,t_i, sigma_i])/denom
dkdpdr=(kdp_resp[iwc_i,alpha_i,bm_i,r_i_plus, dm_i,t_i, sigma_i]-kdp_resp[iwc_i,alpha_i,bm_i,r_i_minus, dm_i,t_i, sigma_i])/denom
dzdrdr=(zdr_resp[iwc_i,alpha_i,bm_i,r_i_plus, dm_i,t_i, sigma_i]-zdr_resp[iwc_i,alpha_i,bm_i,r_i_minus, dm_i,t_i, sigma_i])/denom

endelse

;;;;;;;;;;;;;;;;;;;;;;;;;;;; M-D exponent - bm
if lookup_flag eq 0 then begin

;am_regress=10.^(((bm_guess+(bm_guess/15.))-4.)/0.75)
am_regress=10.^(((bm_guess+(bm_guess/5.))-4.)/0.75)

;zhzv=zhv_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess+(bm_guess/5.), iwc_guess, Dm_guess, alpha_guess, freq_ghz, temperature)
;kdp_plus=kdp_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess+(bm_guess/5.), iwc_guess, Dm_guess, alpha_guess, freq_ghz, temperature)
yy=zhv_forward_model_rayleigh_oblate_canting(r_guess, am_regress, bm_guess+(bm_guess/5.), iwc_guess, Dm_guess, alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)
zh_plus=yy[0]
zv_plus=yy[1]
kdp_plus=yy[2]

;am_regress=10.^(((bm_guess-(bm_guess/15.))-3.7038)/0.75)
am_regress=10.^(((bm_guess-(bm_guess/5.))-4.)/0.75)

;zhzv=zhv_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess-(bm_guess/5.), iwc_guess, Dm_guess, alpha_guess, freq_ghz, temperature)
;kdp_minus=kdp_forward_model_rayleigh_oblate(r_guess, am_regress, bm_guess-(bm_guess/5.), iwc_guess, Dm_guess, alpha_guess, freq_ghz, temperature)

yy=zhv_forward_model_rayleigh_oblate_canting(r_guess, am_regress, bm_guess+(bm_guess/5.), iwc_guess, Dm_guess, alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)
zh_minus=yy[0]
zv_minus=yy[1]
kdp_minus=yy[2]

bm_diff=(bm_guess+(bm_guess/5.))-(bm_guess-(bm_guess/5.))

dzhdbm=(10.*alog10(zh_plus)-10.*alog10(zh_minus))/bm_diff  ; db per bm
dzdrdbm=(10.*alog10(zh_plus/zv_plus)-10.*alog10(zh_minus/zv_minus))/bm_diff
dkdpdbm=(kdp_plus-kdp_minus)/bm_diff

endif else begin
if bm_i gt 0 and bm_i lt n_elements(bm_resp)-1 then begin
  bm_i_plus=bm_i+1 & bm_i_minus=bm_i-1
endif else begin
  if bm_i eq 0 then begin
    bm_i_plus=bm_i+1 & bm_i_minus=bm_i
  endif else begin
    bm_i_plus=bm_i & bm_i_minus=bm_i-1
  endelse
endelse
denom=(bm_resp[bm_i_plus]-bm_resp[bm_i_minus])

dzhdbm=(dbze_h_resp[iwc_i,alpha_i,bm_i_plus,r_i, dm_i,t_i, sigma_i]-dbze_h_resp[iwc_i,alpha_i,bm_i_minus,r_i, dm_i,t_i, sigma_i])/denom
dkdpdbm=(kdp_resp[iwc_i,alpha_i,bm_i_plus,r_i, dm_i,t_i, sigma_i]-kdp_resp[iwc_i,alpha_i,bm_i_minus,r_i, dm_i,t_i, sigma_i])/denom
dzdrdbm=(zdr_resp[iwc_i,alpha_i,bm_i_plus,r_i, dm_i,t_i, sigma_i]-zdr_resp[iwc_i,alpha_i,bm_i_minus,r_i, dm_i,t_i, sigma_i])/denom
endelse

;iwc_guess=alog(first_guess_vector[0])  ; g/cm3
;dm_guess=alog(first_guess_vector[1])
;r_guess=(first_guess_vector[2])
;bm_guess=first_guess_vector[3]
kx=fltarr(3,4)

kx[0,0]=dzhdiwc
kx[0,1]=dzhddm
kx[0,2]=dzhdr
kx[0,3]=dzhdbm

kx[1,0]=dzdrdiwc
kx[1,1]=dzdrddm
kx[1,2]=dzdrdr
kx[1,3]=dzdrdbm

kx[2,0]=dkdpdiwc
kx[2,1]=dkdpddm
kx[2,2]=dkdpdr
kx[2,3]=dkdpdbm

kx[2,*]=kx[2,*]/10.

if min(finite(kx)) eq 0 then return

;;;;;; done filling out Jakobian


;x_hat holds the retrieval quantities as the iteration proceeds.

    x_hat=x_vector

;;;; this loop itertes the x (microphysics) to get y or F(x) (observations).  loop to 30 is arbitrary.  Need to see if this should be more or less.  Convergence test causes exit from the loop,
    for j=0,100 do begin

;;;;;; this tests for convergence of the iteration
      if j gt 0 then begin
        fu=reform(Fx_vector)-reform(Fy_prev)
        m_test=(((fu)))##(((invert(S_delta_y)))##((transpose(fu))))
        fy_prev=Fx_vector
        ;if m_test gt m_test_prev or m_test lt 0. then count_diverge=1000
      endif else begin
        m_test=1.e12
        fy_prev=Fx_vector
      endelse
      
      ;;;;;;; come up with next iteration of the solution

      term1=(invert(double(100.*Sa)))##(reform(double(x_hat-xa_vector)))
      term2=(double(Kx))##(invert(double(S_y)))##(double(y_vector-Fx_vector))
      term3=invert((invert(double(1.*Sa)))+(((double(Kx)))##(invert(double(S_y)))##(transpose(double(Kx)))))
      dxhat=(term3)##((term2)-(term1))

      if min(finite(dxhat)) le 0 then begin
        return
      endif
      
      ;dxhat is the next increment in the iteration.  Add to x_hat
      
        x_hat=x_hat+dxhat
        
        arg = transpose(Kx)#invert(S_y)
        arg2 = arg#(Kx) ;& arg2[0,1:4]=0. & arg2[1,2:4]=0. & arg2[2,3:4]=0. & arg2[3,4]=0.

;;;; compute the covariance of the solution.  The diagonal is the variance of the retrieved quantities.
        Sx = invert( (invert(Sa))+arg2 )

        ;first_guess_vector=[alog(N0_guess),alpha_guess, alog(iwc_guess), rho_guess, r_guess]
        
        ;;;;;;; convert the retrieved quantities in x_hat to the retrieved physical quantities from log space (IWC and Dm).
        
;        iwc_inc=exp(x_hat[0]+(sx[0,0]/2.))  ; this form is for lognormally distributed variables
        iwc_inc=exp(x_hat[0]) ; this returns the median
;        dm_inc=exp(x_hat[1]+(sx[1,1]/2.))
        dm_inc=exp(x_hat[1])  ; this returns the median
        r_inc=x_hat[2] & if r_inc ge 1. then r_inc=0.99
        bm_inc=x_hat[3] & if bm_inc gt 3. then bm_inc=3.
        
        ;bm=3.7038+0.75*alog10(am) ; fit for convective clouds from heysmfield et al. (2010) from Xu and Mace (2017)
        ;am_regress=10.^((bm_guess-4.0)/0.75)
        
        ;am_inc=10.^(((bm_inc-(bm_inc/5.))-3.7038)/0.75)
        am_inc=10.^((bm_inc-4.0)/0.75)
        ;;;; use the new  microphysics to calculate the observables.  The goal is to have the microphysics replicate the observed quantities.  
                
;        zhzv=zhv_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc, iwc_inc, Dm_inc, alpha_assumed, freq_ghz, temperature)
;        kdp=kdp_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc, iwc_inc, Dm_inc, alpha_guess, freq_ghz, temperature)
        
        zhzv=zhv_forward_model_rayleigh_oblate_canting(r_inc, am_inc, bm_inc, iwc_inc, Dm_inc, alpha_assumed, canting_angle_assumed, canting_angle_sdv_assumed, freq_ghz, temperature)
        fx_vector=[10.*alog10(zhzv[0]), 10.*alog10(zhzv[0])-10.*alog10(zhzv[1]) , zhzv[2]] ; simulated measurement - what we are calculating from x.

        
;        fx_vector[0]=10.*alog10(zhzv[0])
;        fx_vector[1]=10.*alog10(zhzv[0])-10.*alog10(zhzv[1])
;        fx_vector[2]=kdp
        
;;;;; computations for convergence test
        term31=invert(S_y)##transpose(Kx)
        s_delta_y=S_y##(invert(Kx##term31)+S_eps)##S_y

; if convergence is achieved, then pop out of the loop.  Otherwise iterate 30 times only.
        if abs(m_test) lt 1.e-15 then begin
          j=1000
        endif
        
      m_test_prev=m_test
      
      
;;;;;;;;;;;;;;;;; redo kx      
;;;;;; use finite diferences to get the first derivatives of the observables in terms of the retrieved quantities.  The Jakobian - Kx matrix

; dzhdiwc
; dzdrdiwc
; dkdpdiwc

;zhzv=zhv_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc, iwc_inc+(iwc_inc/4.), Dm_inc, alpha_guess, freq_ghz, temperature)
;kdp_plus=kdp_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc, iwc_inc+(iwc_inc/4.), Dm_inc, alpha_guess, freq_ghz, temperature)
;zh_plus=zhzv[0]
;zv_plus=zhzv[1]
;zhzv=zhv_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc, iwc_inc-(iwc_inc/4.), Dm_inc, alpha_guess, freq_ghz, temperature)
;kdp_minus=kdp_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc, iwc_inc-(iwc_inc/4.), Dm_inc, alpha_guess, freq_ghz, temperature)
;zh_minus=zhzv[0]
;zv_minus=zhzv[1]
;
;iwc_diff=((iwc_inc+(iwc_inc/4.))-(iwc_inc-(iwc_inc/4.)))/iwc_inc ; fractional change
;
;dzhdiwc=(10.*alog10(zh_plus)-10.*alog10(zh_minus))/iwc_diff  ; db per frational change in iwc
;dzdrdiwc=(10.*alog10(zh_plus/zv_plus)-10.*alog10(zh_minus/zv_minus))/iwc_diff
;dkdpdiwc=(kdp_plus-kdp_minus)/iwc_diff
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;; mass mean D - dm
;
;zhzv=zhv_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc, iwc_inc, Dm_inc+(Dm_inc/4.), alpha_guess, freq_ghz, temperature)
;kdp_plus=kdp_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc, iwc_inc, Dm_inc+(Dm_inc/4.), alpha_guess, freq_ghz, temperature)
;zh_plus=zhzv[0]
;zv_plus=zhzv[1]
;zhzv=zhv_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc, iwc_inc, Dm_inc-(Dm_inc/4.), alpha_guess, freq_ghz, temperature)
;kdp_minus=kdp_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc, iwc_inc, Dm_inc-(Dm_inc/4.), alpha_guess, freq_ghz, temperature)
;zh_minus=zhzv[0]
;zv_minus=zhzv[1]
;
;dm_diff=((Dm_inc+(Dm_inc/4.))-(Dm_inc-(Dm_inc/4.)))/Dm_inc ; fractional change
;
;dzhddm=(10.*alog10(zh_plus)-10.*alog10(zh_minus))/dm_diff  ; db per change in dm
;dzdrddm=(10.*alog10(zh_plus/zv_plus)-10.*alog10(zh_minus/zv_minus))/dm_diff
;dkdpddm=(kdp_plus-kdp_minus)/dm_diff
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;; aspect ratio - r
;
;zhzv=zhv_forward_model_rayleigh_oblate(r_inc+(r_inc/15.), am_inc, bm_inc, iwc_inc, Dm_inc, alpha_guess, freq_ghz, temperature)
;kdp_plus=kdp_forward_model_rayleigh_oblate(r_inc+(r_inc/15.), am_inc, bm_inc, iwc_inc, Dm_inc, alpha_guess, freq_ghz, temperature)
;zh_plus=zhzv[0]
;zv_plus=zhzv[1]
;zhzv=zhv_forward_model_rayleigh_oblate(r_inc-(r_inc/15.), am_inc, bm_inc, iwc_inc, Dm_inc, alpha_guess, freq_ghz, temperature)
;kdp_minus=kdp_forward_model_rayleigh_oblate(r_inc-(r_inc/15.), am_inc, bm_inc, iwc_inc, Dm_inc, alpha_guess, freq_ghz, temperature)
;zh_minus=zhzv[0]
;zv_minus=zhzv[1]
;
;r_diff=(r_inc+(r_inc/15.))-(r_inc-(r_inc/15.))
;
;dzhdr=(10.*alog10(zh_plus)-10.*alog10(zh_minus))/r_diff  ; db per change in r
;dzdrdr=(10.*alog10(zh_plus/zv_plus)-10.*alog10(zh_minus/zv_minus))/r_diff
;dkdpdr=(kdp_plus-kdp_minus)/r_diff
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;; M-D exponent - bm
;
;;am_inc=10.^(((bm_inc+(bm_inc/15.))-3.7038)/0.75)
;am_inc=10.^(((bm_inc+(bm_inc/5.))-3.7038)/0.75)
;
;zhzv=zhv_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc+(bm_inc/5.), iwc_inc, Dm_inc, alpha_guess, freq_ghz, temperature)
;kdp_plus=kdp_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc+(bm_inc/5.), iwc_inc, Dm_inc, alpha_guess, freq_ghz, temperature)
;zh_plus=zhzv[0]
;zv_plus=zhzv[1]
;
;;am_inc=10.^(((bm_inc-(bm_inc/15.))-3.7038)/0.75)
;am_inc=10.^(((bm_inc-(bm_inc/5.))-3.7038)/0.75)
;
;zhzv=zhv_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc-(bm_inc/5.), iwc_inc, Dm_inc, alpha_guess, freq_ghz, temperature)
;kdp_minus=kdp_forward_model_rayleigh_oblate(r_inc, am_inc, bm_inc-(bm_inc/5.), iwc_inc, Dm_inc, alpha_guess, freq_ghz, temperature)
;zh_minus=zhzv[0]
;zv_minus=zhzv[1]
;
;bm_diff=(bm_inc+(bm_inc/5.))-(bm_inc-(bm_inc/5.))
;
;dzhdbm=(10.*alog10(zh_plus)-10.*alog10(zh_minus))/bm_diff  ; db per bm
;dzdrdbm=(10.*alog10(zh_plus/zv_plus)-10.*alog10(zh_minus/zv_minus))/bm_diff
;dkdpdbm=(kdp_plus-kdp_minus)/bm_diff
;
;;iwc_guess=alog(first_guess_vector[0])  ; g/cm3
;;dm_guess=alog(first_guess_vector[1])
;;r_guess=(first_guess_vector[2])
;;bm_guess=first_guess_vector[3]
;kx=fltarr(3,4)
;
;kx[0,0]=dzhdiwc
;kx[0,1]=dzhddm
;kx[0,2]=dzhdr
;kx[0,3]=dzhdbm
;
;kx[1,0]=dzdrdiwc
;kx[1,1]=dzdrddm
;kx[1,2]=dzdrdr
;kx[1,3]=dzdrdbm
;
;kx[2,0]=dkdpdiwc
;kx[2,1]=dkdpddm
;kx[2,2]=dkdpdr
;kx[2,3]=dkdpdbm

;;;;;; done filling out Jakobian      
;;;;;;;;;;;;;;;;;;;; done witht kx redo
    endfor
    
;  redefifne Sy to include the forward model uncertainty to get an accurate error of the retrieval    
    S_y=s_eps+s_forward_model
    arg = transpose(Kx)#invert(S_y)
    arg2 = arg#(Kx) ;& arg2[0,1:4]=0. & arg2[1,2:4]=0. & arg2[2,3:4]=0. & arg2[3,4]=0.

    Sx = invert( (invert(Sa))+arg2 )

    
; pull the quantities for output.    
    
    ret_vector=x_hat
    
;    ret_vector[0]=exp(x_hat[0]+(sx[0,0]/2.))
;    ret_vector[1]=exp(x_hat[1]+(sx[1,1]/2.))
    ret_vector[0]=exp(x_hat[0]) ;+(sx[0,0]/2.))  ; return the median values
   ret_vector[1]=exp(x_hat[1])  ;+(sx[1,1]/2.))
     unc_vector=fltarr(n_elements(ret_vector))
     
     ;lwp_cloud_inc_unc=sqrt((exp((2.*x[0])+sx[0,0]))*((exp(sx[0,0]))-1.)) & fractional_unc_lwp_cloud=lwp_cloud_inc_unc/lwp_cloud_inc

;;;; need to determine correlations, informtion content, number of independent obs, etc.,  

     unc_vector[0]=(sqrt((exp((2.*x_hat[0])+sx[0,0]))*((exp(sx[0,0]))-1.)))/ret_vector[0] ; fractional uncertainties.  
     unc_vector[1]=(sqrt((exp((2.*x_hat[1])+sx[1,1]))*((exp(sx[1,1]))-1.)))/ret_vector[1]
     unc_vector[2]=sqrt(x_hat[2])/ret_vector[2]
     unc_vector[3]=sqrt(x_hat[3])/ret_vector[3]
     
     fx_vector_uncertainty=[s_y[0,0],s_y[1,1],s_y[2,2]]

;;;;;; Need to include convergence test in output. 
;
; 
info_content=(alog(determ((Sa)##(invert(sx, /double)))))/alog(2.)
S_info=((((Sa)##(invert(sx, /double)))))

iwc_info=alog(S_info[0,0])/alog(2.)
dm_info=alog(S_info[1,1])/alog(2.)
r_info=alog(S_info[2,2])/alog(2.)
bm_info=alog(S_info[3,3])/alog(2.)

info_vector=[iwc_info,dm_info,r_info,bm_info]

sqrt_sa_phys=sqrt_matrix(sa)
inverse_sqrt_sa_phys=invert(sqrt_sa_phys)
sqrt_sy=sqrt_matrix(s_y)
inverse_sqrt_sy=invert(sqrt_sy)
k_tilda=inverse_sqrt_sy##transpose(kx)##inverse_sqrt_sa_phys
lambda_matrix=k_tilda##transpose(k_tilda)
svdc, lambda_matrix, W, U, V
num_ret_param=n_elements(where(W gt 1.))
num_ind_obs=0. & for ii=0,n_elements(w)-1 do num_ind_obs=num_ind_obs+(w[ii]/(w[ii]+1.))


return
    end