;
pro peng_fit_error
;
defsysv, '!tkref', 273.16d
defsysv, '!tkmelt', 273.15d
;
;..temp. - aerosol (Method #1) fit coefficients from Tab. 1 in P14
;
fit_lna    = -14.89d
fit_b      = 4.79d
fit_c      = 0.0076d
fit_d      = 0.86d
;
;..statistical error (aka standard deviations) from Tab. 1 in P14
;
sigma_lna   = 3.49d
sigma_b     = 1.06d
sigma_c     = 0.0373d
sigma_d     = 1.07d
;
;..vectorized temperature and a representative value of n05
;
Tlow  = -34 + dindgen(21)
n05   =  16.d
;
ln_N       = fit_lna + fit_b * alog(!tkref-(Tlow+!tkmelt)) + (fit_c*(!tkref-(Tlow+!tkmelt))+fit_d) * alog(n05) 
;
;..squared relative sensitivities; compare with the analysis by Peng
;
srs_a      = 1d / ln_N^2
srs_b      = (alog(!tkref-(Tlow+!tkmelt)))^2 / ln_N^2
srs_c      = (alog(n05)*(!tkref-(Tlow+!tkmelt)))^2 / ln_N^2
srs_d      = alog(n05)^2 / ln_N^2
;
;..contributions to the relative variance; see Equation 13.9 in Young (1962)
;
rv_a      = srs_a * sigma_lna^2
rv_b      = srs_b * sigma_b^2
rv_c      = srs_c * sigma_c^2
rv_d      = srs_d * sigma_d^2
;
;..sum of the squares of the fractional standard deviations, we are in log space
;..this is Equation 13.9 in Young
;..I refer to this as the relative variance
;
rel_var   = rv_a + rv_b + rv_c + rv_d
;
;..square of the standard deviation in log space, see my notes
;
x         = ln_N^2 * rel_var
;
;..the relative standard deviation, see my notes, this is ad hoc
;
rel_error_N   = sqrt(x)
;
;..here are the components
;
rel_error_N_a = sqrt(ln_N^2*rv_a)
rel_error_N_b = sqrt(ln_N^2*rv_b)
rel_error_N_c = sqrt(ln_N^2*rv_c)
rel_error_N_d = sqrt(ln_N^2*rv_d)
;
for i = 0, n_elements(ln_N) - 1 do begin
 print, Tlow[i], exp(ln_N[i]), rv_a[i], rv_b[i], rv_c[i], rv_d[i], rel_error_N[i], rel_error_N_a[i], rel_error_N_b[i], rel_error_N_c[i], rel_error_N_d[i], format='(11f10.1)'
endfor
;
end