%% Load mfet, test tank ctd, and spec pH data then match times and combine
% Define the folder path and file naming convention
folderPath = '\\atlas.shore.mbari.org\ProjectLibrary\901805_Coastal_Biogeochemical_Sensing\TestTank_pH\surpH004\';
filePattern = fullfile(folderPath, 'S4*.txt'); % Adjust the extension if different

% Get a list of files that match the pattern
files = dir(filePattern);

% Initialize an empty table to hold combined data
combinedData = table();

% Loop through each matching file
for i = 1:length(files)
    % Get the full file path
    filepath = fullfile(folderPath, files(i).name);
    
    % Apply the GetMfetData function
    GetMfetData(filepath);
    
    data = data(data.("MM/DD/YYYY HH:MM:SS") >...
    data.("MM/DD/YYYY HH:MM:SS")(1) + day(1),:); % remove warmup time 7 days
    data.("MM/DD/YYYY HH:MM:SS").TimeZone = 'UTC';
    
    % Vertically concatenate the data
    combinedData = vertcat(combinedData, data);
end

data = combinedData;
clear combinedData i;

data = sortrows(data, "MM/DD/YYYY HH:MM:SS");



ctd = readtable('\\atlas.shore.mbari.org\ChemWebData\termserv\testtank\microcat\sbe37.log'); % CTD
vars = {'temp','cond','pres','sal','density','dow','ddmmmyyyy','HHMMSS'};
ctd.Properties.VariableNames = vars;
ctd.datetime = datetime(ctd.ddmmmyyyy, 'InputFormat', 'dd MMM yyyy') + ctd.HHMMSS;
ctd.datetime.TimeZone = 'America/Los_Angeles'; % PST/PDT
ctd.datetime.TimeZone = 'UTC'; % UTC
%%
% Get all filenames in the folder that start with 'TT'
filePattern = fullfile(folderPath, 'TT_*.txt');
fileList = dir(filePattern);

% Initialize arrays to store the output values and datetime values
spec_pH_values = [];
spec_Temp_used_values = [];
dt_values = [];

% Loop through each file and apply the readSpec function
for k = 1:length(fileList)
    % Get the filename
    baseFileName = fileList(k).name;
    fullFileName = fullfile(folderPath, baseFileName);
    
    % Apply the readSpec function
    [spec_pH, spec_Temp_used] = readSpec(fullFileName);
    
    % Store the output values
    spec_pH_values = [spec_pH_values; spec_pH];
    spec_Temp_used_values = [spec_Temp_used_values; spec_Temp_used];
    

    dt_values = [dt_values; datetime(str2double(fileList(k).name(4:7)),...
        str2double(fileList(k).name(8:9)),...
        str2double(fileList(k).name(10:11)),...
        str2double(fileList(k).name(13:14)),...
        30,0,TimeZone='UTC')];
end

% Create a table with appropriate variable names
specpH_tbl = table(dt_values, spec_pH_values, spec_Temp_used_values, ...
                    'VariableNames', {'Datetime','Spec_pH', 'Spec_Temp_Used'});
% % % specpH_tbl.Datetime = specpH_tbl.Datetime - hours(1); % offsetting
% spec sample due to spike when sensor is recovered
%% match indices
ctd.datetime = dateshift(ctd.datetime,"start","minute");
data.("MM/DD/YYYY HH:MM:SS") = dateshift(data.("MM/DD/YYYY HH:MM:SS"),"start","minute");
% Find indices of matching values
[idxctd] = ismember(ctd.datetime, data.("MM/DD/YYYY HH:MM:SS"));
[idxmfet] = ismember(data.("MM/DD/YYYY HH:MM:SS"),ctd.datetime);
ctd = ctd(idxctd,:);
data = data(idxmfet,:);
if length(ctd.temp) > length(data.SampNum)
    ctd = ctd(1:length(data.SampNum),:);
end
%% merge
data = horzcat(data,ctd(:,["temp","cond","pres","sal","density"]));
%% join spec table and data table
% Create a full datetime vector covering the range of both tables
allTimes = union(data.("MM/DD/YYYY HH:MM:SS"), specpH_tbl.Datetime);

% Reindex the first table
T1_full = table(allTimes, 'VariableNames', {'MM/DD/YYYY HH:MM:SS'});
T1_full = outerjoin(T1_full, data, 'Keys', 'MM/DD/YYYY HH:MM:SS', 'MergeKeys', true);

