% Some work to do here: Need to make sure the gridding is working properly.
% Looks a little bit off the raw data right now
clear all; close all;
gdata = readtable("Diel_GPP_log_QC_Revised.csv");
basepath = '\\atlas.shore.mbari.org\ProjectLibrary\901805_Coastal_Biogeochemical_Sensing\Spray_Data\';
mission = 3;
idiel = gdata.StartProf(mission):gdata.EndProf(mission); % before DIC starts decreasing
% idiel = gdata.Start(mission):gdata.End(mission); % before DIC starts decreasing
filename = char(gdata.Filename(mission));
filepath = fullfile(basepath,num2str(gdata.MissionID(mission)),'QC',filename);
load(filepath, 's'); % .mat file
%% recalc DIC with pHCorr
% need everything in vectors
ta = s.ta_canb(:);
pHin = s.pHin(:);
tc = s.tc(:);
psal = s.psal(:);
pres = s.pres(:);
sil = s.sil_canb(:);
po4 = s.po4_canb(:); 
pHcorr = s.pHcorr(:);
% a bit of QC
pHcorr(pHcorr > 8.8 | pHcorr < 7) = NaN;
% 17 Oct 2021 - updated to CO2SYSv3 (J. Sharp) - should have Yui double
% check the input options are correct -JSL
% CO2SYS(PAR1,PAR2,PAR1TYPE,PAR2TYPE,SAL,TEMPIN,TEMPOUT,PRESIN,PRESOUT,SI,PO4,NH4,H2S,pHSCALEIN,K1K2CONSTANTS,KSO4CONSTANT,KFCONSTANT,BORON)
trex = CO2SYSv3(ta, pHcorr, 1, 3, psal, tc, 20, pres, 0, sil, po4, 0, 0, 1, 10, 1, 2, 2);
% extract variables
s.dic_canb = reshape(trex(:,2),size(s.pres)); % CanyonB DIC
s.pH20atm = reshape(trex(:,20), size(s.pres)); % changed from 18 to 20 following CO2SYSv3 update
%% Add new variables
s.dt_UTC = datetime(s.sdn,'ConvertFrom','datenum','TimeZone','UTC');
s.dt_PT = s.dt_UTC;
s.dt_PT.TimeZone = "America/Los_Angeles";
s.sdn = datenum(s.dt_PT); % convert from UTC to local for these calculations. I could also just subtract 7/24 from original datenum. 
s.sigma0 = gsw_sigma0(s.abs_sal,s.cons_tc); % potential density anomaly (-1000 kg/m^3)
s.o2satconc = reshape(calcO2sat(s.tc(:), s.psal(:)), size(s.pres));
s.o2anom = s.doxy - s.o2satconc;
%% index time for diel cycles
time = s.sdn(idiel);
datetime = s.dt_PT(idiel);
depth = abs(s.depth(:,idiel)); % make depth positive
doxy = s.doxy(:,idiel);
o2satconc = s.o2satconc(:,idiel);
o2anom = s.o2anom(:,idiel);
sigma0 = s.sigma0(:,idiel);
tc = s.tc(:,idiel);
psal = s.psal(:,idiel);
dic = s.dic_canb(:,idiel);
dic(dic==-999) = NaN;
%% build data set
data = struct(); % main struct
%% Create scatter table
data.raw = struct();

data.raw.time = time;
data.raw.datetime = datetime;
data.raw.depth = depth;
data.raw.doxy = doxy;
data.raw.o2satconc = o2satconc;
data.raw.o2anom = o2anom;
data.raw.sigma0 = sigma0;
data.raw.tc = tc;
data.raw.psal = psal;
data.raw.dic = dic;

depth_scatter = data.raw.depth(:);
time_scatter = repmat(time,size(data.raw.depth,1),1);
time_scatter = time_scatter(:);
doxy_scatter = data.raw.doxy(:);
o2satconc_scatter = data.raw.o2satconc(:);
o2anom_scatter = data.raw.o2anom(:);
dic_scatter = data.raw.dic(:);
sigma0_scatter = data.raw.sigma0(:);
vars = {'time', 'depth', 'sigma0', 'doxy', 'o2satconc','o2anom', 'dic'};
data.raw.scatterTable = table(time_scatter, depth_scatter, sigma0_scatter, doxy_scatter,...
    o2satconc_scatter, o2anom_scatter, dic_scatter, 'VariableNames',vars);
