function d = calc_phcal_k2(data, hdr, info, tSDN, tPdir, fP_fit_order, k2fP_fit_order)

% Helper function for pH calibration GUI which generates k2 and f(P)
% coefficients, table result info and other needed itemsd info

% CHANGES:
%    11/23/22 JP - some minor updates to deal with incorprration of SBE raw
%       PT chamber cal data. f(P) data are now binned by 5 m bins to reduce
%       weighting issues due to unequal sample density along pressure axis

% ************************************************************************
% ************************************************************************
% TESTING
% %dirs.data = 'C:\Users\jplant\Documents\MATLAB\Apps\pH_Calibration\test data\';
% % % fn        = 'DF#306&VL#306_chB-201104.csv';
% % fn         = 'DF290&VL290chC_200612.csv';
% % fn        = 'GDF6_initial_cal.csv';
% dirs.data = '\\atlas\Chem\DuraFET\PTchamber\RawData\';
% fn        = 'DF#308&VL#308_chB-210209.csv';
% PTfile    = [dirs.data,fn];
% d = parse_PTchamber_file(PTfile);
% % 
% data = d.data;
% hdr  = d.hdr;
% info = d.info;
% fP_fit_order = 6;
% % 
% % 
% I.LSDN = find(strcmp(hdr,'Labview Date/Time') == 1);
% I.T    = find(strcmp(hdr,'Thermistor Temp') == 1);
% I.P    = find(strcmp(hdr,'Pressure (dBar)') == 1);
% I.VRS  = find(strcmp(hdr,'VRS VOLTS') == 1);
% I.VK  = find(strcmp(hdr,'VK VOLTS') == 1);
% % 
% % 
% % % SUBSET TO CAL ZONE
% % t1 = data(:,I.LSDN) >= datenum(2020,11,6,1,20,0) & ...
% %     data(:,I.LSDN) <= datenum(2020,11,10,6,49,0);
% % tSDN = data(:,I.LSDN) >= datenum(2020,6,14,1,31,0) & ...
% %     data(:,I.LSDN) <= datenum(2020,6,18,7,33,0);
% % tSDN = data(:,I.LSDN) >= datenum(2021,4,21,0,0,0) & ...
% %     data(:,I.LSDN) <= datenum(2021,4,21,11,0,0);
% 
% tSDN = data(:,I.LSDN) >= datenum(2021,2,11,3,18,0) & ... % DF308
%     data(:,I.LSDN) <= datenum(2021,2,15,9,46,0);
% 
% tPdir = true(size(data,1),1); % predim all
% % 
% %data = data(tSDN,:);


% ************************************************************************
% ************************************************************************
%                           DO SOME PREP WORK
% ************************************************************************
% ************************************************************************

% DEFINE DATA SUBSET WITH LOGICAL ARRAY BASED ON TIME & PRESSURE DIRECTION
tg     = tSDN & tPdir;

% DEFINE SOME DATA COLUMN INDICES
I.LSDN = find(strcmp(hdr,'Labview Date/Time') == 1);
I.T    = find(strcmp(hdr,'Thermistor Temp') == 1);
I.P    = find(strcmp(hdr,'Pressure (dBar)') == 1);
I.VRS  = find(strcmp(hdr,'VRS VOLTS') == 1);
I.VRSS = find(strcmp(hdr,'VRS STDDEV') == 1);
I.VK   = find(strcmp(hdr,'VK VOLTS') == 1);
I.VKS  = find(strcmp(hdr,'VK STDDEV') == 1);
I.IK   = find(strcmp(hdr,'IK AMPS') == 1);
I.IB   = find(strcmp(hdr,'IB AMPS') == 1);

% ************************************************************************
% EXTRACT SOME DATA ARRAYS FROM DATA MATRIX
Vrs   = data(:,I.VRS);
Vk    = data(:,I.VK);
IK    = data(:,I.IK);
IB    = data(:,I.IB);
Press = data(:,I.P);
Temp  = data(:,I.T);
HCl   = 0.01;

% TT = datetime(data(:,I.LSDN),'ConvertFrom','datenum');
% Pdir = [0;diff(data(:,I.P))];
% tUP  = Pdir < 0;
% tDOWN = Pdir > 0;
% tSTOP = Pdir == 0;

