

clear
close all

% s = load_spray_sat('22602901');

% load('\\atlas\projectlibrary\901805_Coastal_Biogeochemical_Sensing\Spray_Data\22602901\sat_files\22602901_sat.mat', 's');

load('22602901_QC.mat', 's');


% =================================== USER INPUTS ============================================


% index for diel cycle expeirment #1
idiel = 127:323; % 127:323 are 150 m dives
expname = 'Diel67_90_2022_06';


% % index for diel cycle experiment #2
idiel = 360:587; 
expname = 'Diel67_70_2022_06';

% ==========================================================================================


% 100% sat concentration
s.o2satconc = reshape(calcO2sat(s.tc(:), s.psal(:)), size(s.pres));
% oxygen anomaly; [O2] - [O2]_100%
s.o2anom = s.doxy - s.o2satconc;


% extract just surface
dates = s.sdn_(2,idiel)-7/24;
conc = s.o2anom(:,idiel);
lat = mean(s.lat_(2,idiel));
lon = mean(s.lon_(2,idiel));
timezone = -7;

% prange = [0 10; 10 20; 20 30; 30 40];
prange = [0 5; 5 10; 10 15; 15 20; 20 25; 25 30; 30 35; 35 40];

%%

for i = 1:size(prange,1)
    
    % get indicies for pressure range
    ipres = s.pres > prange(i,1) & s.pres < prange(i,2);
    
    % extract useful variables (o2, chl, par)
    o2anom_temp = s.o2anom;
    par_temp = s.par;
    chl_temp = s.chl_npq;
    tc_temp = s.tc;
    psal_temp = s.psal;
    pdens_temp = s.pdens;
    
    % make values outside of pres range NaN for easier averaging
    o2anom_temp(~ipres) = NaN; 
    par_temp(~ipres) = NaN;
    chl_temp(~ipres) = NaN;
    tc_temp(~ipres) = NaN;
    psal_temp(~ipres) = NaN;
    pdens_temp(~ipres) = NaN;
    
    % Get average values in the prsesure range
    o2anom_avg(i,:) = nanmean(o2anom_temp(:,idiel),1);
    par_avg(i,:) = nanmean(par_temp(:,idiel),1);
    chl_avg(i,:) = nanmean(chl_temp(:,idiel),1); 
    tc_avg(i,:) = nanmean(tc_temp(:,idiel),1); 
    psal_avg(i,:) = nanmean(psal_temp(:,idiel),1); 
    pdens_avg(i,:) = nanmean(pdens_temp(:,idiel),1); 
    
    
    % get average GPP cycle over 1 week. 
    
    % compiled GPP and CR
    ftime = floor(dates(1)) + mod(dates,1);

    [rates,residuals,varmat,fitline] = dielFit(ftime,o2anom_avg(i,:),lat,lon,timezone);
    
    % save data
    GPPavg(i,:) = rates.GPP;
    CRavg(i,:) = rates.CR;
    
    % get day to day GPP cycle
    
    % need another for loop
    
    
end

f1 = figure;

% depth profile of GPP
nexttile;
    plot(GPPavg(:,2), nanmean(prange,2), 'ko-');
    % dprof(); % what is this
    xlabel('GOP (mmol m^{-3} d^{-1})');
    ylabel('Pressure (dBar)');
    title(expname, 'interpreter', 'none');

nexttile;
    X = par_avg;
    X(X<1) = NaN;
    plot(nanmean(X,2), GPPavg(:,2), 'ko');
    xlabel('mean PAR');
    ylabel('GOP (mmol m^{-3} d^{-1}');
  
% saveas(f1, [expname,'_GPP_profile_with_PAR.jpg']);


% figure of time series of O2, and PAR

softyel = [244/255 241/255 134/255];

% figure('unit', 'inch', 'position', [0.5 0.5 14 7]);
figure(2)
yyaxis left
    hold on; box on;
    plot(dates, o2anom_avg(1,:), 'c');
    set(gca, 'color', 'none', 'Ycolor', 'w', 'Xcolor', 'w');
    set(gcf, 'color','k');
    ylabel('\color{cyan}O_2 anomaly (O_2 - O_{2,sat}) [\mumol/kg]');
    datetick;
    
 yyaxis right
    hold on; box on;    
    set(gca, 'color', 'none', 'Ycolor', 'w', 'Xcolor', 'w', 'yticklabel', []);   
    ylim([3 6]);

 % export_fig(['O2_anomaly_',expname,'.png'], '-transparent');   
    
    plot(dates, par_avg(1,:), '-', 'color', softyel);
    set(gcf, 'color','k');
    ylabel('PAR', 'color', softyel);
    datetick;
    ylim([-500 2000]);
    
% export_fig(['O2_anomaly_and_PAR_',expname,'.png'], '-transparent');

%%

% ============ calculate NCP based on dO2 =====================================


dO2anom_avg = diff(o2anom_avg,[],2);

% a bit of crude QC
dO2anom_avg(dO2anom_avg < -0.5 | dO2anom_avg > 1) = NaN;

