%% %%%%%%%%%%%% TARA Microbiome %%%%%%%%%%%% %%
cd('/Users/gui/Documents/MATLAB/InLineAnalysis/InLineAnalysis-master/')

% Load InLineAnalysis and the configuration
ila = InLineAnalysis('cfg/TaraMicrobiome_cfg.m');

%%
ila.cfg.days2run = datenum(2022,6,6):datenum(2022,6,18);

%% load ACS57
ila.cfg.instruments2run = {'FLOW','ACS57'};
ila.cfg.qcref.view = 'ACS57';
ila.cfg.parallel = Inf;

%% 
ila.Read('prod');
acs_newlight = ila.instrument.ACS57.prod.p.ap;

%%
ila.DiagnosticPlot('AC',{'prod'}, false); % AC or BB



%% load ACS348
ila.cfg.instruments2run = {'FLOW','ACS348'};
ila.cfg.qcref.view = 'ACS348';
ila.cfg.parallel = Inf;

%% 
ila.Read('prod');

%% keep only common data
ila.instrument.ACS57.prod.p = ila.instrument.ACS57.prod.p(...
  ismember(ila.instrument.ACS57.prod.p.dt, ila.instrument.ACS348.prod.p.dt),:);
ila.instrument.ACS348.prod.p = ila.instrument.ACS348.prod.p(...
  ismember(ila.instrument.ACS348.prod.p.dt, ila.instrument.ACS57.prod.p.dt),:);

ila.instrument.ACS57.prod.p = sortrows(ila.instrument.ACS57.prod.p, 'dt');
ila.instrument.ACS348.prod.p = sortrows(ila.instrument.ACS348.prod.p, 'dt');

%%
ila.DiagnosticPlot('AC',{'prod'},false); % AC or BB



%% extract data for new light
ap_n_nan_newlight = sum(isnan(ila.instrument.ACS57.prod.p.ap))/size(ila.instrument.ACS57.prod.p.ap,1) < 0.3;
cp_n_nan_newlight = sum(isnan(ila.instrument.ACS57.prod.p.cp))/size(ila.instrument.ACS57.prod.p.cp,1) < 0.05;

lambda_a_newlight = ila.instrument.ACS57.lambda_a(ap_n_nan_newlight);
lambda_c_newlight = ila.instrument.ACS57.lambda_c(cp_n_nan_newlight);

ap_newlight = ila.instrument.ACS57.prod.p.ap(:, ap_n_nan_newlight);
cp_newlight = ila.instrument.ACS57.prod.p.cp(:, cp_n_nan_newlight);
ap_sd_newlight = ila.instrument.ACS57.prod.p.ap_sd(:, ap_n_nan_newlight);
cp_sd_newlight = ila.instrument.ACS57.prod.p.cp_sd(:, cp_n_nan_newlight);

ap_sigtonoi_newlight = mean(ap_newlight ./ ap_sd_newlight, 'omitnan');
ap_sigtonoi_sd_newlight = std(ap_newlight ./ ap_sd_newlight, 'omitnan') / sqrt(size(ap_newlight,1));
cp_sigtonoi_newlight = mean(cp_newlight ./ cp_sd_newlight, 'omitnan');
cp_sigtonoi_sd_newlight = std(cp_newlight ./ cp_sd_newlight, 'omitnan') / sqrt(size(cp_newlight,1));

% ap_sigtonoi_newlight = mean(ap_sd_newlight ./ ap_newlight, 'omitnan');
% ap_sigtonoi_sd_newlight = std(ap_sd_newlight ./ ap_newlight, 'omitnan') / sqrt(size(ap_newlight,1));
% cp_sigtonoi_newlight = mean(cp_sd_newlight ./ cp_newlight, 'omitnan');
% cp_sigtonoi_sd_newlight = std(cp_sd_newlight ./ cp_newlight, 'omitnan') / sqrt(size(cp_newlight,1));

