function [rhoc] = PKTcode(rhov,rhfast,rhoc,W)
% rhov=rhoair.*Rv7500_UC(rng7200UC);rhfast=rhfast7500_UC(rng7200UC);rhoc=rhoair.*Rc7500_UC(rng7200UC);W=wj(rngw7200UC);
% 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:
% rhov = Fast response H2O mixing ratio (kg/kg).
% rhfast = Fast response relative humidity (RH/100).
% mrh = mean relative humidity (%).
% mco2 = molecular weight of CO2 (gm).
% rhoc = Fast response CO2 mixing ratio (kg/kg).

plotit=0;

NaNs=find(isnan(rhoc));
if isempty(NaNs)
    rhfast=rhfast*100;
    mrh=mean(rhfast);
    mrhoc=mean(rhoc);
    mco2=44;
    Cp=detrend(rhoc);
    Cwc=cov(W,Cp);
    wc=Cwc(1,2);   %mg/m^2/s
    cflux = wc / mco2 / 1000;            % mols/(m2.s)
    cflux = cflux*60*60*24*365;         % mols/(m2 yr)
    if plotit
        figure(1);clf
        plot(rhfast,rhoc,'.')
        disp(['Prior to PKT: ' num2str(cflux)])
    end

    % ============================================================
    % Calculate RH - q relationship
    % ============================================================
    % Sort rhov and rh before determining fit.
    [P,S,MU] = polyfit(rhov,rhfast,1); % Approximately linear
    % relationship.
    drh_by_dq = P(1)/MU(2);
    % ============================================================
    % Detrend CO2 mixing ratio with respect to RH
    % ============================================================
    % Sort rh and rhoc before determining fit.
    [P,S,MU] = polyfit(rhfast,rhoc,1); %  calculate the polynomial
    slope = P(1) / MU(2);
    rhhat = (rhfast-MU(1))/MU(2);
    response = (P(1).*rhhat);
    t=1:length(rhoc);
    rhoc = rhoc - response; % remove RH trend from co2 mix rat
    if plotit
        figure(2);clf
        plot(rhfast,rhoc,'.')

        figure(4);clf
        plot(t,rhoc,'y.')
        hold
    end

    % ============================================================
    % Calculate detrended CO2 flux
    % ============================================================
    % 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(rhoc);
    Qp=detrend(rhov);
    Cwq=cov(W,Qp);
    Cwc=cov(W,Cp);
    wq=Cwq(1,2);
    wc=Cwc(1,2);   %mg/m^2/s
    dc_by_dq = wc/wq;
    dc_by_drh = dc_by_dq / drh_by_dq;   % determine c*/RH* ratio
    cflux = wc / mco2 / 1000;            % mols/(m2.s)
    cflux = cflux*60*60*24*365;         % mols/(m2 yr)
    % (Equ 2, Prytherch et al, 2010)

    if abs(cflux)<1000

        % =============================================================
        % Iterate CO2 till convergence
        % =============================================================

        cfluxold = cflux+2;
        loop = 0;

        if plotit
            disp([num2str(loop) '   '  num2str(cflux) '  ' num2str(abs(cflux - cfluxold))]);
        end

        while (abs(cflux - cfluxold) > .1 && loop < 10 && ...
                abs(cflux) < 1000 ); % begin loop
            loop = loop + 1;
            cfluxold = cflux;
            mixco2nu = rhoc + ((rhfast - mrh) .* dc_by_drh)/2; % adjust
            %     gradient of CO2 mixing ratio with respect to humidity
            % (Equ 2, Prytherch et al, 2010).

            % Use eddy correlation methods to recalculate cstar
            %     and cflux from the adjusted CO2 mixing ratio.
            Cp=detrend(mixco2nu);
            if plotit
                figure(4)
                plot(t,Cp+mrhoc,'r',t,mixco2nu,'b')
                disp([num2str(loop) '   '  num2str(cflux) '  ' num2str(abs(cflux - cfluxold))])
                pause
            end
            Cwc=cov(W,Cp);
            wc=Cwc(1,2);   %mg/m^2/s
            dc_by_dq = wc/wq;
            dc_by_drh = dc_by_dq / drh_by_dq;  % new c*/rh* ratio

            cflux = wc / mco2 / 1000; % mols/(m2.s)
            cflux = cflux*60*60*24*365; % mols/(m2 yr)
        end
        if plotit
            hold
            figure(3);clf
            plot(rhfast,mixco2nu,'.')
        end
    else
        mixco2nu=rhoc*NaN;
    end

    rhoc=mixco2nu;  %mg/m3
end