par_avg_ncp = par_avg(:,1:end-1) + diff(par_avg,[],2)./2;

hr_elapsed = diff(s.sdn_(2,idiel),[],2).*24; % convert deltaSDN to deltahour.  
% calculate NCP by assuming dO2 = NCP; NCP = dO2 (umol/kg) / dens (kg/L) / hour
ncp_avg = dO2anom_avg ./ ((1000 + pdens_avg(:,2:end))./1000) ./ (ones(size(dO2anom_avg,1),1)*hr_elapsed);

% bin data to PAR
   
PARbin = [0:100:400 600:200:1000 1250 1500];

for i = 1:size(dO2anom_avg,1)
    ncpbin(i,:) = bin_data_to_X(par_avg_ncp(i,:), ncp_avg(i,:), PARbin);
      
end

% surface PI raw data + fit
pink = [250/255 0 214/255]; 

figure
hold on; box on;
set(gca, 'color', 'none', 'Ycolor', 'w', 'Xcolor', 'w');
set(gcf, 'color','k');
plot(par_avg_ncp, ncp_avg(1,:), '.', 'color', [0.5,0.5,0.5]);
plot(PARbin, ncpbin(1,:) - ncpbin(1,1), 'wo-', 'color', pink, 'markerfacecolor', pink, 'MarkerEdgeColor', 'k', 'markersize', 15);
ylim([-0.1 0.3]);
ylabel('GPP [mmol O_2 m^{-3} hr^{-1}]');
xlabel('PAR [uE m^{-2} sec^{-1}]');

 
% export_fig(['PI_surface_',expname,'.png'], '-transparent');

figure
hold on; box on;
% set(gca, 'color', 'none', 'Ycolor', 'w', 'Xcolor', 'w');
% set(gcf, 'color','k');
plot(par_avg_ncp, ncp_avg(1,:), '.', 'color', [0.5,0.5,0.5]);
plot(PARbin, ncpbin(1,:) - ncpbin(1,1), 'ko-', 'color', pink, 'markerfacecolor', pink, 'MarkerEdgeColor', 'k', 'markersize', 15);
ylim([-0.1 0.3]);
ylabel('GPP [mmol O_2 m^{-3} hr^{-1}]');
xlabel('PAR [uE m^{-2} sec^{-1}]');

% save_ms_plot(['PI_surface_',expname,'_whiteback.png'], 100);


% PI curves from different depths

figure
hold on; box on;
plot(PARbin, ncpbin(1:5,:) - ncpbin(1:5,1), 'o-');
xlabel('Mean PAR');
ylabel('NCP (mmol m^{-3} hr^{-1})');
ylim([-0.1 0.3]);


% ================== Calculate GPP per day ======================================

% get GPP rates per day, just based on upper 10 m of data
ipres = s.pres < 10; 
o2anom_ = s.o2anom;
    o2anom_(~ipres) = NaN;
    o2anom_10m = nanmean(o2anom_,1);
    o2anom_10m = o2anom_10m(idiel);
par_ = s.par;
    par_(~ipres) = NaN;
    par_10m = nanmean(par_,1);
    par_10m = par_10m(idiel);
    % max par
    par_10m_max = nanmax(par_,[],1);
    par_10m_max = par_10m_max(idiel);
    
chl_surf_ = s.chl_npq(1,idiel);


    
% loop through each day

sdn_days = floor(dates); % sdn rounded to whole days. will use to parse out day to day
usdns = unique(sdn_days); % unique days, will use to drive for loop. 

% presize matrix
gpp = NaN(3,length(usdns));
cr = gpp; % community respiration
par_avg10m = NaN(1,length(usdns));
par_max = par_avg10m;
par_int = par_avg10m; 
chl_surf = par_avg10m;

for i = 1:length(usdns)
    
    iday = sdn_days == usdns(i);
    
    rates = dielFit(dates(iday),o2anom_10m(iday),lat,lon,timezone);
    
    gpp(:,i) = rates.GPP;
    cr(:,i) = rates.CR;
    
    % average PAR
    par_avg10m(i) = nanmean(par_10m(iday));
    par_max(i) = nanmax(par_10m_max(iday));
    par_int(i) = trapz(par_10m(iday));
    chl_surf(i) = nanmean(chl_surf_(iday));
    
    
end
close all 

% plot GPP and PAR

figure; 

nexttile;
    hold on; box on;
    plot(par_avg10m, gpp(1,:), 'ko', 'markerfacecolor', 'b');
    plot(par_avg10m, gpp(2,:), 'ko', 'markerfacecolor', 'r');
    plot(par_avg10m, gpp(3,:), 'ko', 'markerfacecolor', 'k');
    legend('Linear', 'Sinusoidal', 'PvsE');
    xlabel('Mean PAR surface 10 m');
    ylabel('GPP in upper 10 m');