ap_stderrmean_newlight = mean(ap_sd_newlight ./ ila.instrument.ACS57.prod.p.ap_n, 'omitnan');
ap_stderrmean_sd_newlight = std(ap_sd_newlight ./ sqrt(ila.instrument.ACS57.prod.p.ap_n), 'omitnan') / sqrt(size(ap_newlight,1));
cp_stderrmean_newlight = mean(cp_sd_newlight ./ ila.instrument.ACS57.prod.p.ap_n, 'omitnan');
cp_stderrmean_sd_newlight = std(cp_sd_newlight ./ sqrt(ila.instrument.ACS57.prod.p.cp_n), 'omitnan') / sqrt(size(cp_newlight,1));

%% extract data for old light
ap_n_nan_oldlight = sum(isnan(ila.instrument.ACS348.prod.p.ap))/size(ila.instrument.ACS348.prod.p.ap,1) < 0.3;
cp_n_nan_oldlight = sum(isnan(ila.instrument.ACS348.prod.p.cp))/size(ila.instrument.ACS348.prod.p.cp,1) < 0.05;
lambda_a_oldlight = ila.instrument.ACS348.lambda_a(ap_n_nan_oldlight);
lambda_c_oldlight = ila.instrument.ACS348.lambda_c(cp_n_nan_oldlight);

ap_oldlight = ila.instrument.ACS348.prod.p.ap(:, ap_n_nan_oldlight);
cp_oldlight = ila.instrument.ACS348.prod.p.cp(:, cp_n_nan_oldlight);
ap_sd_oldlight = ila.instrument.ACS348.prod.p.ap_sd(:, ap_n_nan_oldlight);
cp_sd_oldlight = ila.instrument.ACS348.prod.p.cp_sd(:, cp_n_nan_oldlight);

ap_sigtonoi_oldlight = mean(ap_oldlight ./ ap_sd_oldlight, 'omitnan');
ap_sigtonoi_sd_oldlight = std(ap_oldlight ./ ap_sd_oldlight, 'omitnan') / sqrt(size(ap_oldlight,1));
cp_sigtonoi_oldlight = mean(cp_oldlight ./ cp_sd_oldlight, 'omitnan');
cp_sigtonoi_sd_oldlight = std(cp_oldlight ./ cp_sd_oldlight, 'omitnan') / sqrt(size(cp_oldlight,1));

ap_stderrmean_oldlight = mean(ap_sd_oldlight ./ ila.instrument.ACS348.prod.p.ap_n, 'omitnan');
ap_stderrmean_sd_oldlight = std(ap_sd_oldlight ./ sqrt(ila.instrument.ACS348.prod.p.ap_n), 'omitnan') / sqrt(size(ap_oldlight,1));
cp_stderrmean_oldlight = mean(cp_sd_oldlight ./ ila.instrument.ACS348.prod.p.cp_n, 'omitnan');
cp_stderrmean_sd_oldlight = std(cp_sd_oldlight ./ sqrt(ila.instrument.ACS348.prod.p.cp_n), 'omitnan') / sqrt(size(cp_oldlight,1));

%% compare ap signal to noise ratios
% minute binned ap spectra averaged over 8 days
figure()
subplot(2,2,1)
shadedErrorBar(lambda_a_newlight, mean(ap_newlight, 'omitnan'), ...
  std(ap_newlight, 'omitnan') / ...
  sqrt(size(ap_newlight, 1)), 'lineProps', '-r', 'transparent', true)
shadedErrorBar(lambda_a_oldlight, mean(ap_oldlight, 'omitnan'), ...
  std(ap_oldlight, 'omitnan') / ...
  sqrt(size(ap_oldlight, 1)), 'lineProps', '-k', 'transparent', true)
ylabel('ap [m^{-1}]')
xlabel('lambda')
title('                  Minute binned ap \newline (~8 days average Benguela current) N = 13672')
legend('New light ACS57', 'Old light ACS348')
set(gca, 'FontSize', 14)

% minute binned ap_sd spectra averaged over 8 days
subplot(2,2,2)
shadedErrorBar(lambda_a_newlight, mean(ap_sd_newlight, 'omitnan'), ...
  std(ap_sd_newlight, 'omitnan') / sqrt(size(ap_sd_newlight, 1)), 'lineProps', '-r', 'transparent', true)