% Reindex the second table
T2_full = table(allTimes, 'VariableNames', {'Datetime'});
T2_full = outerjoin(T2_full, specpH_tbl, 'Keys', 'Datetime', 'MergeKeys', true);

% Combine the tables horizontally
data = [T1_full, T2_full(:, 2:end)];

%% calc pH insitu
trex = CO2SYSv3(2270, data.Spec_pH, 1, 3, data.sal, data.Spec_Temp_Used, data.temp, 0, 0, 0, 0, 0, 0, 1, 10, 1, 1, 1);
data.spec_pH_insitu = trex(:,20); % pH insitu is 20
data.spec_pH_insitu(data.spec_pH_insitu == -999) = NaN;
%% Calc pH
spec_dt_idx = data.("MM/DD/YYYY HH:MM:SS") == specpH_tbl.Datetime(1);
trex = CO2SYSv3(2400, spec_pH, 1, 3, data.sal(spec_dt_idx), spec_Temp_used, data.temp(spec_dt_idx), 0, 0, 0, 0, 0, 0, 1, 10, 1, 1, 1); % 2400 is TA from TT spec sample
spec_pH_insitu = trex(:,20); % pH insitu is 20
[data.k0, ~] = calc_k0_ext(data.Vrse, data.temp, data.sal, data.spec_pH_insitu, 'UseDefaultk2');
Pcoefs = [];
ik0 = ~isnan(data.k0);
k0start = data.k0(ik0);
k0start = k0start(1);
% k0s at: 337, 669. still waiting to run spec on original TT sample
[phfree,data.phtot]= phcalc_jp(data.Vrse, data.pres, data.temp, data.sal, k0start, -0.001048, Pcoefs); % calc pH
%% Deal with data gaps and warm up times
dt = [0; diff(data.("MM/DD/YYYY HH:MM:SS"))];
idxdaygap = find(dt > days(1));
idxwarmuptime = idxdaygap:idxdaygap+48; % 2 days for warmup
data.phtot(idxwarmuptime) = NaN; % pH in warmup rows are NaN
%% Plot generation
fig = figure;
clf; cla;
tiledlayout(3,1)

% Determine the x-limits based on the datetime range
x_limits = [min(data.("MM/DD/YYYY HH:MM:SS")), max(data.("MM/DD/YYYY HH:MM:SS"))];

% First plot
nexttile
scatter(data.("MM/DD/YYYY HH:MM:SS"), data.phtot,'.')
hold on
scatter(data.("MM/DD/YYYY HH:MM:SS"), data.spec_pH_insitu, 15, 'filled')
title('surpH004 Test Tank Time Series')
legend('pHin', 'spec pHin','Location','southeast')
grid on
xlim(x_limits);  % Set xlim to the range of the datetime data

% Second plot
nexttile
data.pH_residuals = data.spec_pH_insitu - data.phtot;
scatter(data.("MM/DD/YYYY HH:MM:SS"), data.pH_residuals, 15, 'or', 'filled')
hold on
grid on
yline(0, 'k--', 'LineWidth', 1);
ylim([-.05 .05])
yticks(-0.05:0.01:0.05);
ylabel('spec pH - surpH');
title('pH Insitu Residuals vs Time');
xlim(x_limits);  % Set xlim to the range of the datetime data

% Third plot
nexttile
scatter(data.("MM/DD/YYYY HH:MM:SS"), data.k0, 15, 'or', 'filled')
grid on
xlabel('k0')
title('k0 vs Time')
xlim(x_limits);  % Set xlim to the range of the datetime data

% Save the figure
savepath = fullfile(folderPath, 'figures', ...
    ['surpH004_TT_', datestr(datetime('now'), 'YYYYmmDD'), '.png']);
% gdrivepath = ["G:\Shared drives\CarbonGroup\Werb\TestTankpHArray\surpH004"];
% gdrivesavepath = fullfile(gdrivepath, 'figures', ...
%     ['surpH004_TT_', datestr(datetime('now'), 'YYYYmmDD'), '.png']);


set(gcf, 'PaperOrientation', 'landscape');
% Set the paper size and position to match the figure size
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperSize', [11 8.5]); % Standard letter size in landscape
set(gcf, 'PaperPosition', [0 0 11 8.5]); % Fill the whole paper

saveas(fig,savepath); % Project library
% saveas(fig,gdrivesavepath); % Google drive