% DETERMINE IF THERE IS A DATA SUBSET FOR REPORTING INFO
Pdir = median(diff(Press(tg)));
if all(tPdir == 1)
    data_subset = 'ALL';
elseif Pdir > 0
    data_subset = 'DESCENDING';
elseif Pdir < 0
    data_subset = 'ASCENDING';
else
    data_subset = 'STATIONARY';
end

% GET IK & IB STATS
avg_IK = mean(IK(tg),'omitnan'); % amps
std_IK = std(IK(tg),'omitnan');
avg_IB = mean(IB(tg),'omitnan'); 
std_IB = std(IB(tg),'omitnan');

% DEFINE CONSTATNTS
%HCl  = 0.01; % Normal? This could be an input too I guess
% R    = 8.31446; % J/(mol K),                         From ph_calc
% F    = 96485;   % Faraday constant Coulomb / mol,    From ph_calc
R    = 8.3143; % J/(mol K),                         From excel wkb
F    = 96487;   % Faraday constant Coulomb / mol,    From excel wkb
Tk   = 273.15 + Temp; % degrees Kelvin
ln10 = log(10); % natural log of 10

% ************************************************************************
% ************************************************************************
%                  CALCULATE CALIBRATION VALUES
% ************************************************************************
% ************************************************************************

% MAY NOT NEED PRESS

% CALCULATE pH CALIBRATION K2 VALUES
% Vrs = Vc + k0 + k2T + f(P)
%
% Vc  = Vct + Vcp
% Vct = - (R*Tk/F)*ln(?[H]*?[Cl])
% Vcp = - ?Vo/23061*0.0242*P

% ************************************************************************
% CALC TEMPERATURE DEPENDENT PARAMETERS
dVo  = 16.38 + 0.0896*Temp - 0.001264*Temp.^2; % Partial molar vol HCL DI? REF?
aHCl = -0.000113*Temp + 0.9067; % Ref, where from?
Vct  = -(R.*Tk./F) .* log((HCl*aHCl).^2);       % why natural log?
Vcp  = -dVo./23061.*0.0242.*(Press-10)./10*1./1.01; % why Press - 10, 1.01, ref?

% ESTIMATE VRS & VK DRIFT & std dev Vk AT T~20 & -1> P < 10
tT = Temp  > 19.5 & Temp  < 20.5;
tP = Press > -1 & Press < 30;
tDrift = tT & tP & tSDN;
    
if sum(tDrift) > 2 % need at least 2 points for a line
    % Vrs drift
    % lsqfitY has some odd behvior with only 2 populations of X points and
    % and lots of spread compared to the regression slope
    % DF#317&VL#317_chB-210325.csv was the test case switching to poly fit
%     [my,by,ry,smy,sby] = lsqfity(data(tDrift,I.LSDN),data(tDrift,I.VRS));
%     Vrs_Drift = [my,by,ry,smy,sby];
%     dVRSdT      = data(:,I.LSDN)*my+by;
    
    
    Vrs_Drift = polyfit(data(tDrift,I.LSDN),data(tDrift,I.VRS),1); % slope, intercept
    dVRSdT      = data(:,I.LSDN)*Vrs_Drift(1) + Vrs_Drift(2);

    drift_anom = dVRSdT(1) - dVRSdT;
    Vrs_Dcor   = data(:,I.VRS) + drift_anom;
    %Vrs_Dcor   = data(:,I.VRS) + 0; % TESTING
    %Vrs_Vc = Vrs - Vct - Vcp;
    Vrs_Vc = Vrs_Dcor - Vct - Vcp;
    
%     % *****  TESTING   ******
%     figure
%     tstr = sprintf('Vrs drift (mV/day) = %0.4f & B = %0.0f\n N = %0.0f', ...
%         Vrs_Drift(1)*1000,Vrs_Drift(2),sum(tDrift));
%     
%     plot(data(tDrift,I.LSDN),data(tDrift,I.VRS),'bo')
%     title(tstr)
%     hold on
%     plot(data(tDrift,I.LSDN), dVRSdT(tDrift),'k--')
%     hold off
%     TD.data = data(tDrift,:);
%     TD.reg  = Vrs_Drift;
%     TD.I = I;
%     TD.hdr = hdr;
%     save('C:\temp\pHGUI_test_data.mat','TD');
%     
%     % *****  TESTING   ******
    
    
    % VK drift
