Tools and Models:KLidarPhotons2Density
From CedarWiki
- 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
