clear all; close all;
gdata = readtable("Diel_GPP_log_QC_Revised.csv");
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));
basepath = '\\atlas.shore.mbari.org\ProjectLibrary\901805_Coastal_Biogeochemical_Sensing\Spray_Data\';
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
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;
%% Create variables for calculations
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;
%% Calculate Mixed layer depth
% MLD is super important and changes the PQ and RQ significantly
% Methods to calculate MLD
% suga et al 2004: zref = 10m, deltasigma0 > 0.125 kg/m^3
% Thomson and Fine 2003: zref = 2.5, deltasigma0 = 0.03 kg/m^3
[idx_zref, idx_MLD, ML_mask] = GetMLD(depth, sigma0, 'suga');
%% Create scatter plot to show O2 and DIC with Zref and Zsl
% Need to decide if it's better to show Doxy or Anom for this plot
% Showing both in the two types of plots below
depth_scatter = depth(:);
time_scatter = repmat(time,height(s.depth),1);
time_scatter = time_scatter(:);
doxy_scatter = doxy(:);
o2satconc_scatter = o2satconc(:);
o2anom_scatter = o2anom(:);
DIC_scatter = DIC(:);

vars = {'Time', 'Depth', 'Doxy', 'O2satconc','O2anom', 'DIC'};
T_scatter = table(time_scatter, depth_scatter, doxy_scatter,...
    o2satconc_scatter, o2anom_scatter, DIC_scatter, 'VariableNames',vars);
idxnan = isnan(T_scatter.Time);
T_scatter(idxnan,:) = [];

fig = figure;
tiledlayout(2,1)
nexttile
scatter(T_scatter,"Time","Depth","filled","ColorVariable","Doxy","Marker","o")
hold on
plot(time,depth(idx_zref),'k--',LineWidth=2)
plot(time,depth(idx_MLD),'k-',LineWidth=2)
legend('','Z_r_e_f','Z_S_L','Location','southwest')
set(gca,'Ydir','reverse')
ylim([0 60])
ylabel('Depth [m]')
xlabel('Time [Local]')
title('Dissolved Oxygen Concentration')
colormap("jet")
cb = colorbar;
cb.Label.String = '[\mumol/kg]';
clim_min = 240;
clim_max = 260;
% Set color limits
clim([clim_min clim_max])
cb.Ticks = clim_min:2:clim_max;
datetick

nexttile
scatter(T_scatter,"Time","Depth","filled","ColorVariable","DIC","Marker","o")
hold on
plot(time,depth(idx_zref),'k--',LineWidth=2)
plot(time,depth(idx_MLD),'k-',LineWidth=2)
legend('','Z_r_e_f','Z_S_L','Location','southwest')
set(gca,'Ydir','reverse')
ylim([0 60])
ylabel('Depth [m]')
xlabel('Time [Local]')
title('Dissolved Inorganic Carbon')
colormap("jet")
cb = colorbar;
cb.Label.String = '[\mumol/kg]';
clim_min = 2020;
clim_max = 2030;
% Set color limits
clim([clim_min clim_max])
cb.Ticks = clim_min:1:clim_max;
datetick

basepath = 'C:\Users\bwerb\Documents\pH Diel Cycle Project\EDA2';
savefilename = fullfile(basepath, 'figures', ['SurfaceLayer_scatter','_mission_',num2str(mission),'.png']);
exportgraphics(fig, savefilename, 'Resolution', 600); % Adjust resolution as needed


%% We correct for O2 response time so I am going to average this mask to create my diel time series
masked_o2satconc = o2satconc;
masked_doxy = doxy;
masked_o2satconc(ML_mask) = NaN;  % Keep only values < ref + threshhold
masked_doxy(ML_mask) = NaN;  % Keep only values < ref + threshhold

masked_o2anom = o2anom;
masked_o2anom(ML_mask) = NaN;  % Keep only values < ref + threshhold

%% Norms
diel_o2satconc_timeseries = mean(masked_o2satconc,1,"omitmissing");
diel_o2satconc_timeseries_anom = diel_o2satconc_timeseries - mean(diel_o2satconc_timeseries);
diel_o2satconc_timeseries_detrend = detrend(diel_o2satconc_timeseries);

diel_Doxy_timeseries = mean(masked_doxy,1,"omitmissing");
diel_Doxy_timeseries_anom = diel_Doxy_timeseries - mean(diel_Doxy_timeseries);

diel_o2anom_timeseries = mean(masked_o2anom,1,'omitmissing');
diel_o2anom_timeseries_anom = diel_o2anom_timeseries - mean(diel_o2anom_timeseries);
%% O2
dates = s.sdn(idiel);
% conc = diel_o2satconc_timeseries_anom;
lat = mean(s.lat(idiel));
lon = mean(s.lon(idiel));
timezone = -7; % PT  is -7 UTC
conc = diel_o2anom_timeseries_anom;

[rates_O2,fig1] = dielFit(dates,conc,lat,lon,timezone);
title('O2 Anom')
grid on
% ylim([-1 1])

% Save the figure
basepath = 'C:\Users\bwerb\Documents\pH Diel Cycle Project\EDA2';
savefilename = fullfile(basepath, 'figures', ['GOPCRO','_mission_',num2str(mission),'.png']);
exportgraphics(fig1, savefilename, 'Resolution', 600);
%% DIC
masked_DIC = DIC;
masked_DIC(ML_mask) = NaN;