%     [my,by,ry,smy,sby] = lsqfity(data(tDrift,I.LSDN),data(tDrift,I.VK));
%     Vk_Drift = [my,by,ry,smy,sby];
    Vk_Drift = polyfit(data(tDrift,I.LSDN),data(tDrift,I.VK),1); % slope, intercept
    
    % avg std Vk
    avg_std_Vk = mean(data(tDrift,I.VKS),'omitnan');
else
    disp('Not enough data to calculate drift in Vrs & Vk, setting M=NaN & B=NaN');
%     Vrs_Drift = [1,0,NaN,NaN,NaN];
%     Vk_Drift  = [1,0,NaN,NaN,NaN];
%     Vrs_Drift = [1,0];
%     Vk_Drift  = [1,0];
    
    Vrs_Drift = [NaN NaN];
    Vk_Drift  = [NaN NaN];    
    
    
    Vrs_Dcor  = data(:,I.VRS);
    Vrs_Vc = Vrs_Dcor - Vct - Vcp;
    avg_std_Vk = mean(data(:,I.VKS),'omitnan');
end


% [k2,k0,ry,sk2,sk0] = lsqfity(data(tSTOP&tP,I.T), Vrs_Vc(tSTOP&tP));
% 
% F1 = figure(1);
% %plot(TT(tT&tP), data(tT&tP,I.VRS),'bo')
% ax1(1) = subplot(2,1,1);
% scatter(TT(tT&tP&tSTOP), data(tT&tP&tSTOP,I.VRS), 9,data(tT&tP&tSTOP,I.VK),'filled')
% title(sprintf('Drift = %0.4f mV/day',my*1000));
% ylabel('Vrs (T~20C, P<10)')
% box on
% colorbar
% hold on
% plot(TT(tT&tP),data(tT&tP,I.LSDN)*my+by, 'k--')
% hold off
% 
% ax1(2) =subplot(2,1,2);
% scatter(TT(tT&tP&tSTOP), Vrs_Dcor(tT&tP&tSTOP), 9,data(tT&tP&tSTOP,I.VK),'filled')
% ylabel('Drift corrected Vrs')
% ylim(ax1(1).YLim);
% box on
% colorbar
% hold on
% plot(TT(tT&tP&tSTOP),data(tT&tP&tSTOP,I.LSDN)*my+by, 'k--')
% hold off
% 
% 
% F2 = figure(2);
% ax2(1) = subplot(1,2,1);
% scatter(data(tSTOP&tP,I.T), Vrs_Vc(tSTOP&tP), 9,data(tSTOP&tP,I.LSDN),'filled')
% ylabel('VrsCorr - Vc')
% xlabel('Temperature')
% %scatter(data(tSTOP,I.T), Vrs_Vc(tSTOP), 9,data(tSTOP,I.LSDN),'filled')
% %scatter(data(tSTOP,I.T), Vrs_Vc(tSTOP), 9,data(tSTOP,I.P),'filled')
% % scatter(data(tUP,I.T), Vrs_Vc(tUP), 9,data(tUP,I.P),'filled')
% %scatter(data(tDOWN,I.T), Vrs_Vc(tDOWN), 9,data(tDOWN,I.P),'filled')
% %scatter(data(:,I.T), Vrs_Vc, 9,data(:,I.P),'filled')
% box on
% colorbar



% ************************************************************************
% BIN THE DATA & GET K2 ESTIMATES
Pbins = [5,10,20,30,50,75,100,150,200:100:800,1000:200:2200]; % NO ZERO
%Pbins = [0,5,10,20,30,50,75,100,150,200:100:800,1000:200:2200];
%Tbins = 1:2:21;

out = ones(size(Pbins,2)-1,9)*NaN;
% F1 = figure(1); % for testing
% ax = subplot(1,1,1); % for testing
% xlabel('Temperature')
% ylabel('Vrs-Vc')
for i = 1:size(Pbins,2)-1
    t1 = data(:,I.P) >= Pbins(i) & data(:,I.P) < Pbins(i+1) & tg;
    if sum(t1) == 0 % no data for bin
        continue
    end
    TBIN = data(t1,I.T);
    minTBIN = min(TBIN);
    maxTBIN = max(TBIN);
    [k2,k0,ry,sk2,sk0] = lsqfity(TBIN, Vrs_Vc(t1));
    out(i,:) = [mean(data(t1,I.P)),k2,sk2,k0,sk0,ry.^2, sum(t1), ...
        minTBIN, maxTBIN];
    
