function d = calc_phcal_k0(k0data, k0hdr, k0info, tSDN ,tON, k2, color_pH, time_offset, fit_order)
%
% THIS FUNCTION DETERMINES THE KO CALIBRATION CONSTANT FOR A DURAFET pH
% SENSOR.
%
% The intended use is as a helper function for  a Matlab app built in app
% designer
%
% INPUTS:
%   k0data   = structure of parsed k0 pON & POFF seawater calibration 
%              data from txt file
%   k0hdr    = associated hdr cell array identifying column data
%   tSDN     = logical array defining time window subset to use for k0
%              determination
%   tpON     = logical array (1 = pump ON)
%   color_pH = n x 2 array [matlab SDN, colorimetric calibration pH
%              values], ph on Total scale, 20C

% ************************************************************************
% ************************************************************************
% % % TESTING
% % % get calibration data
% % % GENERATE KO TEST DATA
% % % dirs.data = 'C:\Users\jplant\Documents\MATLAB\Apps\pH_Calibration\test data\';
% % dirs.data = 'C:\temp\';
% % % fn          = 'MSC0904_DURA.dat';
% % % MSC         = '0904';
% % fn          = 'MSC0951_DURA.dat';
% % MSC         = '0951';
% time_offset = 7; % DST, GMT-7
% 
% dirs.data = 'C:\temp\';
% MSC         = '0931';
% %fn          = 'MSC0931_DURA.dat';
% fn          = 'MSC0931_DURA (2).dat'; %pump ON
% 
% % 
% % 
% k0file      = [dirs.data,fn];
% d           = parse_k0_SW_file(k0file); % PUMP ON DATA
% k0hdr       = d.hdr;
% k0data.pOFF  = d.data;
% k0info.MSC = MSC;
% k0info.SourceDir = dirs.data;
% k0info.pOFFfile = fn; 
% 
% % pause
% 
% d           = parse_k0_SW_file(regexprep(k0file,'DURA','DURA_2')); %PUMP OFF
% k0data.pON = d.data;
% k0info.pONfile = regexprep(fn,'DURA','DURA_2');
% 
% %GET INDICES
% I.SDN  = find(strcmp(k0hdr,'DATE TIME') == 1);
% I.SCT  = find(strcmp(k0hdr,'SAMPLE COUNTER') == 1);
% I.T    = find(strcmp(k0hdr,'CTD TEMP') == 1);
% I.S    = find(strcmp(k0hdr,'CTD SALINITY') == 1);
% I.VRS  = find(strcmp(k0hdr,'VRS VOLTS') == 1);
% 
% %merge & sort data files
% data   = [k0data.pON; k0data.pOFF];
% 
% tON    = [ones(size(k0data.pON(:,1))); ones(size(k0data.pOFF(:,1)))*0];
% [~,ia] = sort(data(:, I.SDN));
% data   = data(ia,:);
% tON    = tON(ia);
% % tSDN   = data(:,1) >= datenum(2020,8,10,17,42,00) & ...
% %     data(:,1) <= datenum(2020,8,12,18,34,00);
% %tSDN   = true(size(data(:,1)));
% clear k0data d
% 
% tSDN   = data(:,1) >= datenum(2021,7,15,20,0,0) & ...
%     data(:,1) <= datenum(2021,7,19,23,0,0);
% 
% %get discrete pH samples
% fp            = '\\atlas\chem\DuraFET\K0 for DSD Stems\SpecpH_data\';
% color_pH      = merge_discrete_pH(fp, 'K0',0);
% 
% %Supply k2 slope
% %k2 = -1.1622E-03; 
% %k2 = -1.07778e-3; %MSC 0931 JP ALL
% k2 = -1.0776e-3; %MSC 0931 mDF#308+mVL#308calSO_recal_JP.xlsx
% fit_order = 1;
% k0data = data;
% % clearvars -except k0data k0hdr k0info tSDN tON k2 color_pH time_offset fit_order

% ************************************************************************
% ************************************************************************

