pro ciret2_optimal_iwc_Dmass_ap_norm_loglininp, rho_air_in, dbz_in, vdq_in, N0_out, lambda_out, iwc_out, conc_out, $
    count_median_length_out,mass_mean_length_out, v_mass_mean_out, N0_out_err, lambda_out_err, $
    iwc_out_err, conc_out_err, $
    count_median_length_out_err,mass_mean_length_out_err, v_mass_mean_out_err, $
     iwc_exact_out, dmass_exact_out, $
    dge_out, dge_out_err, effective_radius_out, effective_radius_out_err, ext_coeff_out, ext_coeff_out_err,param_fname

;param_fname will hold the values of the parameters and covariance matrices. if it is named zzz then default values will be
; used.

;rho_air_in=0.35 & dbz_in=-25.0 & vdq_in=0.4

; assume rho_air_in input in kg/m^3
; assume dbz_in input in 10*alog10(Ze (mm^6/m^3))
; assume vdq_in input in m/s
; lambda_out and lambda_out_err output in 1/micron
; N0_out and N0_out_err output in 1/l/micron
; iwc and iwc_err output in g/m^3
; conc and conc_err output in 1/liter
; all sizes and their errors output in microns

;pro ciret2_expo_ret_sub_analyt, dbz, vdq2

common pwr_law_coeff, az, av, am, bz_in, bv, bm
common Sb, frac_err_az,frac_err_bz,frac_err_av,frac_err_bv,frac_err_am,frac_err_bm, $
r_am_bm,r_am_av,r_am_bv,r_am_az,r_am_bz,r_bm_av,r_bm_bv,r_bm_az,r_bm_bz,r_av_bv,$
r_bv_az,r_bv_bz,r_az_bz, r_av_az, r_av_bz
common Sa, r_iwc_dmass,frac_var_iwc,frac_var_dmass
common Se, frac_err_z,frac_err_vdq,r_Z_vdq

common micro_params_expo22, iwc_exact,dist_exact,conc_exact,mass_median_length_exact, $
    count_median_length_exact, $
    count_mean_length_exact, Dmass_exact

common micro_params_expo, lambda_exact,N0_exact,rho_air

common micro_params_expo23, v_mass_mean

;common means, lnmean_exp_iwc,lnmean_exp_mmsize,lnmean_z,lnmean_v

; need to calcualte the following arrays and parameters:
;
; 1:  y is the vector of measurements containing Ze and Vdq
; 2:  x_hat is the estimated state vector containing iwc_hat and Dmass_hat
; 3:  b is the vector of model parameters containing, az, av, bz, bv
; 4:  eps_y is the vector measurement errors containing sigma_sqd_Z and sigma_sqd_V
; 5:  S_b is the model parameter covariance matrix
; 6:  F is the vector of forward model equations relating x_hat and y_hat
; 7:  K_x is the matrix of describing the sensitivity of F to x-----------------
; 8:  K_b is the matrix describing the sensitivity of F to b-----------------
; 9:  S_eps is the measurement error covariance matrix
; 10: S_y is the forward model covariance matrix
; 11: x_a is the apriori state vector containing iwc_ap, Dmass_ap
; 12: S_a is the covariance matrix of the apriori data


bz=bz_in

if rho_air_in gt 0.0 and rho_air_in lt 0.6 then begin
    rho_air=rho_air_in
endif else begin
    rho_air=0.35
endelse

if vdq_in lt 0. then vdq_in=0.05

Ze=(10.^((dbz_in)/10.d))*((1.d/10.d)^6)     ; conversion to cm^6/m^3

LI_IWC=0.097*((10.^((dbz_in)/10.d))^0.60)   ; use this below to identify outliers.  The relationship is from Liu Illingworth Jam 2000

vdq=vdq_in*100.d    ; quiet air fall velocity in cm/s.
Ze_l=Ze & vdq_l=vdq
Ze=alog(Ze) & vdq=alog(vdq)
y=transpose([Ze,vdq])
y_l=transpose([Ze_l,vdq_l])

var_az=((az)*frac_err_az)^2 & var_bz=(bz*frac_err_bz)^2
var_av=((av)*frac_err_av)^2 & var_bv=(bv*frac_err_bv)^2
var_am=((am)*frac_err_am)^2 & var_bm=(bm*frac_err_bm)^2