shadedErrorBar(lambda_a_oldlight, mean(ap_sd_oldlight, 'omitnan'), ...
  std(ap_sd_oldlight, 'omitnan') / sqrt(size(ap_oldlight, 1)), 'lineProps', '-k', 'transparent', true)
ylabel('ap std [m^{-1}]')
xlabel('lambda')
title('                  Minute binned ap\_sd \newline (~8 days average Benguela current) N = 13672')
legend('New light ACS57', 'Old light ACS348')
set(gca, 'FontSize', 14)

% minute binned signal to noise spectra averaged over 8 days
subplot(2,2,3)
shadedErrorBar(lambda_a_newlight, ap_sigtonoi_newlight, ap_sigtonoi_sd_newlight, ...
  'lineProps', '-r', 'transparent', true)
shadedErrorBar(lambda_a_oldlight, ap_sigtonoi_oldlight, ap_sigtonoi_sd_oldlight, ...
  'lineProps', '-k', 'transparent', true)
ylabel('ap/ap\_sd [unitless]')
xlabel('lambda')
title('   Signal to noise of Minute binned ap/ap\_sd \newline (~8 days average Benguela current) N = 13672')
legend('New light ACS57', 'Old light ACS348')
set(gca, 'FontSize', 14)

% minute binned ap_sd spectra averaged over 8 days
subplot(2,2,4)
hold on
plot(lambda_a_newlight, ap_stderrmean_newlight, '-r')
plot(lambda_a_oldlight, ap_stderrmean_oldlight, '-k')
% shadedErrorBar(lambda_a_newlight, ap_stderrmean_newlight, ...
%   ap_stderrmean_sd_newlight, 'lineProps', '-r', 'transparent', true)
% shadedErrorBar(lambda_a_oldlight, ap_stderrmean_oldlight, ...
%   ap_stderrmean_sd_oldlight, 'lineProps', '-k', 'transparent', true)
ylabel('ap stdev/sqrt(N) [m^{-1}]')
xlabel('lambda')
title('          Average standard error of the mean \newline (~8 days average Benguela current) N = 13672')
legend('New light ACS57', 'Old light ACS348')
set(gca, 'FontSize', 14)
set(gca, 'YScale', 'log')

%% compare cp signal to noise ratios
% minute binned cp spectra averaged over 8 days
figure()
subplot(2,2,1)
shadedErrorBar(lambda_c_newlight, mean(cp_newlight, 'omitnan'), ...
  std(cp_newlight, 'omitnan') / sqrt(size(cp_newlight, 1)), 'lineProps', '-r', 'transparent', true)
shadedErrorBar(lambda_c_oldlight, mean(cp_oldlight, 'omitnan'), ...
  std(cp_oldlight, 'omitnan') / sqrt(size(cp_oldlight, 1)), 'lineProps', '-k', 'transparent', true)
ylabel('cp [m^{-1}]')
xlabel('lambda')
title('                  Minute binned cp \newline (~8 days average Benguela current) N = 13672')
legend('New light ACS57', 'Old light ACS348')
set(gca, 'FontSize', 14)

% minute binned cp_sd spectra averaged over 8 days
subplot(2,2,2)
shadedErrorBar(lambda_c_newlight, mean(cp_sd_newlight, 'omitnan'), std(cp_sd_newlight, 'omitnan') / ...
  sqrt(size(cp_sd_newlight, 1)), 'lineProps', '-r', 'transparent', true)
shadedErrorBar(lambda_c_oldlight, mean(cp_sd_oldlight, 'omitnan'), std(cp_sd_oldlight, 'omitnan') / ...
  sqrt(size(cp_oldlight, 1)), 'lineProps', '-k', 'transparent', true)
ylabel('cp std [m^{-1}]')
xlabel('lambda')
title('                  Minute binned cp\_sd \newline (~8 days average Benguela current) N = 13672')
legend('New light ACS57', 'Old light ACS348')
set(gca, 'FontSize', 14)

% minute binned signal to noise spectra averaged over 8 days
subplot(2,2,3)
shadedErrorBar(lambda_c_newlight, cp_sigtonoi_newlight, cp_sigtonoi_sd_newlight, ...
  'lineProps', '-r', 'transparent', true)