DIC_timeseries = mean(masked_DIC,1,"omitmissing");
DIC_timeseries_anom = DIC_timeseries - mean(DIC_timeseries);

DIC_timeseries_detrend = detrend(DIC_timeseries);
DICtimeseries_runmean_anom = runmean(DIC_timeseries_detrend,3);

% if mission == 7
%     conc = DIC_timeseries_anom;
% elseif mission ==8
%     conc = detrend(DIC_timeseries_anom);
% end

conc = DIC_timeseries_anom;
[rates_DIC,fig2] = dielFit(dates,-conc,lat,lon,timezone);
title('DIC Anom')
grid on
% ylim([-1 1])
% Save the figure
basepath = 'C:\Users\bwerb\Documents\pH Diel Cycle Project\EDA2';
savefilename = fullfile(basepath, 'figures', ['GCPCRC','_mission_',num2str(mission),'.png']);
exportgraphics(fig2, savefilename, 'Resolution', 600);
%% Final Calculations table
GOP = rates_O2.GPP;
CRO = rates_O2.CR;
GCP = rates_DIC.GPP;
CRC = rates_DIC.CR;
PQ = GOP./GCP;
RQ = CRO./CRC;

vars = {'GOP', 'CRO', 'GCP', 'CRC', 'PQ', 'RQ'};
fits = {'Linear', 'Sin', 'PvE'};
PP = table(GOP, CRO, GCP, CRC, PQ, RQ, 'VariableNames', vars, 'RowNames', fits);

disp(PP)
%% Plot with pcolor
% Create regular grids for contouring
time_min = min(time);
time_max = max(time);
depth_min = 0;
depth_max = 60;

% Create regular grid
[T_grid, D_grid] = meshgrid(linspace(time_min, time_max, 100), linspace(depth_min, depth_max, 100));

% Prepare scattered data (same as before but we'll use it for interpolation)
depth_scatter = depth(:);
time_scatter = repmat(time, height(s.depth), 1);
time_scatter = time_scatter(:);
doxy_scatter = doxy(:);
o2satconc_scatter = o2satconc(:);
o2anom_scatter = o2anom(:);
DIC_scatter = DIC(:);

% Remove NaN values
valid_idx = ~isnan(time_scatter) & ~isnan(depth_scatter) & ~isnan(o2anom_scatter) & ~isnan(DIC_scatter);
time_valid = time_scatter(valid_idx);
depth_valid = depth_scatter(valid_idx);
o2anom_valid = o2anom_scatter(valid_idx);
DIC_valid = DIC_scatter(valid_idx);

% Interpolate scattered data onto regular grid
O2_grid = griddata(time_valid, depth_valid, o2anom_valid, T_grid, D_grid, 'linear');
DIC_grid = griddata(time_valid, depth_valid, DIC_valid, T_grid, D_grid, 'linear');

fig = figure;
tiledlayout(2,1)

% First subplot - O2 saturation concentration with pcolor
nexttile
pcolor(T_grid, D_grid, O2_grid)
shading interp  % This makes it smooth, use 'flat' for discrete blocks
hold on
plot(time, depth(idx_zref), 'k--', 'LineWidth', 2)
plot(time, depth(idx_MLD), 'k-', 'LineWidth', 2)
legend('', 'Z_r_e_f', 'Z_S_L', 'Location', 'southwest')
set(gca, 'Ydir', 'reverse')
ylim([0 50])
ylabel('Depth [m]')
xlabel('Time [Local]')
title('O2 Anom')
colormap("jet")
cb = colorbar;
cb.Label.String = '[\mumol/kg]';
% clim([240 260])
% cb.Ticks = 240:2:260;
% datetick('x')
mn = round(min(masked_o2anom, [], 'all'),0);
mx = round(max(masked_o2anom, [], 'all'),0);
clim([mn mx])
cb.Ticks = mn:1:mx;
datetick('x',15,'keeplimits')

% Second subplot - DIC with pcolor
nexttile
pcolor(T_grid, D_grid, DIC_grid)
shading interp
hold on
plot(time, depth(idx_zref), 'k--', 'LineWidth', 2)
plot(time, depth(idx_MLD), 'k-', 'LineWidth', 2)
legend('', 'Z_r_e_f', 'Z_S_L', 'Location', 'southwest')
set(gca, 'Ydir', 'reverse')
ylim([0 50])
ylabel('Depth [m]')
xlabel('Time [Local]')
title('Dissolved Inorganic Carbon')
colormap("jet")
cb = colorbar;
cb.Label.String = '[\mumol/kg]';
if mission ==3
    clim([2028 2100])
    cb.Ticks = 2028:5:2100;
else
    mn = round(min(masked_DIC, [], 'all'),0);
    mx = round(max(masked_DIC, [], 'all'),0);
    clim([mn mx])
    cb.Ticks = mn:1:mx;
end
datetick('x',15,'keeplimits')

% Save the figure
basepath = 'C:\Users\bwerb\Documents\pH Diel Cycle Project\EDA2';
savefilename = fullfile(basepath, 'figures', ['SurfaceLayer_color','_mission_',num2str(mission),'.png']);
exportgraphics(fig, savefilename, 'Resolution', 600);