function [beta,alpha,LR,model,ab]=fernald_new1(z,pr,ref,Lambda,LidarRatio,Pair,Tair)
% -----------------------------------------------------------------------
% Description:
%
%    Defines important physic's constants, most of which are used
%    by the molecular scattering calculation.
%
%    Then takes temperature and pressure from the reference sounding and
%    calculates the refractive index, king's correction factor,
%    depolarization ratio, molecular phase function and cross-section
%    and, finally, molecular backscatter and extinction coefficients.
%
% Usage:
%
%    [alpha_mol,beta_mol,LR_mol,Constants]=alphabeta(pres,temp,lambda,co2ppmv);
%    Input the desired pressure, temperature, laser wavelength and
%    CO2 concentration and execute the M-file.
%
% Input:
%
%    z              - Height Profile                                    [m]
%    pr             - Lidar Profile                                     [-]
%    ref            - Reference Height(s)                               [m]
%    Lambda         - Lidar Wavelength                                  [m]
%    LidarRatio     - Extinction to Backscatter Ratio                  [sr]
%    Pair           - Atmospheric pressure profile                      [P]
%    Tair           - Atmospheric temperature profile                   [K]
%
% Ouput:
%
%    beta           - Backscatter Structure with fields:
%                     - mol: molecular                               [m^-1] 
%                     - aer: aerosols                                [m^-1]
%                     - tot: total                                   [m^-1]
%    alpha          - Extinction Structure with fields
%                     - mol: molecular                         [m^-1 sr^-1] 
%                     - aer: aerosols                          [m^-1 sr^-1]
%                     - tot: total                             [m^-1 sr^-1]
%    LR             - Lidar Ratio Structure with Fields
%                     - mol: molecular                                 [sr] 
%                     - aer: aerosols                                  [sr]
%    model          - Pure Molecular Lidar Signal from Modelled Extinction
%                       and Backscatter
%    ab             - Fitting Constants 
%                     - a(1,:)      
%                     - a(2,:) 
%
% Authors:
%
%    Pablo Ristori    (pablo.ristori@gmail.com) CEILAP, UNIDEF (CITEDEF-CONICET), Argentina
%    Lidia Otero      (lidia1116@gmail.com)     CEILAP, UNIDEF (CITEDEF-CONICET), Argentina
%
% ------------------------------------------------------------------------
[beta.mol,beta.aer,beta.tot,alpha.mol,alpha.aer,alpha.tot]=...
    deal(zeros(size(pr)));
ab = zeros(2,size(pr,2));

ref=round(ref/(z(2)-z(1)));
if length(ref)>1
    reference=(ref(1):ref(2))';
end
for k=1:size(pr,2)
    if (size(Pair,2)>1)||(k==1)
        [alpha_mol,beta_mol,LR_mol]=alphabeta(Pair(:,k),Tair(:,k),Lambda);
    end
    LR.mol(:,k)=LR_mol;
    LR.aer(:,k)=LidarRatio;
    
    if length(ref)==2 % fit by region
        model = (beta_mol.*...
            exp(-2*cumtrapz(z,alpha_mol)))./z.^2;
        meanmodel = mean(model(reference));
        meanpr = mean(pr(reference,k));
        model=model/meanmodel;
        pr=pr/meanpr;
        ab(:,k) = [model(reference,k) ones(size(reference))]\pr(reference,k);
        pr(:,k) = (pr(:,k)-ab(2,k))/ab(1,k);
        ab(:,k) = [ab(1,k)*meanpr/meanmodel ab(2,k)*meanpr];
        betaref=beta_mol(ref(1))*(pr(ref(1),k)/model(ref(1)));
    else
        betaref=beta_mol(ref(1));
    end
    s = pr.*z.^2;
    pp = exp(-2*cumtrapz(z,(LR.aer-LR.mol)*beta_mol));
    spp = s .* pp;
    sppr = spp ./ spp(ref(1));    
    beta.mol(:,k)  = beta_mol;
    beta.tot(:,k)  = sppr ./ (1/betaref - ...
        (cumtrapz(z,2*LR.aer*sppr)-trapz(z(1:ref(1)),2*LR.aer*sppr(1:ref(1)))));
    beta.aer(:,k)  = beta.tot(:,k) - beta.mol(:,k);
    alpha.mol(:,k) = alpha_mol(:,k);
    alpha.aer(:,k) = beta.aer(:,k) .* LR.aer(:,k);
    alpha.tot(:,k) = alpha.mol(:,k) + alpha.aer(:,k);