shadedErrorBar(lambda_c_oldlight, cp_sigtonoi_oldlight, cp_sigtonoi_sd_oldlight, ...
  'lineProps', '-k', 'transparent', true)
ylabel('cp/cp\_sd [unitless]')
xlabel('lambda')
title('   Signal to noise of Minute binned cp/cp\_sd \newline (~8 days average Benguela current) N = 13672')
legend('New light ACS57', 'Old light ACS348')
set(gca, 'FontSize', 14)

% minute binned cp_sd spectra averaged over 8 days
subplot(2,2,4)
hold on
plot(lambda_c_newlight, cp_stderrmean_newlight, '-r')
plot(lambda_c_oldlight, cp_stderrmean_oldlight, '-k')
% shadedErrorBar(lambda_c_newlight, cp_stderrmean_newlight, cp_stderrmean_sd_newlight, ...
%   'lineProps', '-r', 'transparent', true)
% shadedErrorBar(lambda_c_oldlight, cp_stderrmean_oldlight, cp_stderrmean_sd_oldlight, ...
%   'lineProps', '-k', 'transparent', true)
ylabel('cp stdev/sqrt(N) [m^{-1}]')
xlabel('lambda')
title('          Average standard error of the mean \newline (~8 days average Benguela current) N = 13672')
legend('New light ACS57', 'Old light ACS348')
set(gca, 'FontSize', 14)
set(gca, 'YScale', 'log')

