function [licorval,h2ommol] = licorWebbPTK(co2air,h2oair,pair,tair,RH,tsonic,fs)

%% Calculate the concentration of CO2 using the raw licor data
%% from the 7500.  Uses an approach that is similar to the closed
%% path sensors.  
%%
%% This version corrects for hydrostatic pressure changes using the 
%% mean Vaisala P + dp = P - rho g dz where dz is taken from the heave
%% measurement.  Note that the measured P' and the heave generated dp 
%% are very similar after lagging the Vaisala P by 3 seconds.  That agrees
%% very well with the response time given by Vaisals, i.e., a 1 second
%% average and a 2 second response time.
%%
%% The input is:
%%
%% co2air: the number density of CO2 in moist air (mmol/m^3); this is rho_c in
%%         Webb et al (1980) and c in McGillis et al. (2001)
%% h2oair: the number density of H2O in moist air (mmol/m^3)
%% tair: the air temperature in C from the vaisala or other reference T
%% RH: the relative humidity in % from the vaisala or other reference RH
%% pair: the air pressure in mb from the vaisala or other reference P
%% tsonic: the sonic virtual temperature
%%
%% fs is the sonic sampling frequency
%
%  NOTE: This program assumes that all variables have a common time base,
%  i.e., they are all "sampled" at fs Hz using, e.g., interp1.m

S1=length(tsonic);
t1=(1:S1)/S1;
V1=length(tair);
tV=(1:V1)./V1;
Tair=interp1(tV,tair,t1,'spline','extrap');
%Pair=interp1(tV,pair,t1,'spline','extrap');
Pair=pair;
Exx=6.1121*exp(17.502*Tair./(240.97+Tair));
Exx=Exx.*(1.0007+Pair*3.46E-6);
Qsat=.622*Exx./(Pair-.375*Exx)*1000;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Design a filter to create Tvirt time series
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fc=1/60;                          %60 second cutoff 
nfreq = fs/2.0;   
wp = fc/nfreq;
ws = .7*wp;   
[n,wn] = buttord(wp,ws,10,25);
[bhigh,ahigh] = butter(n,wn,'high');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Determine the median values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rh = median(RH)/100;                       %Relative humidity as a fraction 
P0 = median(pair);                         %Atmospheric pressure (mb)
meanP=P0/10;                               %Atmospheric pressure (kPa)
flucP=pair./10;                            %Includes hydrostatic (kPa)
Press=P0*100;                              %Atmospheric pressure (Pa)
Tairk=median(tair)+273.16;                 %Air temperature (K)
Esat=BUCKESAT(median(tair),P0);            %Saturation vapor pressure
Ems=Esat*rh;                               %Vapor pressure
Qms=.622*Ems/(P0-.378*Ems)*1000;           %Specific humidity (g/Kg)
Tvirt=Tairk*(1+0.61*Qms/1000);             %Virtual temperature
rhoair = Press/287/Tairk;                  %Density of dry air
h2ommol=29000*meanP*Qms/18/8.314/Tvirt;    %mmol/m3
Rms = Qms./(1-Qms/1000);                   %Mixing ratio g/Kg