%     % plot to check linearity of k2
%     plot(ax, TBIN, Vrs_Vc(t1),'bo') % TESTING  
%     title(sprintf('Pbin = %0.1f m  k2 = % 0.4e', out(i,1),k2));
%     hold(ax,'on')
%     plot(ax, xlim,xlim*k2+k0,'r-','LineWidth',3)
%     hold(ax,'off')
%     pause
end

tnan = isnan(out(:,1));
out  = out(~tnan,:);

% IF Tmin - Tmax < 2 for a given P reg, remove row, no spread in data
cols = size(out,2);
tEQ = out(:,cols) - out(:,cols-1) < 2;
if sum(tEQ) > 0
    out = out(~tEQ,:);
end
clear cols tEQ

if ~isempty(out) % NO DATA TO CALC k2 WITH (PROBABLY NO spread IN TEMP)
    %out = ones(1,9)*NaN;
    k0_HCL     = out(1,4); % surface most k2 intercept- what is best HCL k0?
    avg_k2     = mean(out(:,2),1,'omitnan');
    std_avg_k2 = std(out(:,2),1,'omitnan');
    
    % GET K2 PRESSURE RESPONSE FROM BIN AVGERAGED DATA & CALC pH PRESSURE
    % COEFFICIENTS FROM NON BINNED DATA (DOES THIS BIAS RESULT DUE TO
    % DIFFEENTIAL P SAMPLING DENSITY)
    if k2fP_fit_order == 0 % K2 is constant, no pressure response
        k2_Pcoefs = [0,avg_k2];
    else
        tnan = isnan(out(:,2));
        k2_Pcoefs = polyfit(out(~tnan,1), out(~tnan,2), k2fP_fit_order); % fit to all data
    end
else
    disp('NO k2 ESTIMATE! PROBABLY N0 SPREAD IN TEMPERATURE DATA')
    k0_HCL     = NaN;
    avg_k2     = NaN;
    std_avg_k2 = NaN;
    k2_Pcoefs  = NaN;
end



    
% tnan = isnan(out(:,2));
% [M,B,ry,sM,sB] = lsqfity(out(~tnan,1), out(~tnan,2)); % X = Pres, Y = K2
% dk2dP = [M,B,ry*ry,sM,sB];



% ***********************************************************************
% CALC PCOEF'S FROM NON BINNED DATA (DOES THIS BIAS DUE TO DIFFEENTIAL P
% SAMPLING DENSITY

%Vrs_Vc_k2T = Vrs_Vc - data(:,I.T).* avg_k2; % = f(p)
Vrs_Vc_k2T = Vrs_Vc - data(:,I.T).* polyval(k2_Pcoefs, data(:,I.P)); % = f(p)

% NO PRESSURE BINNING
%Pcoefs = polyfit(data(tg,I.P), Vrs_Vc_k2T(tg), fP_fit_order); % fit to all data

% TRY BINNING DATA TO SEE IF F(P) fit improves due to unequal dbar sampling
% weights
pbins    = (0:5:2500)';
binP     = (pbins(1:end-1) + pbins(2:end))/2;
newY     = ones(size(binP))*NaN;
std_newY = newY;
binN     = newY;

rawXY = [data(tg,I.P), Vrs_Vc_k2T(tg)];
for i = 1:size(binP,1)
    t1 = rawXY(:,1) >= pbins(i) & rawXY(:,1) < pbins(i+1);
    if sum(t1) == 0 % no data in bin
        continue
    end
    newY(i)     = median(rawXY(t1,2),'omitnan');
    binN(i)     = sum(t1);
    std_newY(i) = std(rawXY(t1,2),[],'omitnan');
end
tg_bin = ~isnan(newY);
Pcoefs = polyfit(binP(tg_bin), newY(tg_bin), fP_fit_order); % fit to all data
bin_fP = [binP(tg_bin), newY(tg_bin) , std_newY(tg_bin), binN(tg_bin)];