% ************************************************************************
% REALLY CAN ONLY USE k0 DATA WITHIN THE BOUNDS OF THE RELAVENT DISCRETE
% SAMPLE DATA (ACOUNTING FOR ANY TIME OFFSET BETWEEN DISCRETE SAMPLES 
% (DST(7) or PST(8) & kO DATA STREAM (USUALLY GMT)

% % GET INDICES
I.SDN  = find(strcmp(k0hdr,'DATE TIME') == 1);
I.SCT  = find(strcmp(k0hdr,'SAMPLE COUNTER') == 1);
I.T    = find(strcmp(k0hdr,'CTD TEMP') == 1);
I.S    = find(strcmp(k0hdr,'CTD SALINITY') == 1);
I.VRS  = find(strcmp(k0hdr,'VRS VOLTS') == 1);
I.VK   = find(strcmp(k0hdr,'VK VOLTS') == 1);
I.VKS  = find(strcmp(k0hdr,'VK STDDEV') == 1);
I.IK   = find(strcmp(k0hdr,'IK AMPS') == 1);
I.IB   = find(strcmp(k0hdr,'IB AMPS') == 1);


% ADJUST TIME OFFSET BETWEEN DISCRETE SAMPLES & KO SEA WATER DATA
color_pH(:,1) = color_pH(:,1) + time_offset/24;

% SUBSAMPLE DISCRETE pH FIRST TO EXPERIMET TIME
% t1 = color_pH(:,1) >= min(k0data(:, I.SDN),[],'omitnan') & ...
%     color_pH(:,1) <= max(k0data(:, I.SDN),[],'omitnan') & color_pH(:,3) == 0;

t1 = color_pH(:,1) >= min(k0data(tSDN, I.SDN),[],'omitnan') & ...
    color_pH(:,1) <= max(k0data(tSDN, I.SDN),[],'omitnan') & color_pH(:,3) == 0;

% No color pH data available for k0 run
% has it been copied over to chem yet?
if sum(t1) == 0
    disp('WARNING: No color pH data found for k0 time window?')
    color_pH = [min(k0data(tSDN, I.SDN),[],'omitnan'),NaN; ...
        max(k0data(tSDN, I.SDN),[],'omitnan'),NaN];
else
    color_pH = color_pH(t1,:);
end

% NEXT FLAG EXPERIMENT TIME TO DISCRETE SAMPLE BOUNDS
% t1 = k0data(:,I.SDN) >= min(color_pH(:,1),[],'omitnan') & ...
%     k0data(:,I.SDN) <= max(color_pH(:,1),[],'omitnan');
% 
% tg = tSDN & t1; % subset logical for k0 calibration

tg = tSDN; % subset logical for k0 calibration


% DO SOME MAD FILTERING BASED ON Ik & Ib within "tg" range
outlier_threshold = 4;

Ik_tmp = k0data(:, I.IK);
Ik_tmp(~tg) = NaN; % set out of SDN range to Nan
tIk_spike = isoutlier(Ik_tmp,'median','ThresholdFactor',outlier_threshold);

Ib_tmp = k0data(:, I.IB);
Ib_tmp(~tg) = NaN;
tIb_spike = isoutlier(Ib_tmp,'median','ThresholdFactor',outlier_threshold);

tbad           = tIk_spike | tIb_spike;

Ik_bad_ct  = sum(tIk_spike);
Ik_bad_pct = sum(tIk_spike) ./ sum(tg) * 100;

Ib_bad_ct  = sum(tIb_spike);
Ib_bad_pct = sum(tIb_spike) ./ sum(tg) * 100;

clear Ib_tmp Ik_tmp

tg = tg & ~tbad;

% ************************************************************************
SDN   = k0data(tg, I.SDN);
Temp  = k0data(tg, I.T);
Salt  = k0data(tg, I.S);
Vrs   = k0data(tg, I.VRS);
Vk    = k0data(tg, I.VK);
Vkstd = k0data(tg, I.VKS);
Ik    = k0data(tg, I.IK);
Ib    = k0data(tg, I.IB);
tONtg = logical(tON(tg));

% DO SOME MAD FILTERING BASED ON Ik & Ib
% outlier_threshold = 4;
% TV       = [Ik Ib];
% med_IbIk = median(TV,1,'omitnan');
% TV_anom  = abs(TV - ones(size(Ik)) * med_IbIk); % median anomaly for Ik Ib
% MAD      = median(TV_anom, 1, 'omitnan') *1.4826;
% 
% if MAD(1) == 0
%     MAD(1) = mean(TVanom(:,1),1,'omitnan')*1.253314;
% end
% if MAD(2) == 0
%     MAD(2) = mean(TVanom(:,2),1,'omitnan')*1.253314;
% end
% 
% Z_Ik = TV_anom(:,1) ./ MAD(1);
% Z_Ib = TV_anom(:,2) ./ MAD(2);
% 
% tIk_spike =  Z_Ik > outlier_threshold;
% Ik_outlier_pct = sum(tIk_spike)./ size(tIk_spike,1)*100;
% 
% tIb_spike =  Z_Ib > outlier_threshold;
% Ib_outlier_pct = sum(tIb_spike)./ size(tIb_spike,1)*100;
% 
% tbad = tIk_spike | tIb_spike;


% **** HERE ON OUT USING DATA SUBSET DEFINED BY "tg" &"~tbad" ARRAYS ****

% ************************************************************************
% GENERATE INTERPOLATED DISCRETE pH VALUES AT IN SITU TEMPERATURE FOR KO
% SEAWATER CALIBRATION DATA

pp    = polyfit(color_pH(:,1),color_pH(:,2), fit_order); % get coefficients
pHT20 = polyval(pp, SDN);
pHT   = pHT20 - 0.015*(Temp - 20); % USE LIAR + CO2SYS?

% NEXT GET pH AT IN SITU TEMPERATURE

% % plot to check
% Yest = polyval(pp, color_pH(:,1));
% Yanom = color_pH(:,2) - Yest;
% F1 = figure(1);
% F1.Units = 'Normalized';
% F1.Position = [0.1323 0.3492 0.4828 0.5650];
% 
% pH_SDNT = datetime(color_pH(:,1),'ConvertFrom','datenum');
% k0_SDNT = datetime(data(:,I.SDN),'ConvertFrom','datenum');
% 
% ax(1) = subplot(2,1,1);
% plot(pH_SDNT,color_pH(:,2),'bo',pH_SDNT,Yest,'b-')
% hold on
% plot(k0_SDNT(tg), pHT,'r-')
% hold off
% 
% ax(2) = subplot(2,1,2);
% pH_SDNT = datetime(color_pH(:,1),'ConvertFrom','datenum');
% plot(pH_SDNT,Yanom,'ko')


% ************************************************************************
%  SET SOME CONSTANTS
% ************************************************************************
%Universal gas constant, (R) , http://physics.nist.gov/cgi-bin/cuu/Value?r
R    = 8.31446; % J/(mol K) 
F    = 96485; %Faraday constant Coulomb / mol
Tk   = 273.15 + Temp; % degrees Kelvin
ln10 = log(10); % natural log of 10
    
% ************************************************************************
% CALCULATE PHYSICAL AND THERMODYNAMIC DATA
% Dickson, A. G., Sabine, C. L., & Christian, J. R. (2007). Guide to best
% practices for ocean CO2 measurements.
% ************************************************************************

% OK - IONIC STRENGTH OF SEAWATER (mol / kg H2O)
% Varified units by comparing to Dickson et al. 2007: Chap 5, p10 Table 2
% Dickson et al. 2007: Chap 5, p13 Eq 34
IonS = 19.924 .* Salt ./ (1000 - 1.005 * Salt);

% OK - MEAN SEAWATER SULFATE CONCENTRATION (mol / kg solution)
% This wants to be mol/kg-seawater  as KHSO4 is on that scale
% Dickson et al. 2007: Chap 5, p10 Table 2
Stotal = (0.14 / 96.062) .* (Salt / 1.80655);

% !!DIFFERNT!! - MEAN SEAWATER CHLORIDE CONCENTRATION  (mol / kg H20)
% this wants to be mol/kg H2O as activity is on mol/kg H2O scale
% Dickson et al. 2007: Chap 5, p10 Table 2
Cltotal = 0.99889 / 35.453 .* Salt / 1.80655; %(mol / kg solution)
Cltotal = Cltotal ./(1 - 0.001005 .* Salt);  % (mol / kg H20)

% OK - BISULFIDE DISSCIATION CONSTANT AT T,S AND IONIC STRENGTH(mol/kg solution)
% Dickson et al. 2007: Chap 5, p12 Eq 33
Khso4 = exp(-4276.1 ./ Tk + 141.328 - 23.093 .* log(Tk) + ...
        (-13856 ./ Tk + 324.57 - 47.986 .* log(Tk)) .* IonS .^ 0.5 + ...
        (35474 ./ Tk - 771.54 + 114.723 .* log(Tk)) .* IonS - ...
        2698 ./ Tk .* IonS .^ 1.5 + 1776 ./ Tk .* IonS .^ 2 + ...
        log(1 - 0.001005 .* Salt));

% !!! NOT REQUIRED BECAUSE PRESSURE = 0 !!!    
% Millero 1983 Chemical Oceanography vol 8
% partial molar volume and compressibility of HSO4 in seawater. 
% deltaVHSO4 = -18.03 + 0.0466 .* Temp + 0.000316 .* Temp .^ 2;
% KappaHSO4 = (-4.53 + 0.09 .* Temp) / 1000;
% 
% %%%%%%%  Press changed from dbar to bar here by / 10
% lnKhso4fac = (-deltaVHSO4 + 0.5 .* KappaHSO4 .* (Press / 10)) .* ...
%              (Press / 10) ./ (R * 10 .* Tk);
% 
% %  bisulfate association constant at T, S, P
% Khso4TPS = Khso4 .* exp(lnKhso4fac);

% OK - GAMMA +/- HCl, activity coefficient of HCl at T/S, P=1
% ADH is the Debye Huckel constant, calcualted as a polynomial 
% fit to data in Khoo et al. 1977, doi:10.1021/ac50009a016
% See Martz et al. 2010, DOI 10.4319/lom.2010.8.172, p175
% Typo in paper 2nd term should be e-4 not e-6
%  
ADH = (3.4286e-6 .* Temp .^ 2 + 6.7524e-4 .* Temp + 0.49172143); 

% OK
log10gammaHCl = -ADH .* sqrt(IonS) ./ (1 + 1.394 .* sqrt(IonS)) + ...
                (0.08885 - 0.000111 .* Temp) .* IonS;

% !!! NOT REQUIRED BECAUSE PRESSURE = 0 !!! 
% Millero 1983 partial molar volume of HCl in seawater
%deltaVHCl = 17.85 + 0.1044 .* Temp - 0.001316 .* Temp .^ 2;

% effect of pressure on activity coefficient of HCl, divide by 2 because
% its a mean activity coefficient, divide by 10 for units in the cm3 to F
% conversion.

% !!! NOT REQUIRED BECAUSE PRESSURE = 0 !!! 
%log10gammaHCLtP = log10gammaHCl + deltaVHCl.*(Press./10)./(R.*Tk.*ln10)./2./10;


% GET pH_free, mol/kg-seawater scale
phfree = pHT + log10(1 + Stotal ./ Khso4); %

% GET pH_free, mol/kg-H2O scale
phfree = phfree + log10(1 - 0.001005 .* Salt);

k0T = phfree - (log(Cltotal) ./ ln10) - (2 .* log10gammaHCl);
k0T = Vrs - k0T .* (R .* Tk ./ F .* ln10);
k0  = k0T - k2 .* Temp;


k0_pON_med = median(k0(tONtg),'omitnan');
k0_pON_std = std(k0(tONtg),'omitnan');
Vk_pON_avg = mean(Vk(tONtg),'omitnan');
Vkstd_pON_avg = mean(Vkstd(tONtg),'omitnan');

k0_pOFF_med = median(k0(~tONtg),'omitnan');
k0_pOFF_std = std(k0(~tONtg),'omitnan');
Vk_pOFF_avg = mean(Vk(~tONtg),'omitnan');
Vkstd_pOFF_avg = mean(Vkstd(~tONtg),'omitnan');

% GET IK & IB STATS
med_IK_pON = median(Ik(tONtg),'omitnan'); % amps
std_IK_pON = std(Ik(tONtg),'omitnan');
med_IB_pON = median(Ib(tONtg),'omitnan'); 
std_IB_pON = std(Ib(tONtg),'omitnan');

med_IK_pOFF = median(Ik(~tONtg),'omitnan'); % amps
std_IK_pOFF = std(Ik(~tONtg),'omitnan');
med_IB_pOFF = median(Ib(~tONtg),'omitnan'); 
std_IB_pOFF = std(Ib(~tONtg),'omitnan');


% BUILD CELL ARRAY OF STRINGS FOR GUI OUTPUT
date_fmt = 'mm/dd/yyyy HH:MM';
results       = cell(30,3);

results(1,:)  = {'MSC # :',k0info.MSC, ''};
results(2,:)  = {'k0 calibration window', ...
                 datestr(min(k0data(tg,I.SDN)),'mm/dd/yyyy HH:MM'), ...
                 datestr(max(k0data(tg,I.SDN)),'mm/dd/yyyy HH:MM')};
results(3,:)  = {'Color pH time offset', num2str(time_offset,'%0.0f'), 'hours'};
results(4,:)  = {'k0 poly fit', num2str(fit_order,'%0.0f'), ''};
results(5,:)  = {'outlier_threshold', num2str(outlier_threshold,'%0.0f'), ''};
results(6,:)  = {'Ik outlier count', num2str(Ik_bad_ct,'%0.0f'), ''};
results(7,:)  = {'Ik outlier %', num2str(Ik_bad_pct,'%0.2f'), ''};
results(8,:)  = {'Ib outlier count', num2str(Ib_bad_ct,'%0.0f'), ''};
results(9,:)  = {'Ib outlier %', num2str(Ib_bad_pct,'%0.2f'), ''};
results(10,:)  = {'k0 Source directory :',k0info.SourceDir, ''};
results(11,:)  = {'k0 pump ON file', k0info.pONfile, ''};
results(12,:)  = {'k0 pump OFF file',k0info.pOFFfile, ''};
results(13,:) = {'','',''};
results(14,:)  = {'median k0 pump ON',num2str(k0_pON_med,'%0.4f') 'Volts'};
results(15,:)  = {'std k0 pump ON',num2str(k0_pON_std,'%0.4e') 'Volts'};
results(16,:) = {'Avg Vk pump ON', num2str(Vk_pON_avg,'%0.4f'), 'Volts'};
results(17,:) = {'Avg std Vk pump ON', num2str(Vkstd_pON_avg*1e6,'%0.1f'), 'uVolts'};
results(18,:)  = {'median Ik pump ON:',sprintf('%0.2f', med_IK_pON*1e9), 'nano amps'};
results(19,:)  = {'std Ik pump ON:',sprintf('%0.2f', std_IK_pON*1e9), 'nano amps'};
results(20,:)  = {'median Ib pump ON:',sprintf('%0.2f', med_IB_pON*1e9), 'nano amps'};
results(21,:)  = {'std Ib pump ON:',sprintf('%0.2f', std_IB_pON*1e9), 'nano amps'};
results(22,:) = {'','',''};
results(23,:)  = {'median k0 pump OFF',num2str(k0_pOFF_med,'%0.4f') 'Volts'};
results(24,:)  = {'std k0 pump OFF',num2str(k0_pOFF_std,'%0.4e') 'Volts'};
results(25,:)  = {'Avg Vk pump OFF', num2str(Vk_pOFF_avg,'%0.4f'), 'Volts'};
results(26,:)  = {'Avg std Vk pump OFF', num2str(Vkstd_pOFF_avg*1e6,'%0.1f'), 'uVolts'};
results(27,:)  = {'median Ik pump OFF:',sprintf('%0.2f', med_IK_pOFF*1e9), 'nano amps'};
results(28,:)  = {'std Ik pump OFF:',sprintf('%0.2f', std_IK_pOFF*1e9), 'nano amps'};
results(29,:)  = {'median Ib pump OFF:',sprintf('%0.2f', med_IB_pOFF*1e9), 'nano amps'};
results(30,:)  = {'std Ib pump OFF:',sprintf('%0.2f', std_IB_pOFF*1e9), 'nano amps'};



% X =k0_SDNT(tg);
% plot(X(tONtg), k0(tONtg),'r-',X(~tONtg), k0(~tONtg),'b-')

out_fill = ones(size(k0data(:,1)))*NaN;


% FILL OUTPUT STRUCTURE
d.avg_k0_pON     = [mean(k0(tONtg)), std(k0(tONtg))];
d.avg_k0_pOFF    = [mean(k0(~tONtg)), std(k0(~tONtg))];
d.k0_data        = k0data;
d.tON            = tON;
d.tk0            = tg;
d.k0_cal_win     = [min(k0data(tg,I.SDN)), max(k0data(tg,I.SDN))];
d.C_pH           = color_pH;
d.fit_order      = fit_order;
d.coef_C_pH      = pp;

d.int_C_pH20     = out_fill;
d.int_C_pH20(tg) = pHT20;
d.int_C_pH       = out_fill;
d.int_C_pH(tg)   = pHT;
d.k0             = out_fill;
d.k0(tg)         = k0;
d.Ik_bad_ct      = Ik_bad_ct;
d.Ik_bad_pct     = Ik_bad_pct;
d.Ib_bad_ct      = Ib_bad_ct;
d.Ib_bad_pct     = Ib_bad_pct;


d.MSC            = k0info.MSC;
d.SourceDir      = k0info.SourceDir;
d.pONfile        = k0info.pONfile;
d.pOFFfile       = k0info.pOFFfile;
d.Results        = results;
