%input: mfet data, spec pH data, CTD data
% 1) CTD(datetime,S,T)
%    Spec(datetime,pH_tot,temp_used)
%    MFET(datetime,Vrse)
% 2) spec pH_tot -> pH_insitu
% 3) sync datetimes
% 4) k0 for each pH_insitu value
% 5) mean k0
% 6) pH calc for all MFET vrse
%output: averaged k0 and std
clear all; close all;
%% load data
CTD = vertcat(readCTD('2025-06.csv'),readCTD('2025-07.csv'));
CTD = table2timetable(CTD);
%%
% Define a new time vector starting at the closest prior 0,5,10,... minute
startTime = dateshift(CTD.utc_time(1), 'start', 'minute', 'nearest'); % use 'nearest' or 'previous' depending on need
startTime = startTime - minutes(mod(minute(startTime), 5));       % round to nearest multiple of 5

% Retime to 5-minute intervals with linear interpolation
CTD = retime(CTD, startTime:minutes(5):CTD.utc_time(end), "nearest");
%%
Spec = readSpec('20250728_FTk0.txt');
GetMfetData('SB100618a.txt');
MFET = data;
MFET.("MM/DD/YYYY HH:MM:SS").TimeZone = 'UTC'; % annoying naming situation
%% k2
k2_sens = -0.0010626;
%% Calc insitu pH
trex = CO2SYS(2270, Spec.pH_tot,...
    1, 3, Spec.Salinity,...
    20, Spec.Temp_used, 0, 0, 0, 0, 1, 10, 1); % used or samp?
pH_insitu = trex(:,18); % Grab insitu pH
Spec.pH_insitu = pH_insitu;

%% Intersecting times with spec pH

[~, iMFET, iCTD] = intersect(MFET.("MM/DD/YYYY HH:MM:SS"),CTD.utc_time);
MFET = MFET(iMFET,:);
CTD = CTD(iCTD,:);

[~, iCTD, iSPEC] = intersect(CTD.utc_time,Spec.DateTime);
[~, iMFET, ~] = intersect(MFET.("MM/DD/YYYY HH:MM:SS"),Spec.DateTime);

CTD_intersect = CTD(iCTD,:);
MFET_intersect = MFET(iMFET,:);
SPEC_intersect = Spec(iSPEC,:);
%% Calc k0 (one value for each spec pH sample match)
%
k0 = nan(height(SPEC_intersect.DateTime),1);
for i = 1:height(SPEC_intersect.DateTime)
    k0(i) = calc_k0_ext(MFET_intersect.Vrse(i), CTD_intersect.temp(i),...
CTD_intersect.sal(i), SPEC_intersect.pH_insitu(i), k2_sens);
end
%% Stats

k0_std = std(k0);
k0_mean = mean(k0);

%% Calculate pH
[pHext_tot, pHext_free] = calc_dfet_pHext_k0(MFET.Vrse,...
    CTD.temp, CTD.sal, k0_mean, k2_sens);

%% plot k0 data!
figure
hold on;
plot(MFET.("MM/DD/YYYY HH:MM:SS"),pHext_tot)
plot(SPEC_intersect.DateTime,SPEC_intersect.pH_tot,'r+')
title([nanmean(k0) nanstd(k0)])
hold off;
%% Plot mfet vs MLML Intake
figure 
hold on;
plot(MFET.("MM/DD/YYYY HH:MM:SS"),pHext_tot)
% plot(CTD_intersect.utc_time,CTD_intersect.ph, 'bo') % there sensor sucks
plot(SPEC_intersect.DateTime,SPEC_intersect.pH_tot,'r+')
hold off;
legend('pHin [tot]','Spec pHin [tot]')
%% Diagnostics
figure
tl = tiledlayout(5,1);
nexttile
plot(MFET.("MM/DD/YYYY HH:MM:SS"),MFET.Vrse)
ylabel('Vrse')
hold on
title([nanmean(k0) nanstd(k0)])
yyaxis right
plot(MFET.("MM/DD/YYYY HH:MM:SS"),MFET.Vk)
ylabel('Vk')
nexttile
plot(MFET.("MM/DD/YYYY HH:MM:SS"),MFET.Vrse_std)
ylabel('Vrse std')
hold on
yyaxis right
plot(MFET.("MM/DD/YYYY HH:MM:SS"),MFET.Vk_std)
ylabel('Vk std')
nexttile
plot(MFET.("MM/DD/YYYY HH:MM:SS"),MFET.TC_cont)
ylabel('Board Temp')
hold on
yyaxis right
plot(MFET.("MM/DD/YYYY HH:MM:SS"),MFET.Humidity)
ylabel('Humidity')
nexttile
plot(MFET.("MM/DD/YYYY HH:MM:SS"),MFET.Ik)
ylabel('Ik')
hold on
yyaxis right
plot(MFET.("MM/DD/YYYY HH:MM:SS"),MFET.Ib)
ylabel('Ib')
nexttile
plot(MFET.("MM/DD/YYYY HH:MM:SS"),MFET.VbatMain)
ylabel('main batt voltage')
hold on
yyaxis right
plot(MFET.("MM/DD/YYYY HH:MM:SS"),MFET.VbiasPos)
ylabel('bias batt pos')