end
end

function [alpha_mol,beta_mol,LR_mol,Constants]=alphabeta(pres,temp,lambda,co2ppmv)
%
% -----------------------------------------------------------------------
% Description:
%
%    Defines important physic's constants, most of which are used
%    by the molecular scattering calculation.
%
%    Then takes temperature and pressure from the reference sounding and
%    calculates the refractive index, king's correction factor,
%    depolarization ratio, molecular phase function and cross-section
%    and, finally, molecular backscatter and extinction coefficients.
%
% Usage:
%
%    [alpha_mol,beta_mol,LR_mol,Constants]=alphabeta(pres,temp,lambda,co2ppmv);
%    Input the desired pressure, temperature, laser wavelength and
%    CO2 concentration and execute the M-file.
%
% Input:
%
%    pres           - Atmospheric pressure profile                      [P]
%    temp           - Atmospheric temperature profile                   [K]
%    lambda         - wavelength vector                                 [m]
%    co2ppmv        - CO2 concentration                              [ppmv]
%
% Ouput:
%
%    alpha_mol        - extinction coefficient                       [m^-1]
%    beta_mol         - backscatter coefficient                [m^-1 sr^-1]
%    Constants.k      - Boltzmann constant                 [m^2 kg / s^2 K]
%    Constants.Na     - Avogadro's number                               [#]
%    Constants.Mvol   - Molar volume at 1013.25hPa and 273.15K   [L mol^-1]
%    Constants.Rgas   - Universal gas constant              [J K^-1 mol^-1]
%    Constants.Rair   - Dry air gas constant                 [J K^-1 kg^-1]
%    Constants.Tstd   - Temperature of standard atmosphere              [K]
%    Constants.Pstd   - Pressure of standard atmosphere                [Pa]
%    Constants.Nstd   - Molecular density at Tstd and Pstd         [# m^-3]
%
%    Atmospheric Constituents Concentration [ppv] and
%    Molecular weight for: N2, O2, Ar, Ne, He, Kr, H2 Xe         [g mol^-1]
%
%    Constants.lambda - wavelengths                                     [m]
%    Constants.nAir   - index of refraction                             [-]
%    Constants.fAir   - King factor of air                              [-]
%    Constants.rhoAir - depolarization factor of air                    [-]
%    Constants.Pf_mol - phase function at -180deg                       [-]
%    Constants.LR_mol - molecular lidar ratio                          [sr]
%
% Authors:
%
%    H. M. J. Barbosa (hbarbosa@if.usp.br), IF, USP, Brazil
%    B. Barja         (bbarja@gmail.com), GOAC, CMC, Cuba
%    R. Costa         (re.dacosta@gmail.com), IPEN, Brazil
%
% References:
%
%    Seinfeld and Pandis (1998)
%    Handbook of Physics and Chemistry (CRC 1997)
%    Wallace and Hobbs, Atmospheric science (2nd Edition)
%    Bodhaine et al, 1999: J. Atmos. Ocea. Tech, v. 16, p.1854
%
%    Bates, 1984: Planet. Space Sci., v. 32, p. 785
%    Bodhaine et al, 1999: J. Atmos. Ocea. Tech, v. 16, p.1854
%    Bucholtz, 1995: App. Opt., v. 34 (15), p. 2765
%    Chandrasekhar, Radiative Transfer (Dover, New York, 1960)
%    Edlen, 1953: J. Opt. Soc. Amer., v. 43, p. 339
%    McCartney, E. J., 1976: Optics of the Atmosphere. Wiley, 408 pp.
%    Peck and Reeder, 1973: J. Opt. Soc. Ame., v 62 (8), p. 958
%
% ------------------------------------------------------------------------
;
%% Input options and conversions
if nargin == 3
    co2ppmv = 392;                               % CO2 concentration [ppmv]
end

%% -----------------------------------------------------------------------
% Fixed definitions
%------------------------------------------------------------------------

% Atmospheric Constituents Concentration                              [ppv]
% ref: Seinfeld and Pandis (1998)
N2ppv   =    0.78084;
O2ppv   =    0.20946;
Arppv   =    0.00934;
Neppv   =    1.80    * 1e-5;
Heppv   =    5.20    * 1e-6;
Krppv   =    1.10    * 1e-6;
H2ppv   =    5.80    * 1e-7;
Xeppv   =    9.00    * 1e-8;
CO2ppv  = co2ppmv    * 1e-6;

% Atmospheric Constituents Molecular weight [g mol^-1]
% ref: Handbook of Physics and Chemistry (CRC 1997)
N2mwt  =  28.013;
O2mwt  =  31.999;
Armwt  =  39.948;
Nemwt  =  20.18;
Hemwt  =   4.003;
Krmwt  =  83.8;
H2mwt  =   2.016;
Xemwt  = 131.29;
CO2mwt =  44.01;

% Dry air molecular mass [g mol^-1]
Airmwt= (...
    N2ppv*N2mwt + O2ppv*O2mwt + Arppv*Armwt + Neppv*Nemwt + ...
    Heppv*Hemwt + Krppv*Krmwt + H2ppv*H2mwt + Xeppv*Xemwt + ...
    CO2ppv*CO2mwt) / (...
    N2ppv + O2ppv + Arppv + Neppv + Heppv + Krppv + H2ppv + Xeppv + CO2ppv);

% Physic Constants
k  = 1.3806503e-23; % Boltzmann Constant                           [J K^-1]
Na = 6.0221367e+23;  % Avogadro's number                         [# mol^-1]

% Wallace and Hobbs, p. 65
Rgas = k*Na;            % Universal gas constant            [J K^-1 mol^-1]
Rair = Rgas/Airmwt*1e3; % Dry air gas constant               [J K^-1 kg^-1]

% Standard Atmosphere Reference Values
T0   =  273.15;  % zero deg celcius                                     [K]
Tstd =  288.15;  % Temperature                                          [K]
Pstd = 101325;   % Pressure                                            [Pa]

% Molar volume at Pstd and T0
% Bodhaine et al, 1999
Mvol =  22.4141e-3;                                          % [m^3 mol^-1]

% Molecular density at Tstd and Pstd
Nstd = (Na/Mvol)*(T0/Tstd);                                      % [# m^-3]

%% -----------------------------------------------------------------------
%  REFRACTIVE INDEX WITH CO2 CORRECTION
% ------------------------------------------------------------------------

% Calculate (n-1) at 300 ppmv CO2
% Peck and Reeder (1973), Eqs (2) and (3)
% or Bucholtz (1995), Eqs (4) and (5)
% or Bodhaine et al (1999), Eq (4)
if ((lambda*1e6) > 0.23)
    dn300 = (5791817 ./ (238.0185 - 1./(lambda*1e6).^2) + ...
        167909 ./ (57.362 - 1./(lambda*1e6).^2))*1e-8;
else
    dn300 = (8060.51 + 2480990 ./ (132.274 - 1./(lambda*1e6).^2) + ...
        14455.7 ./ (39.32957 - 1./(lambda*1e6).^2))*1e-8;
end

%  Correct for different concentration of CO2
%  Bodhaine et al (1999), Eq (19)
dnAir = dn300 .* (1 + (0.54 * (co2ppmv*1e-6 - 0.0003)));

%  Actual index of refraction at 300 ppmv CO2
nAir = 1 + dnAir;

%% -----------------------------------------------------------------------
%  KING FACTOR AND DEPOLARIZATION RATIO
% ------------------------------------------------------------------------

% Bates (1984), Eqs (13) and (14)
% or Bodhaine et al (1999), Eqs (5) and (6)

% Nitrogen
fN2  = 1.034 + ...
    ( 3.17e-4 ./ ((lambda*1e6).^2));
% Oxygen
fO2  = 1.096 + ...
    (1.385e-3 ./ ((lambda*1e6).^2)) + ...
    (1.448e-4 ./ ((lambda*1e6).^4));
% Argon
fAr  = 1;
% Carbon dioxide
fCO2 = 1.15;

% Bodhaine et al (1999) Eq (23)
% NOTE: numerator and denominator are not written in percent

% Standard dry air mixture with co2ppmv
fAir = ...
    (N2ppv*fN2 + O2ppv*fO2 + Arppv*fAr + co2ppmv*1e-6*fCO2)./...
    (N2ppv     + O2ppv     + Arppv     + co2ppmv*1e-6);
%
% Depolarization ratio estimated from King's factor
%
% hmjb - What is the correct way?
%rhoN2  = (6*fN2 -6)./(3+7*fN2 );
%rhoO2  = (6*fO2 -6)./(3+7*fO2 );
%rhoCO2 = (6*fCO2-6)./(3+7*fCO2);
%rhoAr  = (6*fAr -6)./(3+7*fAr );
%rhoAir = (0.78084*rhoN2+0.20946*rhoO2+0.00934*rhoAr+co2ppmv*1e-6*rhoCO2)/...
%	  (0.78084      +0.20946      +0.00934      +co2ppmv*1e-6)
%
% hmjb - What is the correct way?
rhoAir = (6*fAir-6)./(3+7*fAir);


%% -----------------------------------------------------------------------
%  RAYLEIGH PHASE FUNCTION AT 180deg
% ------------------------------------------------------------------------

% Chandrasekhar (Chap. 1, p. 49)
% or Bucholtz (1995), eqs (12) and (13)
gammaAir = rhoAir./(2-rhoAir);
Pf_mol = 0.75*((1+3.*gammaAir)+(1-gammaAir)*(cos(pi)^2))./(1+2.*gammaAir);

%% -----------------------------------------------------------------------
%  RAYLEIGH TOTAL SCATERING CROSS SECTION
%  -----------------------------------------------------------------------

% McCartney (1976)
sigma_std = 24 * (pi^3) * ((nAir.^2-1).^2) .* fAir ./...
    ((lambda.^4) .* (Nstd^2) .* (((nAir.^2)+2).^2));   %Cross Section [m^2]

%% -----------------------------------------------------------------------
%  RAYLEIGH VOLUME-SCATTERING COEFFICIENT
%  -----------------------------------------------------------------------

% In traditional lidar notation, Bucholtz (1995) eqs (2), (9) and (10)
% defines the scattering part of the molecular extinction coeficient.
% Therefore, here the usual greek letter 'alpha' is used instead of
% 'beta' as in Bucholtz.

% Bucholtz (1995), eq (9)
alpha_std = Nstd * sigma_std;      % extinction at std pres and temp [m^-1]

% Bucholtz (1995), eq (10), units [m^-1]
% scaling for each P and T in the column
alpha_mol = pres./temp*Tstd/Pstd * alpha_std; % molecular extinction [m^-1]

%% -----------------------------------------------------------------------
%  RAYLEIGH ANGULAR VOLUME-SCATTERING COEFFICIENT
% ------------------------------------------------------------------------

% Rayleigh extinction to backscatter ratio
LR_mol = (4*pi)./Pf_mol;                         %Rayleigh lidar ratio [sr]

% In traditional lidar notation, Bucholtz (1995) eq (14) defines the
% backscattering coeficient. Here the usual greek letter 'beta' is
% used as in Bucholtz.

% Multiply by phase function for -180deg and divide by 4pi steradians
% Units: [m]-1 [sr]-1

beta_mol = alpha_mol*diag(LR_mol.^-1);       % molecular backscatter [m^-1]

%% -----------------------------------------------------------------------
%  CREATE VARIABLE "Constants" AS A STRUCTURE
% ------------------------------------------------------------------------
variables=who;
for k=1:length(variables);
    if exist(variables{k},'var')
        eval(['Constants.(''' variables{k} ''')=' variables{k} ';']);
    end
end
Fields={'k','alpha_mol','beta_mol','LR_mol','pres','temp'};
Constants = rmfield(Constants, Fields);
end