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'.
PARAMETERS
data
the raw data profile
hts
the raw heights profile
dsray
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
KEYWORDS
UNITS
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)
KRANGE  
Altitude range in km over which to compute the K density.
KHTS  
Set to a named variable in which to return the height array for the K profile.
RAYRANGE
Altitude range in km from which to compute the Rayleigh profile to subtract.
DENSERR  
Set to a named variable in which to return the error bars for the K densities.
DENSGOOD
Set to a named variable in which to returns 1 if the density calculation is good, 0 if not
History
Written originally in Spring 2002 as part of the larger Arecibo K lidar data processing
system.
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)

darr=DBLARR(nalt)

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

dwork=double(data)

 ;
 ; 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) ;

RETURN, darr END