nexttile;
    hold on; box on;
    plot(par_int, gpp(1,:), 'ko', 'markerfacecolor', 'b');
    plot(par_int, gpp(2,:), 'ko', 'markerfacecolor', 'r');
    plot(par_int, gpp(3,:), 'ko', 'markerfacecolor', 'k');
    legend('Linear', 'Sinusoidal', 'PvsE');
    xlabel('Daily integrated PAR 10 m');
    ylabel('GPP in upper 10 m');


figure(2);

nexttile;
hold on; box on;
    % plot(par_avg10m, gpp(1,:), 'ko', 'markerfacecolor', 'b');
    plot(par_avg10m(2:end-1), gpp(2,2:end-1), 'ko', 'markerfacecolor', 'r');
    xlabel('Mean PAR surface 10 m');
    ylabel('GPP in upper 10 m (Sinusoidal)');
    ylim([1.3 2.3]);
    xlim([100 400]);
    
nexttile;
    plot(par_int(2:end-1), gpp(2,2:end-1), 'ko', 'markerfacecolor', 'b');
%     plot(par_int, gpp(2,:), 'ko', 'markerfacecolor', 'b');

    % plot(par_avg10m, gpp(3,:), 'ko', 'markerfacecolor', 'k');
    xlabel('Daily integrated PAR surface 10 m');
    ylabel('GPP in upper 10 m (Sinusoidal)');

    
%% depth profiles

par_ = s.par(:,idiel);
pres_ = s.pres(:,idiel);
par_(par_ < 1) = NaN;

figure_2panel;
    nexttile;
    hold on; box on;
    set(gca, 'color', 'none', 'Ycolor', 'w', 'Xcolor', 'w');
    set(gcf, 'color','k');
    plot(par_, pres_, '-', 'color', ones(1,3).*0.3, 'linewidth', 0.1);
    plot(nanmean(par_,2), nanmean(pres_,2), '-', 'color', softyel, 'linewidth', 4);
    xlim([0 1000]);
    dprof();
    ylim([0 50]);
    xlabel('PAR');
    ylabel('Pressure (dBar)');

    export_fig(['PAR_profile_and_mean_',expname,'.png'], '-transparent');

% depth profile of GPP
nexttile;
    hold on; box on;
    set(gca, 'color', 'none', 'Ycolor', 'w', 'Xcolor', 'w');
    plot(GPPavg(:,2), nanmean(prange,2), 'o-', 'color', green);
    dprof();
    xlabel('GPP (mmol m^{-3} d^{-1})');
    ylabel('Pressure (dBar)');
    ylim([0 50]);
    xlim([0 4]);
 
export_fig(['PAR_and_GPP_profile_',expname,'.png'], '-transparent');   



   
% figure for comparing daily integrated GPP values
X = par_avg10m(2:end-1);
Y = gpp(2,2:end-1);
cgpp = polyfit(X,Y,1);
[m,b,r,sm,sb] = lsqfitma(X,Y);
cgpp2 = [m,b];

green = [12/255 250/255 0];

% for depth profile daily integrated GPP


figure
hold on; box on;
set(gca, 'color', 'none', 'Ycolor', 'w', 'Xcolor', 'w');
set(gcf, 'color','k');
plot([150 350], polyval(cgpp2,[150 350]), ':', 'color', green);
plot(X, Y, 'o', 'markerfacecolor', green, 'markeredgecolor', 'none');
ylim([1.5 2.3]);
xlim([100 400]);
xlabel('Daily mean PAR');
ylabel('GPP [mmol O_2 m^{-3} hr^{-1}]');

export_fig(['Daily_integrated_GPP_and_meanPAR_',expname,'.png'], '-transparent');

%add depth profile GPP
X2 = par_avg;
X2(X2<1) = NaN;
plot(nanmean(X2(1:8,:),2), GPPavg(1:8,2), 'o', 'color', green);
ylim([0.3 2.5]);
xlim([0 400]);


export_fig(['Daily_integrated_GPP_and_meanPAR_with_profile_data',expname,'.png'], '-transparent');


% nexttile;
%     X = par_avg;
%     X(X<1) = NaN;
%     plot(nanmean(X,2), GPPavg(:,2), 'ko');
%     xlabel('mean PAR');
%     ylabel('GOP (mmol m^{-3} d^{-1}');
% 

% rough DIC cycle

s.sdic_canb = s.dic_canb .* (35./s.psal);


dic_diel = parse_to_diel(dates, s.dic_canb(1,idiel),24);
sdic_diel = parse_to_diel(dates, nanmean(s.sdic_canb(1:5,idiel)),24);
o2anom_diel = parse_to_diel(dates, s.o2anom(1,idiel),24);
ta_diel = parse_to_diel(dates, s.ta_canb(1,idiel),24);

dic_anom = dic_diel - nanmedian(dic_diel,2);
sdic_anom = sdic_diel - nanmedian(sdic_diel,2);

YO2 = nanmedian(o2anom_diel,1); 


figure
hold on; box on;
plot(1:24, nanmedian(dic_anom,1), 'ro');
plot(1:24, nanmedian(sdic_anom,1), 'bo');
plot(1:24, YO2 - mean(YO2), 'ko');