;here:  Fill out S_b with r's and take variances defined above into S_e and S_a
S_b=fltarr(6,6)
S_b[0,0]=(frac_err_am*am)^2
S_b[0,1]=(r_am_az)*(sqrt(var_az))*(sqrt(var_am)) & S_b[1,0]=S_b[0,1]    ;covar_az_am
S_b[0,2]=(r_am_av)*(sqrt(var_av))*(sqrt(var_am)) & S_b[2,0]=S_b[0,2]    ;covar_am_av
S_b[0,3]=(r_am_bm)*(sqrt(var_bm))*(sqrt(var_am)) & S_b[3,0]=S_b[0,3]    ;covar_am_bm
S_b[0,4]=(r_am_bz)*(sqrt(var_bz))*(sqrt(var_am)) & S_b[4,0]=S_b[0,4]    ;covar_bz_am
S_b[0,5]=(r_am_bv)*(sqrt(var_bv))*(sqrt(var_am)) & S_b[5,0]=S_b[0,5]    ;covar_am_bv

S_b[1,1]=(frac_err_az*az)^2     ;var_az
S_b[1,2]=(r_av_az)*(sqrt(var_az))*(sqrt(var_av)) & S_b[2,1]=S_b[1,2]    ;covar_az_av
S_b[1,3]=(r_bm_az)*(sqrt(var_az))*(sqrt(var_bm)) & S_b[3,1]=S_b[1,3]    ;covar_az_bm
S_b[1,4]=(r_az_bz)*(sqrt(var_az))*(sqrt(var_bz)) & S_b[4,1]=S_b[1,4]    ;covar_az_bz
S_b[1,5]=(r_bv_az)*(sqrt(var_az))*(sqrt(var_bv)) & S_b[5,1]=S_b[1,5]    ;covar_az_bv

S_b[2,2]=(frac_err_av*av)^2 ;var_av
S_b[2,3]=(r_bm_av)*(sqrt(var_bm))*(sqrt(var_av)) & S_b[3,2]=S_b[2,3]    ;covar_bm_av
S_b[2,4]=(r_av_bz)*(sqrt(var_bz))*(sqrt(var_av)) & S_b[4,2]=S_b[2,4]    ;covar_bz_av
S_b[2,5]=(r_av_bv)*(sqrt(var_bv))*(sqrt(var_av)) & S_b[5,2]=S_b[2,5]    ;covar_av_bv

S_b[3,3]=(frac_err_bm*bm)^2 ;var_bm
S_b[3,4]=(r_bm_bz)*(sqrt(var_bm))*(sqrt(var_bz)) & S_b[4,3]=S_b[3,4]    ;covar_bz_bm
S_b[3,5]=(r_bm_bv)*(sqrt(var_bm))*(sqrt(var_bv)) & S_b[5,3]=S_b[3,5]    ;covar_bm_bv

S_b[4,4]=(frac_err_bz*bz)^2 ;var_bz
S_b[4,5]=(r_bv_bz)*(sqrt(var_bv))*(sqrt(var_bz)) & S_b[4,5]=S_b[4,5]

S_b[5,5]=(frac_err_bv*bv)^2

S_b=transpose(S_b)

;av=(550.+av)/2.

ciret2_expo_ret_sub_analyt, (dbz_in+6.0), vdq_in
bz=abs(6.+bz)
iwc_ap=iwc_exact  & Dmass_ap=Dmass_exact  ; set the apriori to the exact solution
iwc_ap_l=iwc_exact & Dmass_ap_l=Dmass_exact
iwc_exact_out=iwc_ap & dmass_exact_out=dmass_ap & iwc_exact_out_l=iwc_ap & dmass_exact_out_l=dmass_ap
iwc_hat=iwc_ap & Dmass_hat=Dmass_ap & iwc_hat_l=iwc_ap & Dmass_hat_l=Dmass_ap
iwc_hat=alog(iwc_hat) & Dmass_hat=alog(Dmass_hat)
x_hat_l=transpose([iwc_hat_l, Dmass_hat_l])
x_hat=transpose([iwc_hat, Dmass_hat])
x_a=x_hat & x_a_l=x_hat_l

;;;;;;;;;;;;;;;;;;;
var_iwc_ap=((iwc_ap*frac_var_iwc)^2)
var_Dmass_ap=((Dmass_ap*frac_var_dmass)^2)
S_a=fltarr(2,2) & S_a_l=fltarr(2,2)
S_a[0,0]=alog(var_iwc_ap) & S_a_l[0,0]=var_iwc_ap
S_a[1,1]=alog(var_Dmass_ap) & S_a_l[1,1]=(var_Dmass_ap)
S_a_l[0,1]=(r_iwc_dmass*(sqrt(var_iwc_ap))*(sqrt(var_Dmass_ap))) & S_a_l[1,0]=S_a[0,1]
S_a[0,1]=alog(r_iwc_dmass*(sqrt(var_iwc_ap))*(sqrt(var_Dmass_ap))) & S_a[1,0]=S_a[0,1]
S_a=transpose(S_a) & S_a_l=transpose(S_a_l);

