pro gaussian_sdist, n,rm,dr,nb,cutf, numbins,wid,r,dcon,vol,sfc,re ;---------------------------------------------------------------------- ; ; purpose: ; Return the concentration vs. radius, and the total volume and surface ; area densities for a gaussian particle size distribution ; ; input: ; n......total concentration, #/cm-3 ; rm.....median radius, units = L (e.g., microns) ; dr.....distribution width, units of L (same as rm) ; nb.....max # bins in the size distribution, 50 to 100 is OK (test it) ; cutf...consider only bins with > (cutf) * (the max bin concentration) ; ~1e-4 seems to be OK (test it), set cutf = zero to ignore this step ; ; output: ; numbins....the number of radius bins accepted for the distribution (depends on cutf), scalar ; wid........width of each size bin, L, fltarr(numbins) ; r..........center radius of each size bin, L, fltarr(numbins) ; dcon.......differential concentration in each bin, per cm3 per L, fltarr(numbins) ; vol........total volume density, L3/cm3, scalar ; sfc........total surface area density, L2/cm3, scalar ; re.........effective radius (units of rm), scalar ; ; source: Mark Hervig ; ; Revisions: ; 9/ 9/2009: MH revised how the max & min radii are defined ; 12/30/2010: MH added cutf as input variable, was set in here ; ;---------------------------------------------------------------------- nb = float(nb) ;- set the max and min radii to consider ; fac = 100.0 ; rmin = rm / fac ; smallest r ; rmax = rm * fac ; largest r delr = sqrt( -2. * dr^2. * alog(frac) ) ; the +/-radius step to get to dcon = dcon*frac rmax = rm + delr ; largest r rmin = rm - delr ; smallest r if (delr gt rm) then rmin = rm *0.01 esr = alog10(rmin) ; exponent of smallest r elr = alog10(rmax) ; exponent of largest r ;- define the radius bins nd = elr - esr ; number of decades, i.e., log(r) ra = 10.d0^(esr +nd*(findgen(nb)+1.)/nb) ; bin upper radii, microns rb = 10.d0^(esr +nd*(findgen(nb)+2.)/nb) ; bin lower radii r = (ra + rb) / 2.0 ; bin center radii, microns wid = rb - ra ; bin widths, microns ;print,'min,max ',min(r),max(r) ;- calculate differential concentration (dconc) in each radius bin dcon = (N / (sqrt(2.*!pi) * dr)) * exp( -1.*(r - rm)^2. / (2.*dr^2.) ) con = dcon * wid ; #/cm-3 ;- trim down the radius range to only those with con > max(con) * cutf, ; this is usefull when doing integrals of optical cross section numbins = nb if (cutf gt 0) then begin ok = where(con gt max(con) * cutf , numbins) dcon = dcon(ok) con = con(ok) r = r(ok) wid = wid(ok) endif ;- scale the concentration vs. radius so that the integral = N ; this is required when the median radius approaches the width and ; some of the concentration would have been out at negative radii con = con * N / total(con) dcon = con / wid ;- compute surface area, volume density, and effective radius vol = total( ((4./3.) * !pi * r^3 ) * con) sfc = total( (4.0 * !pi * r^2 ) * con) re = 3*vol/sfc ;- done return end