% The following equation computes the mole fractions (or mass concentrations) using 
% expressions given by the Licor 7500 manual pages 2-9 and 2-10.  The problem with 
% these expressions is that they are appropriate for number density computed in dry 
% air, while the number densities output by the Licor are the moles per unit volume
% of moist air.  
%
%For future reference, note that the following steps are equivalent
%
%q = 18*h2oair*8.314/29*Tvirt/meanP/1000;      %Specific Humidity g/Kg
%
%or
%
%wcell = h2oair*8.314*Tairk/meanP/1000;     %H2O mass concentration mmol/mol
%M = 29 - 11/1000*wcell;                    %Molecular weight of moist air
%q = 18*wcell./M;                           %Specific humidity g/Kg
%
%M*q = 18*wcell;
%29*q - 11/1000*wcell*q = 18*wcell
%29*q = 11/1000*wcell*q + 18*wcell
%29*q = 18*wcell*(11/18*q/1000 + 1)
%29*q = 18*wcell*(1 + 0.61*q/1000)
%29*q = 18*h2oair*8.314*Tairk/meanP/1000*(1 + 0.61*q/1000)
%29*q = 18*h2oair*8.314/meanP/1000*Tvirt
%q = 18*h2oair*8.314/meanP/1000*Tvirt/29
%
%
%Now let's check
%
%co2 = (44*8.314/29/meanP)*co2air.*Tvirt;  % Specific content of CO2
%
%versus
%
%ccell = co2air*8.314*Tairk/meanP;     %CO2 mass concentration mmol/mol
%co2 = 44*ccell./M;                         % Specific content of CO2
%co2 = 44*co2air*8.314*Tairk/meanP/M;
%co2 = 44*co2air*8.314*Tairk/meanP/29/(1 - 11/29/1000*wcell);
%co2 = 44*co2air*8.314*Tairk/meanP/29/(1 - 11/29/1000*q*M/18);
%co2 = 44*co2air*8.314*Tairk/meanP/29/(1 - 11/29/1000*q/18*(29-11*wcell/1000));
%co2 = 44*co2air*8.314*Tairk/meanP/29/(1 - 11/18*q/1000+11/29/1000*q/18*11*wcell/1000));
%co2 ~ 44*co2air*8.314*Tairk/meanP/29/(1 - 0.61*q/1000));
%co2 ~ 44*co2air*8.314*Tairk/meanP/29*(1 + 0.61*q/1000));
%co2 ~ (44*8.314/29/meanP)*co2air.*Tvirt;  % Specific content of CO2
%
%Need to check Chris's email to see which one is more exact.
%
%First approximation

%h2oair=h2ommol+h2oair-mean(h2oair);            %Use mean from Vaisala
q = 18*h2oair*8.314/29*Tvirt/flucP/1000;        %Specific Humidity g/Kg from Licor
                                                %with other means from
                                                %Vaisala

%Use sonic temperature to get better Tvirt
clear Tairk

Tairk=Tair+273.16;                             %Vaisala or other reference temp

Tvirt=tsonic-mean(tsonic)+mean(Tvirt);         %Let's keep it simple

qair = (18*8.314/29/1000)*h2oair.*Tvirt./flucP;%Recompute Specific Humidity g/Kg
                                               %from the Licor, using Tvirt
                                               %from the scaled sonic Ts.
rhfast=qair./Qsat;                             %fast rh from Licor and Vaisala
rair = qair./(1-qair/1000);                    %Mixing ratio g/Kg
rhov = rhoair*rair;                            %Corrected density

Tair = Tvirt./(1+0.61*qair./1000);             %Fluctuating air temperature

%***********************************************
% PKT Approach to remove cross talk
% Webb correct and then detrend
%***********************************************
co2 = (44*8.314/29)*co2air.*Tair.*(1+29/18.*rair/1000)./flucP;  %mg/kg

useQ=1;
if ~isnan(co2)
    co2save=co2;
    if useQ
        p1=polyfit(qair,co2,2);                        %Remove mean trend
        co2=mean(co2) + (co2-polyval(p1,qair));
    else
        p1=polyfit(rhfast,co2,2);                      %Remove mean trend
        co2=mean(co2) + (co2-polyval(p1,rhfast));
    end
end

co2air = rhoair * co2;            % This represents rho_a * r_c (mg/m^3),
                                  % which is what should be correlated
                                  % w'.
% rair = qair./(1-qair/1000);                    %Mixing ratio g/Kg
rco2 = co2./(1-qair/1000);                    %Mixing ratio mg/Kg
% rco2 = co2./(1-1.6*qair/1000);                   %done above I think 29/18=1.61

licorval = [rco2; rair; Tvirt; Tair;qair; rhfast];