var_Ze=(((10.^((dbz_in+6.0)/10.d))*((1.d/10.d)^6))*frac_err_z)^2
var_vdq=(((vdq_in*100.d)*frac_err_vdq))^2
S_eps1=fltarr(2,2) & S_eps1_l=fltarr(2,2)
S_eps1[0,0]=alog(var_Ze) & S_eps1[1,1]=alog(var_vdq)
S_eps1_l[0,0]=(var_Ze) & S_eps1_l[1,1]=(var_vdq)
S_eps1_l[0,1]=( r_Z_vdq*(sqrt(var_Ze))*(sqrt(var_vdq)) ) & S_eps1_l[1,0]=S_eps1_l[0,1]
S_eps1[0,1]=alog( r_Z_vdq*(sqrt(var_Ze))*(sqrt(var_vdq)) ) & S_eps1[1,0]=S_eps1[0,1]
S_eps1=transpose(S_eps1) & S_eps1_l=transpose(S_eps1_l)


iter_count=0
while iter_count lt 1000 do begin

iter_count=iter_count+1
; For a given retrieval, the above values should not change.  The loop iteration should
; start here.

; now calculate the K and K_b matrices
dZldiwc=((az*bz)/(am*bm))*((gamma(bz))/(gamma(bm)))*(( (bm+1.)/Dmass_hat_l )^(bm-bz))
dZldDmass=((az*bz)/(am*bm))*((gamma(bz))/(gamma(bm)))*iwc_hat_l*((bm+1.)^(bm-bz))*(bz-bm)*(Dmass_hat_l^(bz-bm-1))
dvdqldDmass=av*(1.+(bv/bz))*(gamma(bz+bv)/gamma(bz))*((bm+1.)^(-bv))*(bv)*(Dmass_hat_l^((bv-1.)))
dvdqldiwc=0.0

K_x_l=fltarr(2,2)
K_x_l[0,0]=dZldiwc & K_x_l[0,1]=dZldDmass
K_x_l[1,0]=dvdqldiwc & K_x_l[1,1]=dvdqldDmass
K_x_l=transpose(K_x_l)

dZdiwc=1.0
dZdDmass=bz-bm
dvdqdDmass=bv
dvdqdiwc=0.0

K_x=fltarr(2,2)
K_x[0,0]=dZdiwc & K_x[0,1]=dZdDmass
K_x[1,0]=dvdqdiwc & K_x[1,1]=dvdqdDmass
K_x=transpose(K_x)

; now calculate K_b

dZldam=((az*bz)/bm)*((gamma(bz))/(gamma(bm)))*iwc_hat_l*(((bm+1)/Dmass_hat_l)^(bm-bz))*(-1./(am^2))
dZldaz=(bz/(am*bm))*((gamma(bz))/(gamma(bm)))*iwc_hat_l*(((bm+1)/Dmass_hat_l)^(bm-bz))
dZldav=0.0

dZdam=-1./am
dZdaz=1./az
dZdav=0.0

;dzdbm
term1=(az*bz/am)*(gamma(bz))*iwc_hat_l
term2=(1./(gamma(bm)))
term3=((bm+1.)/Dmass_hat_l)^(bm-bz)
term4=(1./bm)*term2*term3
term5=(-1./((gamma(bm))^2))*(1./bm)*term3
term6=(bm-bz)*(((bm+1.)/Dmass_hat_l)^(bm-bz-1))*(1./Dmass_hat_l)
term7=term3*alog((bm+1/Dmass_hat_l))
term8=(1./(bm*gamma(bm)))*(term6+term7)
dZldbm=term1*(-term4+term5+term8)

term1=1./bm
term2=(1./gamma(bm))*gamma_deriv_bz(bm)
term3=alog(bm+1.)
term4=(bm-bz)/(bm+1.)
term5=Dmass_hat
dZdbm=-term1-term2+term3+term4-term5

;dzdbz
term1=(az/(am*bm))*iwc_hat_l/(gamma(bm))
term2=((bm+1.)/Dmass_hat_l)^(bm-bz)
term3=term2*(gamma_deriv_bz(bz))
term4=(gamma(bz))*term2*(alog((bm+1.)/Dmass_hat_l))*(-1.)
term5=term2*gamma(bz)
dZldbz=term1*(term3+term4+term5)

dZdbz=1./bz

