function [rhoc] = PKTcode(mixrat,rh,mixco2,W,rhoair)

% mixrat=Rv7500_UC(rng7200UC)/1000;  %kg/kg
% rh=rhfast7500_UC(rng7200UC)*100;  %(%)
% mixco2=Rc7500_UC(rng7200UC)./1e6;
% W=wj(rngw7200UC);

% mixrat=LI7500_H2O(rngw7200PSD)/1000;
% rh=rhfast7500_PSD(rngw7200PSD)*100;
% mixco2=LI7500_CO2mf(rngw7200PSD)/1e6;
% W=wj(rngw7200PSD);
 
mixco2nu=mixco2;  %initialization



mrh=nanmean(rh);
mrhoc=nanmean(mixco2);
mco2=44;
    
    
Cp=detrend(mixco2);
Cwc=nancov1(W,Cp);
wc=Cwc(1,2);   %m/s kg/kg
% cflux = wc * 1000 / mco2; % mols/(m2.s)
cflux = wc *1e6*rhoair/ mco2 / 1000; % mols/(m2.s)
cflux = cflux*60*60*24*365; % mols/(m2 yr)
% PKT method for correcting humidity cross sensitivity in
% open-path IR CO2 analysers,
% from Prytherch et al, 2010.
% Code by Peter K Taylor and John Prytherch, NOCS.
% Variables:
% mixrat = 20 hz H2O mixing ratio (kg/kg).
% rh = 20 hz relative humidity (%).
% mrh = mean relative humidity (%).
% mco2 = molecular weight of CO2 (gm).
% mixco2 = 20 hz CO2 mixing ratio (kg/kg). 'Webb' correction
% has been used in conversion from mas density on a
% point by point basis as per Fairall et al 2000, equ 62.
% ustar = friction velocity (m/s).
% qstar = H2O flux, normalised by ustar (kg/m^3).
% ============================================================
% Calculate RH - q relationship
% ============================================================
% Sort mixrat and rh before determining fit.
[P,S,MU] = polyfit(mixrat,rh,1); % Approximately linear
% relationship.
drh_by_dq = P(1)/MU(2);
% ============================================================
% Detrend CO2 mixing ratio with respect to RH
% ============================================================
% Sort rh and mixco2 before determining fit.
[P,S,MU] = polyfit(rh,mixco2,3); % calculate the polynomial
slope = P(1) / MU(2);
rhhat = (rh-MU(1))/MU(2);
response = (P(1).*rhhat.^3 + P(2).*rhhat.^2 + P(3).*rhhat);
mixco2 = mixco2 - response; % remove trend from co2 mix rat

% ============================================================
% Calculate detrended CO2 PSD and flux
% ============================================================
% calculate CO2 power spectrum density
% Use eddy corelation techniques to calculate the flux of the
% detrended CO2 (cflux: kg/(m^2.s)) and hence cstar (kg/m^3).
Cp=detrend(mixco2);
Qp=detrend(mixrat);
Cwq=nancov1(W,Qp);
Cwc=nancov1(W,Cp);
wq=Cwq(1,2);
wc=Cwc(1,2);   %m/s kg/kg
dc_by_dq = wc/wq;
    
% dc_by_dq = cstar/qstar;
dc_by_drh = dc_by_dq / drh_by_dq; % determine c*/RH* ratio
% (Equ 2, Prytherch et al, 2009)
cflux = wc *1e6*rhoair/ mco2 / 1000; % mols/(m2.s)
cflux = cflux*60*60*24*365; % mols/(m2 yr)

% =============================================================
% Iterate CO2 till convergence
% =============================================================
t=1:length(mixco2);
cfluxold = cflux+2;
loop = 0;
while (abs(cflux - cfluxold) > 1 && loop < 10 &&  abs(cflux) < 1000 ); % begin loop
% while (abs(cflux - cfluxold) > .1 && loop < 200 &&  abs(cflux) < 1000 ); % begin loop
% for ll=1:200
loop = loop + 1;
cfluxold = cflux;
mixco2nu = mixco2 + ((rh - mrh) .* dc_by_drh)/2; % adjust
% mixco2nu = mixco2 - ((rh - mrh) .* dc_by_drh)/2; % adjust
% gradient of CO2 mixing ratio with respect to humidity
% (Equ 2, Prytherch et al, 2009).
% Use eddy correlation methods to recalculate cstar
% and cflux from the adjusted CO2 mixing ratio.
Cp=detrend(mixco2nu);
           
Cwc=nancov1(W,Cp);
wc=Cwc(1,2);  
dc_by_dq = wc/wq;
dc_by_drh = dc_by_dq / drh_by_dq;  % new c*/rh* ratio
            
cflux = wc *1e6*rhoair/ mco2 / 1000; % mols/(m2.s)
cflux = cflux*60*60*24*365; % mols/(m2 yr)

       if 0
                figure(4)
                plot(t,Cp+mrhoc,'r',t,mixco2nu,'b')
                disp([num2str(loop) '   '  num2str(cflux) '  ' num2str(abs(cflux - cfluxold))])
%                 pause
       end

end

% rhoc=mixco2nu*1e6*rhoair;  %mg/m3
rhoc=mixco2nu*1e6;  %ppm