idxnan = any(ismissing(data.raw.scatterTable), 2);
data.raw.scatterTable(idxnan,:) = [];

[data.raw.idx_zref, data.raw.idx_MLD, data.raw.ML_mask] = GetMLD(data.raw.depth, data.raw.sigma0, 'suga');

data.raw.meta = struct();
data.raw.meta.info = 'data for diel cycle time window ungridded';
data.raw.meta.scatterTableInfo = 'data as 1xn array';
%% Grid data
% Create regular grids for contouring
time_min = min(data.raw.scatterTable.time);
time_max = max(data.raw.scatterTable.time);
depth_min = 1; %round(min(data.Raw.scatterTable.depth), 0);
depth_max = round(max(data.raw.scatterTable.depth), 0) - 10;
time_grid_size = round((time_max-time_min) * 24, 0); % 1 hour resolution
depth_grid_size = length(depth_min:depth_max); % 1 m resolution

% Create regular grid
data.gridded = struct();
data.gridded.time = linspace(time_min, time_max, time_grid_size);
data.gridded.depth = linspace(depth_min, depth_max, depth_grid_size);

[data.gridded.time_grid, data.gridded.depth_grid] = meshgrid(data.gridded.time, data.gridded.depth);

% List of variables to interpolate
interp_vars = {'sigma0', 'doxy', 'o2satconc', 'o2anom', 'dic'};

% Loop over each variable and interpolate using scatteredInterpolant
for i = 1:numel(interp_vars)
    varname = interp_vars{i};
    x = data.raw.scatterTable.time;
    y = data.raw.scatterTable.depth;
    v = data.raw.scatterTable.(varname);

    % Create interpolant
    F = scatteredInterpolant(x, y, v, 'linear', 'none');

    % Evaluate on grid
    data.gridded.(varname) = F(data.gridded.time_grid, data.gridded.depth_grid);
end

[data.gridded.idx_zref, data.gridded.idx_MLD, data.gridded.ML_mask] = GetMLD(data.gridded.depth_grid, data.gridded.sigma0, 'suga');


% Add metadata
data.gridded.meta = struct();
data.gridded.meta.gridinfo = 'linearly interpolated grid: 1 m depth x 1 hour time';
%% Calculate timeseries
data.raw.timeseries = struct();
masked_o2anom = data.raw.o2anom;
masked_o2anom(data.raw.ML_mask) = NaN;  % Keep only values < ref + threshhold
data.raw.timeseries.o2anom = mean(masked_o2anom,1,'omitmissing');
data.raw.timeseries.o2anom_anom = data.raw.timeseries.o2anom - mean(data.raw.timeseries.o2anom);
masked_dic = data.raw.dic;
masked_dic(data.raw.ML_mask) = NaN;
data.raw.timeseries.dic = mean(masked_dic,1,'omitmissing');
data.raw.timeseries.dic_anom = mean(masked_dic,1,'omitmissing') - mean(data.raw.timeseries.dic);

data.gridded.timeseries = struct();
masked_o2anom = data.gridded.o2anom;
masked_o2anom(data.gridded.ML_mask) = NaN;  % Keep only values < ref + threshhold
data.gridded.timeseries.o2anom = mean(masked_o2anom,1,'omitmissing');
data.gridded.timeseries.o2anom_anom = data.gridded.timeseries.o2anom - mean(data.raw.timeseries.o2anom);

masked_dic = data.gridded.dic;
masked_dic(data.gridded.ML_mask) = NaN;
data.gridded.timeseries.dic = mean(masked_dic,1,'omitmissing');
data.gridded.timeseries.dic_anom = mean(masked_dic,1,'omitmissing') - mean(data.raw.timeseries.dic);
%% Save
datasetname = ['DielData_mission',num2str(mission),'.mat'];
fsavename = fullfile('C:\Users\bwerb\Documents\pH Diel Cycle Project\EDA2\data\',datasetname);
save(fsavename,"data",'-mat');