dZldbv=0.0 & dvdqldam=0.0 & dvdqldaz=0.0
dZdbv=0.0 & dvdqdam=0.0 & dvdqdaz=0.0
dvdqldav=(((bm+1.)/Dmass_hat_l)^(-bv))*(1.+(bv/bz))*((gamma(bz+bv))/gamma(bv))
dvdqdav=1./av

dvdqldbm=av*(1.+(bv/bz))*((gamma(bz+bv))/gamma(bz))*(-bv)*(((bm+1.)/Dmass_hat_l)^(-(bv+1)))*((1./Dmass_hat_l)^(-bv))

dvdqdbm=bv/(bm+1.)

;dvdqdbz
term1=av*(((bm+1.)/Dmass_hat_l)^(-bv))
term2=(gamma(bz+bv))/gamma(bz)
term3=1.+(bv/bz)
term4=term2*bv*(-1./(bz^2))
term5=(term3/gamma(bz))*gamma_deriv_bvv_bzc(bv,bz)
term6=term3*(gamma(bv+bz))*(-1./((gamma(bz))^2))*gamma_deriv_bz(bz)
dvdqldbz=term1*(term4+term5+term6)

term1=bv/(bz*(bz+bv))
term2=(1./gamma(bz+bv))*gamma_deriv_bvv_bzc(bv,bz)
term3=(1./gamma(bz))*gamma_deriv_bz(bz)
dvdqdbz=-term1+term2-term3

;dvdqdbv
term1=av/(gamma(bz))
term2=1.+(bv/bz)
term3=(((bm+1.)/Dmass_hat_l)^(-bv))
term4=term2*(gamma(bz+bv))*term3*alog(((bm+1.)/Dmass_hat_l))*(-1.)
term5=term3*term2*(gamma_deriv_bvv_bzc(bz, bv))
term6=term3*gamma(bz+bv)*(1./bz)
dvdqldbv=term1*(term4+term5+term6)

term1=alog(bm+1.)
term2=Dmass_hat
term3=1./(bz+bv)
term4=(gamma(bz)/gamma(bz+bv))*gamma_deriv_bvv_bzc(bz,bv)
dvdqdbv=-term1+term2+term3+term4



K_b_l=fltarr(2,6)
K_b_l[0,0]=dZldam & K_b_l[0,1]=dZldaz & K_b_l[0,2]=dZldav & K_b_l[0,3]=dZldbm & K_b_l[0,4]=dZldbz & K_b_l[0,5]=dZldbv
K_b_l[1,0]=dvdqldam & K_b_l[1,1]=dvdqldaz & k_b_l[1,2]=dvdqldav & K_b_l[1,3]=dvdqldbm &
    K_b_l[1,4]=dvdqldbz & K_b_l[1,5]=dvdqldbv
K_b_l=transpose(K_b_l)


K_b=fltarr(2,6)
K_b[0,0]=dZdam & K_b[0,1]=dZdaz & K_b[0,2]=dZdav & K_b[0,3]=dZdbm & K_b[0,4]=dZdbz & K_b[0,5]=dZdbv
K_b[1,0]=dvdqdam & K_b[1,1]=dvdqdaz & k_b[1,2]=dvdqdav & K_b[1,3]=dvdqdbm &
    K_b[1,4]=dvdqdbz & K_b[1,5]=dvdqdbv
K_b=transpose(K_b)

; now the F vector which holds Z and vdq calculated from x_hat using the forward model equations
Z_hat_l=((az*bz)/(am*bm))*(gamma(bz)/gamma(bm))*iwc_hat_l*(((bm+1.)/Dmass_hat_l)^(bm-bz))
vdq_hat_l=av*(((bm+1)/Dmass_hat_l)^(-bv))*(1.+(bv/bz))*(gamma(bz+bv)/gamma(bz))
F_l=fltarr(2) & F_l[0]=Z_hat_l & F_l[1]=vdq_hat_l & F_l=transpose(F_l)

Z_hat=iwc_hat+alog((az*bz)/(am*bm))+alog(gamma(bz)/gamma(bm))+(bm-bz)*alog(bm+1.)-(bm-bz)*Dmass_hat
vdq_hat=alog(av)-bv*alog(bm+1.)+bv*Dmass_hat+alog(1.+(bv/bz))+alog(gamma(bz+bv)/gamma(bz))
F=fltarr(2) & F[0]=Z_hat & F[1]=vdq_hat & F=transpose(F)

;Z_test=(az*N0_exact*bz*(gamma(bz)))/(lambda_exact^(bz+(1.d)))
;vdq_test=av*(1./(lambda_exact^(bv)))*((1.d)+(bv/bz))*((gamma(bz+bv))/gamma(bz))

