pro lognorm, n1,r1,s1,n2,r2,s2,nb, $ ; input numbins,width,radi,dconc,vol,sfc ; output ;---------------------------------------------------------------------- ; ; purpose: ; Return the concentration vs. radius, and the total volume and surface ; area densities for a lognormal size distribution of 1 or 2 modes ; ; input: ; nb = max # bins in size distribution, 50 to 100 is good, but test it ; ; n1,r1,s1,n2,r2,s2...the lognormal size distribution parameters, set ; n2 = 0 for a unimodal distribution. ; ; n1 = mode 1 concentration, cm-3 ; r1 = mode 1 median radius, microns ; s1 = mode 1 distribution width, unitless ; n2 = mode 2 concentration, cm-3 ; r2 = mode 2 median radius, microns ; s2 = mode 2 distribution width, unitless ; ; output: ; numbins....the number of radius bins determined for the distribution ; width......array, the width of each size bin, microns ; radi.......array, the center radius of each size bin, microns ; dconc......array, the differential concentration in each bin, cm-3 um-1 ; vol........total volume density, um3/cm3 ; sfc........total surface area density, um2/cm3 ; ; source: Mark Hervig ; ;---------------------------------------------------------------------- ;---- setup some vars pi = 3.14159 ; yes, it is ; nb = 100. ; max # of bins in size distribution scon = n1 + n2 ; total concentration cmin = 1.e-9 ; concentration cutoff, skip bins with n < cmin*scon ;---- find the max and min radii using the distribution parameters esr = alog10(r1)-s1 ; exponent of smallest r (r = 10^esr) elr = alog10(r1)+s1 ; exponent of largest r if (s2 gt 0.0) then elr = alog10(r2)+s2 ;---- define the radius bins nd = elr - esr ; number of decades, i.e., log(r) ra = 10.0^( esr + nd * ( findgen(nb) +1) /nb ) ; bin upper radii, microns rb = 10.0^( esr + nd * ( findgen(nb) +2) /nb ) ; bin lower radii r = (ra + rb) / 2.0 ; bin center radii, microns wid = rb - ra ; bin widths, microns ;---- calculate differential concentrations in each bin s1l = alog(s1) dcon = n1 * exp(-1.*(alog(r/r1))^2 / (2.0*s1l^2) ) / (r *s1l *sqrt(2.0*pi)) if (n2 gt 0.0) then begin s2l = alog(s2) dcon2 = n2 * exp(-(alog(r/r2))^2 / (2.0*s2l^2)) / (r *s2l *sqrt(2.0*pi)) dcon = dcon + dcon2 endif ;---- throw out bins with conc less than total conc * cmin conc = dcon * wid ok = where(conc gt (scon*cmin),numbins) if (numbins gt 0) then begin conc = conc(ok) radi = r(ok) dconc = dcon(ok) width = wid(ok) vol = total( (1.33333 * pi * radi^3 ) * conc) sfc = total( (4.0 * pi * radi^2 ) * conc) con = total(conc) if (con lt .95*(n1+n2)) then $ print,'integral concentration ',100.*con/(s1+s2),' % of expected in lognorm.pro' endif return end