Tools and Models:KLidarPhotons2Density

From CedarWiki
Jump to: navigation, search
FUNCTION KPhotons2Dens
This function takes a raw data profile from the K lidar and returns a density profile.
The density returned is in units of atoms/cc unless the "UNITS" keyword is set to 'm'.
the raw data profile
the raw heights profile
a data structure with the mean rayleigh profile to subtract from the data,
the Rayleigh and background subtraction range, and the 40 km atmospheric
density given by MSIS
set this to a string in order to specify the units of length to use for density
'm' means meters (atmos/m^3), 'cm' => centimeters (atoms/cc)
Altitude range in km over which to compute the K density.
Set to a named variable in which to return the height array for the K profile.
Altitude range in km from which to compute the Rayleigh profile to subtract.
Set to a named variable in which to return the error bars for the K densities.
Set to a named variable in which to returns 1 if the density calculation is good, 0 if not
Written originally in Spring 2002 as part of the larger Arecibo K lidar data processing
Modifications periodically, particularly in 2005 adding the MSIS atmospheric density profile
to remove the Rayleigh scattering profiles from the data.
This June 2011 version is a modification specifically for the CDB, with reduced options and
more automation.

Function KPhotons2Dens, $

           data, hts,  $
           KRANGE=krange, KHTS=khts, UNITS=units, $
           Denserr=denserr, DensGood=densgood

km = 1.0e3 cm = 1.0e-2

 ; MSIS-90 atmospheric number densities converted to atoms/m^3 and altitudes follow ...

msisdens = DOUBLE( $

          [  3.72723e+17,  3.18297e+17,  2.71909e+17,  2.32388e+17,  1.99005e+17,  1.70796e+17,  1.46938e+17,  1.26694e+17,  1.09494e+17, $
             9.48517e+16,  8.23568e+16,  7.16717e+16,  6.25160e+16,  5.46474e+16,  4.78763e+16,  4.20318e+16,  3.69744e+16,  3.25824e+16, $
             2.87583e+16,  2.54165e+16,  2.24867e+16,  1.99120e+16,  1.76438e+16,  1.56394e+16,  1.38657e+16,  1.22923e+16,  1.08947e+16, $
             9.65082e+15,  8.54226e+15,  7.55312e+15,  6.67004e+15,  5.88186e+15,  5.17836e+15,  4.55071e+15,  3.99166e+15,  3.49433e+15, $
             3.05269e+15,  2.66131e+15,  2.31520e+15,  2.00978e+15,  1.74087e+15,  1.50490e+15,  1.29824e+15,  1.11822e+15,  9.62444e+14, $
             8.28047e+14,  7.12039e+14,  6.11922e+14,  5.25509e+14,  4.50931e+14,  3.86593e+14,  3.31111e+14,  2.83277e+14,  2.42079e+14, $
             2.06606e+14,  1.76099e+14,  1.49875e+14,  1.27373e+14,  1.08078e+14,  9.15541e+13,  7.74231e+13,  6.53508e+13,  5.50368e+13, $
             4.62283e+13,  3.87163e+13,  3.23256e+13,  2.69040e+13,  2.23208e+13,  1.84613e+13,  1.52253e+13,  1.25213e+13,  1.02735e+13, $
             8.41695e+12,  6.89212e+12,  5.64516e+12,  4.62849e+12,  3.80111e+12,  3.12857e+12,  2.58200e+12,  2.13777e+12,  1.77654e+12, $
             1.48221e+12,  1.24198e+12,  1.04515e+12,  8.83492e+11,  7.50299e+11,  6.40232e+11,  5.49036e+11,  4.73267e+11,  4.10168e+11, $
             3.57520e+11,  3.13533e+11,  2.76745e+11,  2.46124e+11,  2.20071e+11,  1.97688e+11,  1.78333e+11,  1.61496e+11,  1.46767e+11, $
             1.33823e+11,  1.22387e+11  ] ) / cm^3

msisalt = findgen(N_ELEMENTS(msisdens)) + 30.0

 ; Interpolate the atmospheric density to match the altitudes of the lidar data

msisdenswork = INTERPOL(msisdens,msisalt,hts) ; interpolate to fit to the data range resolution IF (WHERE(msisdenswork LT 0))[0] NE -1 THEN msisdenswork[where(msisdenswork LT 0)] = 0.0 msisdenswork[WHERE(hts GT MAX(msisalt))] = 0.0

 ; rayrange is for the Rayleigh scatter. It starts above most aerosols and ends below the metal layer.

rayrange = [38,70]

 ; bkrange is limited to under 175 km because some data, mostly in 2009, 
 ; some in 2005-06 have odd behavior above that altitude

bkrange = [135,175] IF NOT KEYWORD_SET(KRANGE) THEN krange = [70,130] indkrange = WHERE(hts GE krange[0] and hts LE krange[1]) indray = WHERE(hts GE rayrange[0] and hts LE rayrange[1]) khts = hts[indkrange] nalt = N_ELEMENTS(khts) nhts = N_Elements(hts)


if KEYWORD_SET(units) then $ case STRLOWCASE(units) of 'm' : ufctr=1.0 'cm' : ufctr=cm^3 else : ufctr=cm^3 endcase $ else ufctr=cm^3

resback= DOUBLE(7.51e-17) ; m^2/ster, differential backscatter cross-section for K @ 200 K (von Zahn, 1996) lambda = 770.0  ; K lidar wavelength in nm rayback=5.45*(550.0^4/lambda^4)*1e-32 ; m^2/sr Rayleigh backscatter coefficient (Measures, p. 42)

 ; Remove Background using our mean Rayleigh profile with an MSIS density profiling through the metal range
 ; Establish the index range for the Rayleigh part and calculate the Rayleigh scaling factor
 ; create a work data array, remove the Rayleigh+noise background and range squared parts


 ; The following assumes that the data are "good". The user may want to put in some checking for 
 ; data quality. Some possible checks would be:
 ; 1. Mean background counts (> 130 km) under 10.
 ; 2. # of counts in the metal layer (minus background) > (600 shots) * (1 photon/shot) > 600
 ; First, fit a line to the background (data above K layer) and subtract that from the data

inds = where(hts ge bkrange[0] and hts le bkrange[1]) parms = linfit(hts[inds],dwork[inds]) ybgd = parms[0] + parms[1]*hts

 ; 1. Subtract noise background 

dwork = dwork - ybgd

 ; 2. Range correct the data

dwork = dwork*hts^2*km^2  ; x z^2 & convert hts from km to m

 ; 3. Subtract Rayleigh scaled to the background-subtracted, range-corrected data. 
 ;    The remainder is the metal layer.
 ;    SCLFCTR is the ratio of the range and rayleigh-corrected data to the atmospheric number density.
 ;    it is used to scale the data properly to a known quantity (the atmospheric density at 40 km). 

sclfctr = TOTAL(dwork[indray])/TOTAL(msisdenswork[indray]) dwork = dwork - msisdenswork*sclfctr darr = dwork[indkrange] / sclfctr * (rayback/resback) * ufctr

 ;  density error, make certain denominator >0, (data is longint)

denserr = darr/Sqrt(Double(data[indkrange])>0.01) ;