S_y_l=(double(S_eps1_l))+((((double(K_b_l)))##(double(S_b)))##(((transpose(double(K_b_l))))))

S_y=(double(S_eps1))+((((double(K_b)))##(double(S_b)))##(((transpose(double(K_b))))))
;+(alog(abs(K_b##S_b##(transpose(K_b)))))
;S_y=alog((exp(S_eps1))+((abs(K_b))##S_b##(abs((transpose(K_b))))))
;S_y=alog((abs(K_b))##S_b##(abs((transpose(K_b)))))
S_eps=S_y & S_eps_l=S_y_l

; expression 5.8 in Rogers


term1=(invert(double(S_a_l)))##(double(x_hat_l-x_a_l))
term2=transpose(double(K_x_l))##(invert(double(S_eps_l)))##(double(y_l-F_l))
term3=invert((invert(double(S_a_l)))+((transpose(double(K_x_l)))##(invert(double(S_eps_l)))##(double(K_x_l))))
x_hat_58_l=x_hat_l+(term3)##((term2)-(term1))

term1=(invert(double(S_a)))##(double(x_hat-x_a))
term2=transpose(double(K_x))##(invert(double(S_eps)))##(double(y-F))
term3=invert((invert(double(S_a)))+((transpose(double(K_x)))##(invert(double(S_eps)))##(double(K_x))))
x_hat_58=x_hat+(term3)##((term2)-(term1))

;x_hat_58=alog((exp((term3)##((term2)-(term1))))*(exp(x_hat)))

x_hat_l=x_hat_58_l
x_hat=x_hat_58

;x_hat=x_new
iwc_hat_l=x_hat_l[0]
Dmass_hat_l=x_hat_l[1]

iwc_hat=x_hat[0]
Dmass_hat=x_hat[1]


;no use zyy
;if Dmass_hat gt 0.0 then begin
;
;if n_elements(iwc_hat_out) eq 0 then begin
;   iwc_hat_out=iwc_hat
;   Dmass_hat_out=Dmass_hat
;   z_out=10.*alog10(F[0]*1.e6)
;   vdq_out=F[1]
;endif else begin
;   iwc_hat_out=[iwc_hat_out,iwc_hat]
;   Dmass_hat_out=[Dmass_hat_out, Dmass_hat]
;   z_out=[z_out, 10.*alog10(F[0]*1.e6)]
;   vdq_out=[vdq_out,F[1]]
;endelse
;end zyy

; calculate the covariance of the solution using equation 5.14 from Rogers
S_hat_l=invert((((transpose(double(K_x_l)))##(invert(double(S_eps_l)))##(double(K_x_l))))+ invert(double(S_a_l)))
S_hat=invert((((transpose(double(K_x)))##(invert(double(S_eps)))##(double(K_x))))+ invert(double(S_a)))
matrix1=invert(S_eps+((K_x##S_a)##(transpose(K_x))))
matrix2=S_a##(transpose(K_x))
matrix3=matrix2##matrix1
matrix4=matrix3##K_x
matrix5=matrix4##S_a
S_hat2=S_a-matrix5
;S_hat=invert((((transpose(exp(K_x)))##(invert(exp(S_eps)))##(exp(K_x))))+ invert(exp(S_a)))

; calcualte a convergence parameter using equation 5.33 from Rogers

S_delta_y=S_eps##(invert( (K_x##S_a##(transpose(K_x))+S_eps   )))##S_eps

F_prev=F

; now the F vector which holds Z and vdq calculated from x_hat using the forward model equations
Z_hat_l=((az*bz)/(am*bm))*(gamma(bz)/gamma(bm))*iwc_hat_l*(((bm+1.)/Dmass_hat_l)^(bm-bz))
vdq_hat_l=av*(((bm+1)/Dmass_hat_l)^(-bv))*(1.+(bv/bz))*(gamma(bz+bv)/gamma(bz))
F_l[0]=Z_hat_l & F_l[1]=vdq_hat_l

Z_hat=iwc_hat+alog((az*bz)/(am*bm))+alog(gamma(bz)/gamma(bm))+(bm-bz)*alog(bm+1.)-(bm-bz)*Dmass_hat
vdq_hat=alog(av)-bv*alog(bm+1.)+bv*Dmass_hat+alog(1.+(bv/bz))+alog(gamma(bz+bv)/gamma(bz))
F[0]=Z_hat & F[1]=vdq_hat
;F=fltarr(2) & & F=transpose(F)

;m_test=fltarr(1)
m_test=(transpose(F-F_prev))##(invert(S_delta_y))##(F-F_prev)

;print, m_test
;print,  sqrt(S_hat[0,0]), sqrt(S_hat[1,1])
if abs(m_test[0]) lt 1.e-02 then iter_count=1111

if max(F-F_prev) eq 0.0 then iter_count=1111
    ;endif ; check for Dmass_hat gt 0
if exp(Dmass_hat) lt 0.0 then iter_count=1112
;if iter_count ge 10 then iter_count=1112


;endfor
endwhile

;S_hat[0,0]=exp(S_hat[0,0]) & S_hat[1,1]=exp(S_hat[1,1])

;factor=float(fix(max([(abs(iwc_hat)),(abs(S_hat[0,0]))])))+1.
factor=float(fix(abs(S_hat[0,0])))+1.
;factor=10.

var_iwc=S_hat_l[0,0]
var_dmass=S_hat_l[1,1]

;var_iwc=(exp((2.*iwc_hat)+(2.*S_hat[0,0])))-(exp((2.*iwc_hat)+S_hat[0,0]));*((exp(logvariwc))-1.)
;var_iwc=(((exp((2.*(iwc_hat+(factor/2.)))+(2.*(S_hat[0,0]+(factor)))))-(exp((2.*(iwc_hat+(factor/2.)))+(S_hat[0,0]+(factor)))))*exp(-factor))
;var_iwc=(exp((2.*dmass_hat)+(S_hat[0,0]^2)))*(exp(S_hat[0,0]^2)-1.)
;var_dmass=(exp((2.*dmass_hat)+S_hat[1,1]))
factor=float(fix(abs(S_hat[1,1])))+1.
;var_dmass=(((exp((2.*(dmass_hat+(factor/2.)))+(2.*(S_hat[1,1]+(factor)))))-(exp((2.*(dmass_hat+(factor/2.)))+(S_hat[1,1]+(factor)))))*exp(-factor))
; if the solution converged, calculate the error in the derived quantities using a monte carlo approach
;print, 'vdq and dbz',dbz_in, vdq_in, s_hat[0,0],iwc_hat,s_hat[1,1],dmass_hat
;print, exp(iwc_hat), iwc_exact
;print, exp(dmass_hat), dmass_exact
;print, (sqrt(var_iwc))/exp(iwc_hat),(sqrt(var_dmass))/exp(dmass_hat)
if (sqrt(var_iwc))/exp(iwc_hat) gt 2 then begin
    print, 'shit'
endif
;print, sigma_iwc(exp(iwc_hat),S_hat[0,0])/exp(iwc_hat),((exp(S_hat[1,1]/2.)))/exp(dmass_hat)
;print, sqrt((S_hat[0,0]))/exp(iwc_hat),sqrt((S_hat[1,1]))/exp(dmass_hat)
if iter_count eq 1111 and exp(Dmass_hat) gt 0.0 and exp(iwc_hat) gt 0.0 then begin

    if var_iwc gt 0.0 then begin

    iwc_out=exp(iwc_hat)

    if (abs(iwc_out-LI_IWC)/LI_IWC) gt 1.5 then begin
       var_iwc=(10.*iwc_out)^2
       ;print, 'found outlier ',iwc_out, LI_IWC
    endif

    mass_mean_length_out=exp(Dmass_hat)
    iwc_out_err=sqrt(var_iwc);iwc_out_err=(exp(S_hat[0,0]/2.))
    mass_mean_length_out_err=sqrt(var_dmass);mass_mean_length_out_err=(exp(S_hat[1,1]/2.))
    Dmass_out=mass_mean_length_out

;if Dmass_out lt mass_mean_length_out_err then mass_mean_length_out_err=Dmass_out
;if iwc_out lt iwc_out_err then iwc_out_err=iwc_out

;the idea here is to get 50 random numbers and using the sigma from above calculate a distribution of
;retrieved quantities.  Then calculate the variance of each quantitiy using the
;distribution.

       init=randomn(some_seed,1000)
       init2=randomn(some_seed,1000)
       init3=randomn(some_seed,1000)
       effective_radius_perturb=fltarr(50)
       dge_perturb=fltarr(50)
       ext_coeff_perturb=fltarr(50)

         iwc_perturb=exp(iwc_hat)+(iwc_out_err*((init))) & iwc_perturb=iwc_perturb[where(iwc_perturb gt 0.)]
         Dmass_perturb=exp(Dmass_hat)+(mass_mean_length_out_err*((init2)))  & Dmass_perturb=Dmass_perturb[where(Dmass_perturb gt 0.)]

         lambda_perturb=(bm+1.)/Dmass_perturb
         N0_perturb=((lambda_perturb^(bm+1.))*iwc_perturb)/(am*bm*gamma(bm))


         conc_perturb=N0_perturb/lambda_perturb
         count_mean_length_perturb=1./lambda_perturb
         count_median_length_perturb=abs((alog(0.5))/(lambda_perturb))
         v_mass_mean_perturb=av*(1+(bv/bm))*(lambda_perturb^(-bv))*(gamma(bm+bv)/gamma(bm))


          for kkk=0,49 do begin
         effective_radius_perturb[kkk]=effective_radius((1.e-7)*N0_perturb[kkk],(1.e-4)*lambda_perturb[kkk])  ; input liter-microns output microns
         dge_perturb[kkk]=dge((1.e-7)*N0_perturb[kkk],(1.e-4)*lambda_perturb[kkk]) ; input liter-microns output microns output microns
         ext_coeff_perturb[kkk]=ext_coeff_expo_single_rex( (1.e-5)*N0_perturb[kkk],lambda_perturb[kkk], am, bm, av, bv, rho_air) ; input/output cgs
          endfor


    lambda_out=(bm+1.)/Dmass_out & N0_out=((lambda_out^(bm+1.))*iwc_out)/(am*bm*gamma(bm)) ; N0 has units here of #/cm/m3 with lambda in 1/cm and iwc in g/m3
    conc_out=N0_out/lambda_out ; 1/m3
    count_median_length_out=abs((alog(0.5))/(lambda_out))
    v_mass_mean_out=av*(1+(bv/bm))*(lambda_out^(-bv))*(gamma(bm+bv)/gamma(bm))
    effective_radius_out=effective_radius((1.e-7)*N0_out,(1.e-4)*lambda_out) ; input liter-microns output microns
    dge_out=dge((1.e-7)*N0_out,(1.e-4)*lambda_out) ; input liter-microns output microns
    ext_coeff_out=ext_coeff_expo_single_rex( (1.e-5)*N0_out,lambda_out, am, bm, av, bv, rho_air) ; input/output cgs

    if max(where(N0_perturb/N0_out gt 0.0 and N0_perturb/N0_out lt 5.)) ne -1 then result=moment(N0_perturb[where(N0_perturb/N0_out gt 0.0 and N0_perturb/N0_out lt 5.)], sdev=N0_out_err)
    if max(where(lambda_perturb/lambda_out gt 0.0 and lambda_perturb/lambda_out lt 5.)) ne -1 then result=$
       moment(lambda_perturb[where(lambda_perturb/lambda_out gt 0.0 and lambda_perturb/lambda_out lt 5.)], sdev=lambda_out_err)
    if max(where(conc_perturb/conc_out gt 0.0 and conc_perturb/conc_out lt 5.)) ne -1 then $
       result=moment(conc_perturb[where(conc_perturb/conc_out gt 0.0 and conc_perturb/conc_out lt 5.)], sdev=conc_out_err)
    if max(where(count_median_length_perturb/count_median_length_out gt 0.0 and count_median_length_perturb/count_median_length_out lt 22222.)) ne -1 then result=$
       moment(count_median_length_perturb[where(count_median_length_perturb/count_median_length_out gt 0.0 and count_median_length_perturb/count_median_length_out lt 22222.)], $
       sdev=count_median_length_out_err)
    if max(where(v_mass_mean_perturb/v_mass_mean_out gt 0.0 and v_mass_mean_perturb/v_mass_mean_out lt 22222.)) ne -1 then result= $
       moment(v_mass_mean_perturb[where(v_mass_mean_perturb/v_mass_mean_out gt 0.0 and v_mass_mean_perturb/v_mass_mean_out lt 22222.)], sdev=v_mass_mean_out_err)

    if max(where(dge_perturb/dge_out gt 0.0 and dge_perturb/dge_out lt 22222.)) ne -1 then result= $
       moment(dge_perturb[where(dge_perturb/dge_out gt $
       0.0 and dge_perturb/dge_out lt 22222.)], sdev=dge_out_err)

    if max(where(effective_radius_perturb/effective_radius_out gt 0.0 and effective_radius_perturb/effective_radius_out lt 22222.)) ne -1 then result= $
       moment(effective_radius_perturb[where(effective_radius_perturb/effective_radius_out gt $
       0.0 and effective_radius_perturb/effective_radius_out lt 22222.)], sdev=effective_radius_out_err)

       if n_elements(where(ext_coeff_perturb/ext_coeff_out gt $
       0.0 and ext_coeff_perturb/ext_coeff_out lt 2.)) gt 2 then begin
    if max(where(ext_coeff_perturb/ext_coeff_out gt 0.0 and ext_coeff_perturb/ext_coeff_out lt 2.)) ne -1 then result= $
       moment(ext_coeff_perturb[where(ext_coeff_perturb/ext_coeff_out gt $
       0.0 and ext_coeff_perturb/ext_coeff_out lt 2.)], sdev=ext_coeff_out_err)
       endif


endif else begin ;  if S_hat[0,0] gt 0.0 then begin
    print, 'setting N0_out to iwc_exact'
    N0_out=    iwc_exact & lambda_out=Dmass_exact
    iwc_out=iwc_exact
    conc_out=conc_exact
    mass_mean_length_out=Dmass_exact
    count_median_length_out=count_median_length_exact
    v_mass_mean_out=-9999.
    dge_out=-9999.
    effective_radius_out=-9999.
    ext_coeff_out=-9999.

    N0_out_err=    -9999. & lambda_out_err=-9999.
    iwc_out_err=-9999.
    conc_out_err=-9999.
    mass_mean_length_out_err=-9999.
    count_median_length_out_err=-9999.
    v_mass_mean_out_err=-9999.
    dge_out_err=-9999.
    effective_radius_out_err=-9999.
    ext_coeff_out_err=-9999.


endelse

endif else begin

;iwc,dist,conc,mass_median_length, count_median_length, $
;   count_mean_length, mass_mean_length
    print, 'setting N0_out to iwc_exact2'

    N0_out=    iwc_exact & lambda_out=Dmass_exact
    iwc_out=iwc_exact
    conc_out=conc_exact
    mass_mean_length_out=Dmass_exact
    count_median_length_out=count_median_length_exact
    v_mass_mean_out=v_mass_mean
    dge_out=-9999.
    effective_radius_out=-9999.
    ext_coeff_out=-9999.


    N0_out_err=    -9999. & lambda_out_err=-9999.
    iwc_out_err=-9999.
    conc_out_err=-9999.
    mass_mean_length_out_err=-9999.
    count_median_length_out_err=-9999.
    v_mass_mean_out_err=-9999.
    dge_out_err=-9999.
    effective_radius_out_err=-9999.
    ext_coeff_out_err=-9999.

endelse

if iwc_out_err eq -9999. or mass_mean_length_out_err eq -9999. or N0_out_err eq -9999. $
or conc_out_err eq -9999. or count_median_length_out_err eq -9999. $
or v_mass_mean_out_err eq -9999. then begin

    N0_out_err=    -9999. & lambda_out_err=-9999.
    iwc_out_err=-9999.
    conc_out_err=-9999.
    mass_mean_length_out_err=-9999.
    count_median_length_out_err=-9999.
    v_mass_mean_out_err=-9999.
    dge_out_err=-9999.
    effective_radius_out_err=-9999.
    ext_coeff_out_err=-9999.

endif

;print, 'N0_out ', N0_out, N0_out_err
;print, 'lambda_out', lambda_out, lambda_out_err
;print, 'iwc_out ', iwc_out, iwc_out_err
;print, 'conc_out ', conc_out, conc_out_err
;print, 'mass_mean_length_out ', mass_mean_length_out, mass_mean_length_out_err
;print,  'count_median_length_out ', count_median_length_out, count_median_length_out_err
;print, 'v_mass_mean_out ', v_mass_mean_out, v_mass_mean_out_err

if N0_out gt 0. then N0_out=N0_out*1.e-7 & N0_out_err=N0_out_err*1.e-7 ; converts to 1/liter/micron from 1/cm/m3
if lambda_out gt 0. then lambda_out=lambda_out*1.e-4 & lambda_out_err=lambda_out_err*1.e-4  ; converts from 1/cm to 1/micron
; iwc_out is in g/m^3
if conc_out gt 0. then conc_out=conc_out/1000. & conc_out_err=conc_out_err/1000.    ;1/liter
if mass_mean_length_out gt 0. then mass_mean_length_out=mass_mean_length_out*1.e4 & mass_mean_length_out_err=mass_mean_length_out_err*1.e4  ; converts to microns
if count_median_length_out gt 0. then count_median_length_out=count_median_length_out*1.e4 & count_median_length_out_err=count_median_length_out_err*1.e4   ; converts to microns
if v_mass_mean_out gt 0. then v_mass_mean_out=v_mass_mean_out*1.e-2 & v_mass_mean_out_err=v_mass_mean_out_err*1.e-2 ; converts to m/s
if ext_coeff_out gt 0. then ext_coeff_out=ext_coeff_out*1.e2 ;& ext_coeff_out_err=ext_coeff_out_err*1.e2 ; converts from 1/cm to 1/m
; dge and effective_radius are calculated and output in microns



end