%% correlation ap & cp
% ap
wl = 400:20:710;
ap_newl = interp1(lambda_a_newlight, ap_newlight', wl, 'linear')';
ap_oldl = interp1(lambda_a_oldlight, ap_oldlight', wl, 'linear')';

col = reshape(spectrumRGB(wl), max(size(wl)),  3);
  
figure()
subplot(1,2,1)
hold on
for i = 1:size(ap_oldl, 2)
  sc = scatter(ap_newl(:,i), ap_oldl(:,i), 20, col(i,:), 'filled', ...
    'MarkerFaceAlpha', 0.5, 'MarkerEdgeAlpha', 0.5);
end
rline = refline(1, 0);
rline.Color = 'k';
rline.LineStyle = '--';
rline.LineWidth = 2;
xlim([min(min(ap_newl)) max(max(ap_newl))]);
ylim([min(min(ap_oldl)) max(max(ap_oldl))]);

legend(cellstr(strcat(num2str(wl'), ' nm')), 'location', 'southeast')
title('a_p in the Benguela current ~8 days average N = 13672')
xlabel('a_p new light ACS57 [m^{-1}]')
ylabel('a_p old light ACS348 [m^{-1}]')
set(gca, 'FontSize', 12)


% cp
wl = 400:20:710;
cp_newl = interp1(lambda_c_newlight, cp_newlight', wl, 'linear')';
cp_oldl = interp1(lambda_c_oldlight, cp_oldlight', wl, 'linear')';

col = reshape(spectrumRGB(wl), max(size(wl)),  3);
  
subplot(1,2,2)
hold on
for i = 1:size(cp_oldl, 2)
  scatter(cp_newl(:,i), cp_oldl(:,i), 20, col(i,:), 'filled', ...
    'MarkerFaceAlpha', 0.5, 'MarkerEdgeAlpha', 0.5);
end
rline = refline(1, 0);
rline.Color = 'k';
rline.LineStyle = '--';
rline.LineWidth = 2;
xlim([min(min(cp_newl)) max(max(cp_newl))]);
ylim([min(min(cp_oldl)) max(max(cp_oldl))]);

legend(cellstr(strcat(num2str(wl'), ' nm')), 'location', 'southeast')
title('c_p in the Benguela current ~8 days average N = 13672')
xlabel('c_p new light ACS57 [m^{-1}]')
ylabel('c_p old light ACS348 [m^{-1}]')
set(gca, 'FontSize', 12)



%% %%%%%%%%%%%% TARA Pacific %%%%%%%%%%%% %%
cd('/Users/gui/Documents/MATLAB/InLineAnalysis/InLineAnalysis-master/')

% Load InLineAnalysis and the configuration
ila = InLineAnalysis('cfg/TaraPacific_cfg.m');

%%
ila.cfg.days2run = datenum(2016,6,8,1,30,0):datenum(2016,6,8,14,57,0);
ila.cfg.instruments2run = {'FLOW','ACS007'};
ila.cfg.qcref.view = 'ACS007';
ila.cfg.parallel = Inf;
%% 
ila.Read('prod');
%% get N
ila.instrument.ACS007.prod.p = ila.instrument.ACS007.prod.p(1:N_newlight,:);
N_oldlight = size(ila.instrument.ACS007.prod.p, 1);
%%
ila.DiagnosticPlot('AC',{'prod'}); % AC or BB
%% extract data for old light
% id676 = abs(ila.instrument.ACS007.lambda_a - 676) == min(abs(ila.instrument.ACS007.lambda_a - 676));
lambda_a_oldlight = ila.instrument.ACS007.lambda_a;
% ap_oldlight = ila.instrument.ACS007.prod.p.ap - ila.instrument.ACS007.prod.p.ap(:,id676);
% ap_sd_oldlight = ila.instrument.ACS007.prod.p.ap_sd - ila.instrument.ACS007.prod.p.ap_sd(:,id676);
% ap_oldlight = ila.instrument.ACS007.prod.p.ap ./ ila.instrument.ACS007.prod.p.ap(:,id676);
% ap_sd_oldlight = ila.instrument.ACS007.prod.p.ap_sd ./ ila.instrument.ACS007.prod.p.ap_sd(:,id676);
ap_oldlight = ila.instrument.ACS007.prod.p.ap;
ap_sd_oldlight = ila.instrument.ACS007.prod.p.ap_sd;
ap_sigtonoi_oldlight = mean(ila.instrument.ACS007.prod.p.ap_sd ./ ...
  ila.instrument.ACS007.prod.p.ap, 'omitnan') * 100;
ap_sigtonoi_sd_oldlight = std(ila.instrument.ACS007.prod.p.ap_sd ./ ...
  ila.instrument.ACS007.prod.p.ap, 'omitnan') * 100;

%% compare signal to noise ratios
% minute binned ap spectra averaged over 2 hours
figure()
subplot(1,2,1)
shadedErrorBar(lambda_a_newlight, mean(ap_newlight, 'omitnan'), std(ap_newlight, 'omitnan'), ...
  'lineProps', '-r', 'transparent', true)
shadedErrorBar(lambda_a_oldlight, mean(ap_oldlight, 'omitnan'), std(ap_oldlight, 'omitnan'), ...
  'lineProps', '-k', 'transparent', true)
ylabel('ap [m^{-1}]')
xlabel('lambda')
title('Minute binned ap average over ~8h (equatorial Atlantic) N = 496')
legend('New light (Tara Microbiome 2021)', 'Old light (Tara Pacific 2016)')


% minute binned ap_sd spectra averaged over 2 hours
subplot(1,2,2)
shadedErrorBar(lambda_a_newlight, mean(ap_sd_newlight, 'omitnan'), std(ap_sd_newlight, 'omitnan'), ...
  'lineProps', '-r', 'transparent', true)
shadedErrorBar(lambda_a_oldlight, mean(ap_sd_oldlight, 'omitnan'), std(ap_sd_oldlight, 'omitnan'), ...
  'lineProps', '-k', 'transparent', true)
ylabel('ap std [m^{-1}]')
xlabel('lambda')
title('Minute binned ap\_sd average over ~8h (equatorial Atlantic) N = 496')

% % minute binned ap_sd spectra averaged over 2 hours
% subplot(1,3,3)
% shadedErrorBar(lambda_a_newlight, sigtonoi_newlight, sigtonoi_sd_newlight, 'lineProps', '-r', 'transparent', true)
% shadedErrorBar(lambda_a_oldlight, sigtonoi_oldlight, sigtonoi_sd_oldlight, 'lineProps', '-k', 'transparent', true)
% ylabel('ap/ap_sd [unitless]')
% xlabel('lambda')
% title('signal to noise of Minute binned ap\_sdap/ average over 2h (equatorial Atlantic)')






