%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 = readCTD('2024-05.csv'); % if statement here: 
% if ctd is present on the board then use that data. else pull the latest data from mlml ctd and use that instead.
Spec = readSpec('20240509_T51_pH_k0_spec.txt');
GetMfetData('20240501.txt');
MFET = data;
MFET.("MM/DD/YYYY HH:MM:SS").TimeZone = 'UTC'; % annoying naming situation
%% k2
k2_T51 = -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, ~] = 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,:);
%% Calc k0 (one value for each spec pH sample)
%
k0 = nan(height(Spec.DateTime),1);
for i = 1:height(Spec.DateTime)
    k0(i) = calc_k0_ext(MFET_intersect.Vrse(i), CTD_intersect.temp(i),...
        CTD_intersect.sal(i), Spec.pH_insitu(i), k2_T51);
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_T51);

%% plot!
figure
plot(MFET.("MM/DD/YYYY HH:MM:SS"),pHext_tot)
hold on
plot(Spec.DateTime,Spec.pH_insitu,'.')
legend('MBARI MFET','Spec pH Insitu')
%% more plot
figure
tl = tiledlayout(6,1);

nexttile
plot(MFET.("MM/DD/YYYY HH:MM:SS"),MFET.Vrse)
ylabel('Vrse')
hold on
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
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')

nexttile
plot(MFET.("MM/DD/YYYY HH:MM:SS"),MFET.VbiasPos)
ylabel('bias batt pos')
%%
[pHext_tot_compare, ~] = calc_dfet_pHext_k0(MFET_intersect.Vrse,...
    CTD_intersect.temp, CTD_intersect.sal, k0_mean, k2_T51);
%%
plot(MFET.("MM/DD/YYYY HH:MM:SS"),pHext_tot)
hold on
plot(CTD.utc_time,CTD.ph)
legend('MBARI MFET','MLML')
%% 
% I thought it was easy to read, and gave what was requested. The only bits I have for a 'to-do' involve plotting the spec samples as circles over the mFET pH plot, that adds the complexity of the QF being back in, but helps 'cherry pick' the samples if one is wildly off and skewing average
% In terms of stuff that we might want for diagnostics, I was going to add to it and then let you do the cleanup/streamlining, but basically:
% 4 panel plot, x axis is always time
% Plot 1: Vrs left axis, Vk right axis
% Plot 2: Vrs_std Left, Vk_std Right
% Plot 3: Board temp left,  and Hum right
% Plot 4: Ik and Ib
% Battery voltage plot? Main batt left, bias batt right (Note: This will be helpful for our self-deployed loggers where they require battery, will be less useful for nano modular style)
% A plot with avg K0 and plotted cutoffs at 1 and 2 STD to see what values are out of line