LABELA-NET

LABELA-NET: Lidar Algorithm of Backscatter and Extinction for LALINET

Based on developed algorithms: 
LEBEAR- Lidar Elastic Backscatter and Extinction Analysis Routine, by: Fábio Lopes - USP/IPEN-Brazil
ROCEBSELI - Rutina paraObtener los Coeficientes de Extinctión y Backscatter a partir de Señal Elástica LIdar, by: Daniel Nisperuza - UNAL/GLEO-Colombia

Initial set up

SetDirectory[NotebookDirectory[]];

version="v01";

λ=0.532; (*in micrometers*)
lidarfiletype=".dat";
radiosoundfilename="**********.***"; (*radiosounding name file.type*)
fileres=7.5 ;(*Lidar file resolution in meters*)
r0=7.5;(*First altitude of lidar signal measured*)

Unit factor: Choose the final unit of the Extinction and Backscatter coefficients  
1 for (m^(-1))
1000 for (km^(-1))
(10^6) for (km^(-1))

unitfactor=1000;

Atmospheric information - radiosounding data

Physics constants

<<PhysicalConstants`
<<StandardAtmosphere`
<<PlotLegends`
<<Units`

molarmass=28.9644 (10^(-3)); (*molecular mass for dry air (kg (m^(-3)))*)
Avogadro=6.02214179`*^23;
Rm=MolarGasConstant Kelvin Mole/Joule; (* J/K/mole *)
tempstandard=N[ConvertTemperature[15,Celsius,Kelvin]];
CO2=392;

Radiosounding data importing

fileradiosound=Import[FileNames[radiosoundfilename][[1]],"Table"];

pressurefile=Table[{fileradiosound[[i]][[2]]-fileradiosound[[1]][[2]]+r0,fileradiosound[[i]][[1]]},{i,1,Length[fileradiosound]}];
temperaturefile=N[Table[{fileradiosound[[i]][[2]]-fileradiosound[[1]][[2]]+r0,ConvertTemperature[fileradiosound[[i]][[3]],Celsius,Kelvin]},{i,1,Length[fileradiosound]}]];
presshPa=Interpolation[pressurefile];
temperaturek=Interpolation[temperaturefile];
z0=temperaturefile[[1]][[1]];
maxrangerayleigh=temperaturefile[[-1]][[1]];

Plot: Atmospheric temperature and pressure - radiosounding profiles

temperaturegraph=ListPlot[temperaturefile,Frame->True,FrameLabel->{"ALtitude (m)","Temperature (K)"},FrameStyle->Directive[{Thickness[1/300]}],BaseStyle->{FontFamily->"Tahoma", FontSize->18,FontWeight->"Bold"}, PlotLabel->"Temperature Profile",PlotStyle->{Darker[Green],Thickness[0.004]},GridLines->Automatic, GridLinesStyle->Directive[White,Dashed],PlotRange->{All},
ImageSize->900];

pressuregraph=ListPlot[pressurefile,Frame->True,FrameLabel->{"ALtitude (m)","Pressure (hPa)"},FrameStyle->Directive[{Thickness[1/300]}],BaseStyle->{FontFamily->"Tahoma", FontSize->18,FontWeight->"Bold"}, PlotLabel->"Pressure Profile",PlotStyle->{Darker[Green],Thickness[0.004]},GridLines->Automatic, GridLinesStyle->Directive[White,Dashed],PlotRange->{All},
ImageSize->900];

Refractive index - assuming the (Subscript[CO, 2])global concentration

(* Assuming the (Subscript[CO, 2]) Concentration of 392 ppm according to global monthly mean measured at Mauna Loa Observatory,Hawaii By NOAA - http://www.esrl.noaa.gov/gmd/ccgg/trends/global.html*)

Real refractive index for a standard temperature and pressure condictions (288.15 K, 1013.25 mb) - λ > 0.23 μm

MrealCO2[λ_,CO2PPM_]:=1+(10^(-8))((5791817/(238.0185-((((1/λ)))^2)))+(167909/(57.362-((((1/λ)))^2))))(1+0.54(CO2PPM (10^(-6))-0.0003));

King Factor

The King factor is used in order to correct dispersion and scattering of light due the effects of molecular anisotropy

kingN2[λ_]:=1.034+3.17 (10^(-4))((((1/λ)))^2)

kingO2[λ_]:=1.096+1.385 (10^(-3))((((1/λ)))^2)+1.448 (10^(-4))((((1/λ)))^4)

For the pure air the King factor is equal to 1. The king factor for (Subscript[CO, 2])is equal to 1.15 
 However, it can be calculate the King Factor for air with varying concentrations of (Subscript[CO, 2])

kingAir[λ_]:=1(λ^0)

kingCO2[λ_]:=1.15 (λ^0)

Since the king factor is expressed as concentration by volume in percent,     
   it is written as a factor of 100 larger than the value used in the equation for the refractive index -> ppm:(10^(-6))is written as (10^(-4))

kingAirCO2[λ_,CO2PPM_]:=(((78.084 kingN2[λ]+20.946 kingO2[λ]+0.934kingAir[λ]+CO2PPM (10^(-4))kingCO2[λ]))/((78.084+20.946+0.934+CO2PPM (10^(-4)))))

Depolarization of the air

ρ[F_]:=((6F-6)/(3+7F))

ρn[λ_,CO2PPM_]:=(((78.084ρ[ kingN2[λ]]+20.946ρ[kingO2[λ]] +0.934ρ[kingAir[λ]]+CO2PPM (10^(-4))ρ[kingCO2[λ]]))/((78.084+20.946+0.934+CO2PPM (10^(-4)))))

Pdepol[θ_,λ_,CO2PPM_]:=(3/(4(1+2((ρn[λ,CO2PPM])/(2-ρn[λ,CO2PPM]))))) ((1+3((ρn[λ,CO2PPM])/(2-ρn[λ,CO2PPM])))+(1-((ρn[λ,CO2PPM])/(2-ρn[λ,CO2PPM])))((Cos[θ])^2))

(Subscript[σ, ray]) - Rayleigh scattering cross section - T and P standard

airdensitynumber[Avogadro_,press_,temp_,Rm_,Cx_]:=( ((Avogadro press)/(Rm temp)))Cx      (*unit in (m^3)*)

(Subscript[σ, ray])[λ_,CO2PPM_]:=((24 (π^3)(((((MrealCO2[λ,CO2PPM])^2)-1))^2))/((((λ (10^(-6))))^4)(((airdensitynumber[Avogadro,Pressure[0](10^2)/(Bar Milli),tempstandard,Rm,(0.78+0.21+0.0093+CO2PPM (10^(-6)))] ))^2)(((((MrealCO2[λ,CO2PPM])^2)+2))^2)))kingAirCO2[λ,CO2PPM](*unit in (m^(-2))*)

(Subscript[α, mol]) - Rayleigh extinction coefficient - radiosouding

αstd=airdensitynumber[Avogadro,Pressure[0](10^2)/(Bar Milli) ,tempstandard,Rm,(0.78+0.21+0.0093+CO2 (10^(-6)))](Subscript[σ, ray])[λ,CO2];

αmol[z_]:=αstd (tempstandard/(Pressure[0]/(Bar Milli) ))((presshPa[z])/(temperaturek[z]));(*in (m^(-1))*)

αmolnorm=Table[αmol[i],{i,z0,maxrangerayleigh,fileres}];

αmolgraf2=Plot[ unitfactor αmol[z],{z,z0,maxrangerayleigh},Frame->True,FrameLabel->{"Altitude (m)","Rayleigh Scattering Coefficient (Subscript[α, s]) ((km^(-1)))"},FrameStyle->Directive[{Thickness[1/300]}],BaseStyle->{FontFamily->"Tahoma", FontSize->16,FontWeight->"Bold"}, PlotLabel->"Rayleigh Extinction Coefficient",PlotStyle->{Darker[Green],Thickness[0.004]},GridLines->Automatic, GridLinesStyle->Directive[Gray,Dashed],PlotRange->All,ImageSize->700];

(Subscript[β, mol,θ]) - Rayleigh backscattering coefficient - using depolarization function

βmolθ[z_]:=((αmol[z])/(4π))Pdepol[π,λ,CO2]

betamolecularexport=Table[βmolθ[i],{i,z0,maxrangerayleigh,fileres}];

βmolθgraf=Plot[unitfactor βmolθ[z],{z,z0,maxrangerayleigh},Frame->True,FrameLabel->{"Altitude (m)","Rayleigh Backcattering Coefficient (Subscript[β, s]) (sr (km^(-1)))"},FrameStyle->Directive[{Thickness[1/300]}],BaseStyle->{FontFamily->"Tahoma", FontSize->16,FontWeight->"Bold"}, PlotLabel->"Rayleigh Backscattering Coefficient",PlotStyle->{Darker[Red],Thickness[0.004]},GridLines->Automatic, GridLinesStyle->Directive[Gray,Dashed],PlotRange->All,ImageSize->700];

(Subscript[S, mol]) - Rayleigh extinction-to-backscatter ratio

Smol=((αstd)/(((αstd)/(4π))Pdepol[π,λ,CO2]));(*unit in sr*)

Range corrected lidar signal

molecularrminimum=12000;(*minimum range for molecular fitting - extracting the background value*)
maxrangebg=15000; (*maximum range to apply into background calculation*)

Reading lidar file

signalfile=Flatten[Import[FileNames[ToString["*"<>StringDrop[ToString[λ],2]<>"*"<>lidarfiletype]][[1]],"Table"]];
r=Table[i fileres,{i,1,Length[signalfile]}];(*Lidar data trange using vertical space resolution and first altitude*)

Molecular fitting

integralmol=Table[SetPrecision[fileres*((αmolnorm[[1]])/2)+fileres*((αmolnorm[[i]])/2)+fileres*((∑_(j=1)^(i-1))αmolnorm[[1+j]]),6],{i,1,Length[αmolnorm]}];
AOTmol=Table[Exp[-2 integralmol[[i]]],{i,1,Length@integralmol}];
smolsim=Table[((betamolecularexport[[i]]AOTmol[[i]])/((r[[i]])^2)),{i,1,IntegerPart[(maxrangebg/fileres)]}];
noise=Table[{smolsim[[i]],signalfile[[i]]},{i,IntegerPart[(molecularrminimum/fileres)],IntegerPart[(maxrangebg/fileres)]}];

noisesignalgraph=ListLinePlot[noise,Frame->True,FrameLabel->{"Synthetic Signal","Measured lidar signal (A.U.)"},FrameStyle->Directive[{Thickness[1/300]}],BaseStyle->{FontFamily->"Tahoma", FontSize->13,FontWeight->"Bold"}, PlotLabel->"Lidar signal vs. Synthetic signal -molecular region",PlotStyle->{Darker[Blue],Thickness[0.005]},GridLines->Automatic, GridLinesStyle->Directive[Gray,Dashed],PlotRange->All,AspectRatio->0.7,ImageSize->600];

noisefitting =SetPrecision[LinearModelFit[noise, {1, x}, x],10]; noisefitting[{"BestFit", "ParameterTable"}]

bg=SetPrecision[ToExpression[noisefitting["BestFit"][[1]]],12];
scale=ToExpression[noisefitting["BestFit"][[2]][[1]]];
simulatedsignalscaled=Table[{r[[i]],scale smolsim[[i]]},{i,1,Length[smolsim]}];

Background Extraction and Range Corrected Signal

bgextractedsignal=Table[{r[[i]],signalfile[[i]]-bg},{i,1,Length[r]}];
RCS=Table[{fileres i,(signalfile[[i]]-bg)(((fileres i))^2)},{i,1,Length[signalfile]}];

bgsimulatedsignalgraph=ListLogPlot[{bgextractedsignal,simulatedsignalscaled},Frame->True,FrameLabel->{"Altitude (m)","Offset corrected signal (U.A.)"},FrameStyle->Directive[{Thickness[1/300]}],BaseStyle->{FontFamily->"Tahoma", FontSize->18,FontWeight->"Bold"}, PlotLabel->"Lidar signal and Synthetic signal",PlotStyle->{{Darker[Blue],Thickness[0.006]},{Darker[Red],Thickness[0.004]}},GridLines->Automatic, GridLinesStyle->Directive[Lighter[Gray],Dashed],Joined->True,PlotRange->{{0,maxrangebg},All},ImageSize->700];

rangecorrectedgraf=ListLogPlot[RCS,Frame->True,FrameLabel->{"Altitude (m)","Range corrected signal (A.U.)"},FrameStyle->Directive[{Thickness[1/300]}],BaseStyle->{FontFamily->"Tahoma", FontSize->18,FontWeight->"Bold"},Joined->True, PlotLabel->"Backscatter Range Corrected Signal - BG extracted",PlotStyle->{Darker[Blue],Thickness[0.001]},GridLines->Automatic, GridLinesStyle->Directive[Lighter[Gray],Dashed],PlotRange->{{0,maxrangebg},All},ImageSize->700];

Aerosol backscatter and extinction retrieval

Backscatter and Extinction coefficients retrieval

Saer=39;
βaerzo=0; (*backscatter reference value for molecular range*)
rmax=15000; (*maximum range to retrieval backscatter and extinction profile*) 
rrefinitial=12000; (*minimum altitude of reference range*)
rreffinal=14000; (*maximum altitude of reference range*)
rminklett=180; (*minimum range to start the calculation of AOD*) 
rmaxklett=8000; (*maximum range to start the calculation of AOD*)

Backscatter and extinction calculation

auxmin=Ceiling[(r0/fileres)];
auxmax=Ceiling[(rmax/fileres)];
auxrefinitial=Ceiling[(rrefinitial/fileres)];
auxreffinal=Ceiling[(rreffinal/fileres)];
calibration=Mean[Table[((RCS[[i]][[2]])/((βaerzo+betamolecularexport[[i]]))),{i,auxrefinitial,auxreffinal}]];
βmolθ1=Table[SetPrecision[(Saer-Smol)*βmolθ[r[[i]]],6],{i,1,auxreffinal}];
int=Table[SetPrecision[fileres*((βmolθ1[[i]])/2)+fileres*((βmolθ1[[Length[βmolθ1]]])/2)+fileres*((∑_(j=i)^(Length@βmolθ1-1))βmolθ1[[1+j]]),6],{i,1,auxreffinal}];
Trans=Table[SetPrecision[Exp[2int[[i]]],6],{i,1,auxreffinal}];Tnorm=Table[SetPrecision[Trans[[i]]*RCS[[i,2]],6],{i,1,auxreffinal}];
integ=Table[SetPrecision[fileres*Saer*((Tnorm[[i]])/2)+fileres*Saer*((Tnorm[[Length[Tnorm]]])/2)+fileres*Saer*((∑_(j=i)^(Length@Tnorm-1))Tnorm[[1+j]]),6],{i,1,auxreffinal}];
betacoeff=Table[{SetPrecision[(((Tnorm[[i]])/(calibration+2*integ[[i]]))-βmolθ[r[[i]]]),6],r[[i]]},{i,1,auxreffinal}];
 alphacoeff=Table[{SetPrecision[Saer betacoeff[[i]][[1]],6],r[[i]]},{i,1,auxreffinal}];
betacoefftotal=Table[{SetPrecision[betacoeff[[i]][[1]]+( βmolθ[fileres i]),6],r[[i]]},{i,1,auxreffinal}];
betamolgraph=Table[{SetPrecision[ βmolθ[fileres i],6],r[[i]]},{i,1,auxreffinal}];
particulateresults=Table[{r[[i]],unitfactor betacoeff[[i]][[1]],unitfactor alphacoeff[[i]][[1]],Saer},{i,1,Length[betacoeff]}];
finalresults=Prepend[particulateresults,{"Altitude (m)","Backscatter (sr (km^(-1)))","Extinction ((km^(-1)))","Lidar ratio (sr)"}];

TableForm[finalresults,TableHeadings->{None,{"Aerosol optical properties from elastic lidar signal"<>" "<>StringDrop[ToString[λ],2]<>" nm","",DateString[]}}];

Particulate backscatter and extinction profiles

backscattergraph=ListLinePlot[MovingAverage[unitfactor betacoeff,1],Frame->True,FrameLabel->{"Backscatter Coefficient ((km^(-1))(sr^(-1)))","Altitude (m)"},FrameStyle->Directive[{Thickness[1/300]}],BaseStyle->{FontFamily->"Tahoma", FontSize->16,FontWeight->"Bold"}, PlotLabel->"Particulate Backscatter Profile",PlotStyle->{Darker[Green],Thickness[0.007]},GridLines->Automatic, GridLinesStyle->Directive[LightGray,Dashed],PlotRange->All,
ImageSize->500,AspectRatio->1.3]

totalbackscattergraph=ListLinePlot[{unitfactor betacoefftotal,unitfactor betamolgraph,1},Frame->True,FrameLabel->{"Backscatter coefficient ((km^(-1))(sr^(-1)))","Altitude (m)"},FrameStyle->Directive[{Thickness[1/300]}],BaseStyle->{FontFamily->"Tahoma", FontSize->16,FontWeight->"Bold"}, PlotLabel->"Total backscatter profile",PlotStyle->{{Darker[Green],Thickness[0.009]},{Darker[Red],Thickness[0.004]}},GridLines->Automatic, GridLinesStyle->Directive[LightGray,Dashed],PlotRange->All,
ImageSize->500,AspectRatio->1.3]

extinctiongraph=ListLinePlot[unitfactor alphacoeff,Frame->True,FrameLabel->{"Extinction Coefficient ((km^(-1)))","Altitude (m)"},FrameStyle->Directive[{Thickness[1/300]}],BaseStyle->{FontFamily->"Tahoma", FontSize->16,FontWeight->"Bold"}, PlotLabel->"Particulate Extinction Profile",PlotStyle->{Darker[Blue],Thickness[0.006]},GridLines->Automatic, GridLinesStyle->Directive[LightGray,Dashed],PlotRange->All,
ImageSize->500,AspectRatio->1.3]

Aerosol optical depth calculation

Interalpha:=Interpolation[Table[{r[[i]],alphacoeff[[i]][[1]]},{i,1,IntegerPart[(rmaxklett/fileres)]}],InterpolationOrder->1];
Intalpha[ii_,iii_]:=Integrate[Interalpha[i],{i,r[[ii]],r[[iii]]}];
auxminklett=IntegerPart[(rminklett/fileres)];
auxmaxklett=IntegerPart[(rmaxklett/fileres)];
AODlidar=Intalpha[auxminklett,auxmaxklett]/(1 (10^3));

Exporting profile files

Export["output_01_"<>StringTake[ToString[λ],-3]<>"nm"<>"_"<>version<>".dat",results];
resp=Join[{{"LR","overlap","ref_initial", "ref_final" ,"calibration_coeff","BG","AOD","AOD initial altitude","AOD final altitude"}},{{Saer,rrefinitial,rreffinal,calibration,bg,AODlidar, rminklett,rmaxklett}}];
Export["output_02_"<>StringTake[ToString[λ],-3]<>"nm"<>"_"<>version<>".dat",resp];
Export["backscatter_profile_"<>StringTake[ToString[λ],-3]<>"nm_"<>version<>".tiff",backscattergraph];
Export["extinction_profile_"<>StringTake[ToString[λ],-3]<>"nm_"<>version<>".tiff",extinctiongraph];
Export["total_backscatter_profile_"<>StringTake[ToString[λ],-3]<>"nm_"<>version<>".tiff",totalbackscattergraph];