clear i t1 newY pbins binP rawXY binN std_newY





% PRESSURE BIN THE DATA BEFORE PCOEF FIT
% ************************************************************************
% BIN THE DATA & THEN GET PCOEFS
% PCbins = [0:5:2020];
% %PCbins = [0:5:200, 210:10:2000];
% %PCbins = [0:5:600, 610:10:1000, 1010:20:2000];
% %PCbins = [0:5:400, 410:10:1000, 1010:20:2000];
% 
% PCout    = ones(size(PCbins,2)-1,4)*NaN; % P,T, f(p), nbin
% for i = 1:size(PCbins,2)-1
%     t1 = data(:,I.P) >= PCbins(i) & data(:,I.P) < PCbins(i+1) & tg;
%     if sum(t1) == 0 % no data for bin
%         continue
%     end
%     PBIN = [data(t1,I.P), Vrs_Vc_k2T(t1)];
%     PCout(i,:) = [mean(PBIN,1,'omitnan'), std(PBIN(:,2),1,'omitnan'), sum(t1)];
%     PCout(i,:) = [median(PBIN,1,'omitnan'), std(PBIN(:,2),1,'omitnan'), sum(t1)];
% end
% tnan  = isnan(PCout(:,1));
% PCout = PCout(~tnan,:);
% Pcoefs = polyfit(PCout(:,1), PCout(:,2), fP_fit_order);

ssr    = sum((Vrs_Vc_k2T(tg) - polyval(Pcoefs, data(tg,I.P))) .^2,'omitnan');
nssr   = ssr ./ sum(tg);

Y     = Vrs_Vc_k2T(tg);
yfit  = polyval(Pcoefs, data(tg,I.P));         % Estimated  Regression Line
SStot = sum((Y-mean(Y,'omitnan')).^2);         % Total Sum-Of-Squares
SSres = sum((Y-yfit).^2);                      % Residual Sum-Of-Squares
RR    = 1-SSres/SStot;                         % R^2
clear Y yfit SStot SSres

cal_win     = [min(data(tg, I.LSDN)), max(data(tg, I.LSDN))];

% BUILD CELL ARRAY OF STRINGS FOR GUI OUTPUT
date_fmt = 'mm/dd/yyyy HH:MM';
results       = cell(40,3);
results(1,:)  = {'Process Date:',datestr(now, date_fmt), ''};
results(2,:)  = {'Durafet SN:',info.DF_SN, ''};
results(3,:)  = {'pH stem mfg:','MBARI', info.DF_TYPE};
if regexp(info.SourceDir,'K2_SBE_PTChamber_data','once')
    results(3,:)  = {'pH stem mfg:','SBE', info.DF_TYPE};
end

results(4,:)  = {'PT Source directory :',info.SourceDir, ''};
results(5,:)  = {'PT Source file :',info.SourceFile, ''};
results(6,:)  = {'PT test channel :',info.DF_CH, ''};
results(7,:)  = {'HCl concentration:', ...
    sprintf('%0.5f',HCl), 'Molar'};
results(8,:)  = {'k2 Calibration window:', datestr(cal_win(1), date_fmt), ...
    datestr(cal_win(2), date_fmt)};
results(9,:)  = {'Data subset:',data_subset, ''};
if ~isnan(Vrs_Drift(1))
    results(10,:)  = {'Vrs drift rate:', ...
        sprintf('%0.5f', Vrs_Drift(1)*1000), 'mV/day'};
    results(11,:)  = {'Vk drift rate:', ...
        sprintf('%0.5f', Vk_Drift(1)*1000), 'mV/day'};
else
    results(10,:)  = {'Vrs drift rate:','NO DATA', 'mV/day'};
    results(11,:)  = {'Vk drift rate:', 'NO DATA', 'mV/day'};
end

results(12,:)  = {'Avg std Vk', ...
    sprintf('%0.5f', avg_std_Vk*1e6), 'uVolts'};
results(13,:)  = {'Average k2:', ...
    sprintf('%0.4e', avg_k2), 'V/degC'};
results(14,:)  = {'Std dev k2:', ...
    sprintf('%0.4e', std_avg_k2), 'V/degC'};


% ADD k2 f(P) 
r = size(k2_Pcoefs,2);
tmp = flip(k2_Pcoefs);
for i = 1:r
    coef_str = sprintf('%0.0f',i-1);
    results(size(results,1)+i,:) ={['k2 f(P) c',coef_str], ...
        sprintf('%0.4e',tmp(i)),''};
end

% results(15,:)  = {'dk2dP M:', ...
%     sprintf('%0.4e', k2_Pcoefs(1)), 'V/degC/P'};
% results(16,:)  = {'dk2dP B:', ...
%     sprintf('%0.4e', k2_Pcoefs(2)), 'V/degC'};

results(size(results,1)+1,:)  = {'k0 HCL:', sprintf('%0.4f', k0_HCL),'Volts'};
results(size(results,1)+1,:)  = {'Avg Ik:', ...
    sprintf('%0.2f', avg_IK*1e9), 'nano amps'};
results(size(results,1)+1,:)  = {'Std dev Ik:', ...
    sprintf('%0.2f', std_IK*1e9), 'nano amps'};
results(size(results,1)+1,:)  = {'Avg Ib:', ...
    sprintf('%0.2f', avg_IB*1e9), 'nano amps'};
results(size(results,1)+1,:)  = {'Std dev Ib:', ...
    sprintf('%0.2f', std_IB*1e9), 'nano amps'};

results(size(results,1)+1,:)  = {'f(P) NSSR:',sprintf('%0.4e', nssr), ''};
results(size(results,1)+1,:)  = {'f(P) R*R:', sprintf('%0.4f', RR), ''};

r = size(Pcoefs,2);
tmp = flip(Pcoefs);
for i = 1:r
    coef_str = sprintf('%0.0f',i-1);
    results(size(results,1)+i,:) ={['f(P) c',coef_str], ...
        sprintf('%0.4e',tmp(i)),''};
end

t1 = ~cellfun(@isempty, results(:,1));
results = results(t1,:);

% ************************************************************************
% ************************************************************************
%  OUTPUT TO STRUCTURE
d.DF_SN       = info.DF_SN;
d.DF_CH       = info.DF_CH;
d.DF_TYPE     = info.DF_TYPE;
d.SourceDir   = info.SourceDir;
d.SourceFile  = info.SourceFile;
d.HCl         = HCl;
d.tDrift      = tDrift; % array to sub set data to determine ddrift
d.Vrs_Drift   = Vrs_Drift; % drift regression stats
d.Vrs_Dcor    = Vrs_Dcor;
d.Vrs_Vc      = Vrs_Vc;
d.Vk_Drift    = Vk_Drift; % drift regression stats
d.avg_std_Vk  = avg_std_Vk;
d.avg_Ik      = avg_IK;
d.std_Ik      = std_IK;
d.avg_Ib      = avg_IB;
d.std_Ib      = std_IB;

d.BinK2_hdr   = {'Avg Press' 'k2' 'std k2' 'k0' 'std k0' 'R*R' 'N' , ...
                 'Tmin' 'Tmax'};
d.BinK2_data  = out;
d.k2_Pcoefs   = k2_Pcoefs; %poly coef's
d.binfP_data  = bin_fP; %[bin P, bin Vrs_Vc_k2T, std Vrs_Vc_k2T, N]
d.binfP_hdr   = {'Avg Press' 'fP' 'std fP' ' bin N'};
d.Results     = results;
d.avg_k2      = avg_k2;
d.std_avg_k2  = std_avg_k2;
d.Vrs_Vc_k2T  = Vrs_Vc_k2T;
d.fP_SSR      = ssr;
d.fP_NSSR     = nssr;
d.fP_RR       = RR;
d.Pcoefs      = Pcoefs;
d.stem_mfg    = 'MBARI';
if regexp(info.SourceDir,'K2_SBE_PTChamber_data','once')
    d.stem_mfg    = 'SBE';
end
d.cal_win     = cal_win;
d.k0_HCL      = k0_HCL;

% OPTIMIZE K & f(P) models
d.k2_opt = get_opt_poly(d.BinK2_data(:,1), d.BinK2_data(:,2), 0:6, 5);
d.fP_opt = get_opt_poly(d.binfP_data(:,1), d.binfP_data(:,2), 0:8, 7);

%save('C:\temp\k2fP_test_data.mat', 'd') % TESTING

clearvars -except d

