%% Cruise
% IOP package
clear
addpath utils/
addpath utils/ac_meters/
addpath cruise
addpath ../utils/
addpath ../utils/data_ie/
addpath ../utils/misc_lab/
data_path = 'cruise/data/';

%% 1. Load data
update=false;
if update
  for i=[1:22 36 38:43];
    fprintf('Loading cast %03d... ', i);
    [acs(i).jd, acs(i).wl, acs(i).a, acs(i).c] =...
      ac_import_profile(sprintf('%srun_21_ACS.%03d', data_path, i));
    [bb9(i).jd, bb9(i).wl, bb9(i).bb9] =...
      bb9_import_profile(sprintf('%srun_22_ASCII.%03d', data_path, i)); 
    [bb3(i).jd, bb3(i).scat_angle, bb3(i).bb3] =...
      ecoVSF_import_profile(sprintf('%srun_23_ASCII.%03d', data_path, i));
    [ctd(i).jd, ctd(i).p, ctd(i).t, ctd(i).s, ctd(i).sig] =...
      ctd_import_profile(sprintf('%srun_24_ASCII.%03d', data_path, i));
    [fchl(i).jd, fchl(i).fchl, fchl(i).fchl_delta] =...
      fchl_import_profile(sprintf('%srun_29_ANALOG.%03d', data_path, i)); % WS3S
    fprintf('Done\n');
  end;
  clearvars('-except', 'acs', 'bb9', 'bb3', 'ctd', 'fchl');
  save('cruise/iop_l0');
else
  load('cruise/iop_l0');
end;

clearvars('-except', 'acs', 'bb9', 'bb3', 'ctd', 'fchl');

%% 2. Merge with CTD according to field jd
sel = [1:22 36 38:43];
% file 6 ACS has no time 
for i=sel;
  % Meging IOPs
  iop1(i) = merge_profiles(ctd(i), acs(i), 0.000008);
  iop2(i) = merge_profiles(iop1(i), bb9(i), 0.000008);
  iop3(i) = merge_profiles(iop2(i), bb3(i), 0.000008);
  iop(i) = merge_profiles(iop3(i), fchl(i), 0.000008);
end;

% Adding wavelength
meta.wl_acs = acs(i).wl;
meta.wl_bb9 = bb9(i).wl;
meta.scat_angle = bb3(i).scat_angle;

clearvars('-except', 'acs', 'bb9', 'bb3', 'ctd', 'fchl', 'iop', 'meta', 'sel');

%% 3. Synchronise visualization
fp=fig(1, 'sync');
i = 3;
n = iop(i).jd;

% Normalize
foo=iop(i).sig(:,1); sig = (foo-min(foo))/max(foo-min(foo));
foo=iop(i).c(:,10); acs = (foo-min(foo))/max(foo-min(foo)); % 3.5e-5
foo=iop(i).bb9(:,5); bb9 = (foo-min(foo))/max(foo-min(foo)); % 2.4e-5
foo=iop(i).bb3(:,2); bb3 = (foo-min(foo))/max(foo-min(foo)); % 2.9e-5
foo=iop(i).fchl(:,1); fchl = (foo-min(foo))/max(foo-min(foo));

% Select var
var1 = sig;
var2 = bb9;
%var2(acs > 0.9) = NaN; 

% plot
plot(n, var1, n, var2);

% UI
fui=fig(2, 'controls');
% shift
global shift_; shift_=0;
global scale_; scale_=1;
ts = 0.000001;
b1 = uicontrol('Parent',fui,'Style','slider','Position',[81,54,419,23],...
              'value',0, 'min',-40*ts, 'max',40*ts, 'SliderStep',[1,1]/160);
uicontrol('Parent',fui,'Style','text','Position',[50,54,23,23],...
                'String','-20','BackgroundColor',fui.Color);
uicontrol('Parent',fui,'Style','text','Position',[500,54,23,23],...
                'String','20','BackgroundColor',fui.Color);
uicontrol('Parent',fui,'Style','text','Position',[240,35,100,23],...
                'String','Shift','BackgroundColor',fui.Color);
set(b1,'Callback', @(es,ed) callback(n, var1, var2, 'shift', -es.Value, fp));
b2 = uicontrol('Parent',fui,'Style','slider','Position',[81,104,419,23],...
              'value',-1, 'min',-20, 'max',2, 'SliderStep',[1,1]/40);
uicontrol('Parent',fui,'Style','text','Position',[50,104,23,23],...
                'String','-2','BackgroundColor',fui.Color);
uicontrol('Parent',fui,'Style','text','Position',[500,104,23,23],...
                'String','2','BackgroundColor',fui.Color);
uicontrol('Parent',fui,'Style','text','Position',[240,85,100,23],...
                'String','Scale','BackgroundColor',fui.Color);
set(b2,'Callback', @(es,ed) callback(n, var1, var2, 'scale', -es.Value, fp));

clearvars('-except', 'acs', 'bb9', 'bb3', 'ctd', 'fchl', 'iop', 'meta', 'sel');

%% Shift values
% RUN ONLY ONCE
shift.bb9_filter = -2.7*10^(-5);
shift.bb9 = -1.4*10^(-5);
shift.bb3_filter = -2.7*10^(-5);
shift.bb3 = -1.4*10^(-5);
shift.acs_filter = -1.3*10^(-5);
shift.acs = -0.4*10^(-5);
shift.fchl_filter = -0.4*10^(-5);
shift.fchl = -0.2*10^(-5);
for i=sel;
  if i == 1 || i == 2 || i == 7 || i == 8 || i == 36
    % filtered
    s_bb9=[]; s_bb3=[]; s_acs=[]; s_fchl=[];
    for j = 1:size(iop(i).jd, 1)
      foo = find(abs((iop(i).jd(j) - shift.bb9_filter) - iop(i).jd) < 0.000005);
      if ~isempty(foo); s_bb9(j) = j-foo; end;
      foo = find(abs((iop(i).jd(j) - shift.bb3_filter) - iop(i).jd) < 0.000005);
      if ~isempty(foo); s_bb3(j) = j-foo; end;
      foo = find(abs((iop(i).jd(j) - shift.acs_filter) - iop(i).jd) < 0.000005);
      if ~isempty(foo); s_acs(j) = j-foo; end;
      foo = find(abs((iop(i).jd(j) - shift.fchl_filter) - iop(i).jd) < 0.000005);
      if ~isempty(foo); s_fchl(j) = j-foo; end;
    end;
  else
    % non filtered
    s_bb9=[]; s_bb3=[]; s_acs=[]; s_fchl=[];
    for j = 1:size(iop(i).jd, 1)
      foo = find(abs((iop(i).jd(j) - shift.bb9) - iop(i).jd) < 0.000005);
      if ~isempty(foo); s_bb9(j) = j-foo; end;
      foo = find(abs((iop(i).jd(j) - shift.bb3) - iop(i).jd) < 0.000005);
      if ~isempty(foo); s_bb3(j) = j-foo; end;
      foo = find(abs((iop(i).jd(j) - shift.acs) - iop(i).jd) < 0.000005);
      if ~isempty(foo); s_acs(j) = j-foo; end;
      foo = find(abs((iop(i).jd(j) - shift.fchl) - iop(i).jd) < 0.000005);
      if ~isempty(foo); s_fchl(j) = j-foo; end;
    end;
  end;
  % Apply correction
  s_bb9 = floor(median(s_bb9));
  s_bb3 = floor(median(s_bb3));
  s_acs = floor(median(s_acs));
  s_fchl = floor(median(s_fchl));
  iop(i).bb9 = [iop(i).bb9(1-s_bb9:end, :); nan(-s_bb9, size(iop(i).bb9,2))];
  iop(i).bb3 = [iop(i).bb3(1-s_bb3:end, :); nan(-s_bb3, size(iop(i).bb3,2))];
  iop(i).a = [iop(i).a(1-s_acs:end, :); nan(-s_acs, size(iop(i).a,2))];
  iop(i).c = [iop(i).c(1-s_acs:end, :); nan(-s_acs, size(iop(i).c,2))];
  iop(i).fchl = [iop(i).fchl(1-s_fchl:end, :); nan(-s_fchl, size(iop(i).fchl,2))];
  iop(i).fchl_delta = [iop(i).fchl_delta(1-s_fchl:end, :); nan(-s_fchl, size(iop(i).fchl_delta,2))];
end;

clearvars('-except', 'acs', 'bb9', 'bb3', 'ctd', 'fchl', 'iop', 'meta', 'sel');

%% 4. Cut profiles visualization
fp = fig(1, 'p');
i=36;
n=1:size(iop(i).jd,1);
plot(n, iop(i).p);
cut(1).b = 155; cut(1).e = 440;
cut(2).b = 7; cut(2).e = 278;
cut(3).b = 38; cut(3).e = 264;
cut(4).b = 52; cut(4).e = 247;
cut(5).b = 31; cut(5).e = 238;
cut(6).b = 75; cut(6).e = 524;
cut(7).b = 82; cut(7).e = 186; % Bug arroud 10m
cut(8).b = 4; cut(8).e = 88;
cut(9).b = 27; cut(9).e = 110;
cut(10).b = 3; cut(10).e = 49; % Got a lobster trap
cut(11).b = 88; cut(11).e = 376;
cut(12).b = 9; cut(12).e = 264;
cut(13).b = 1; cut(13).e = 58;
cut(14).b = 11; cut(14).e = 297;
cut(15).b = 23; cut(15).e = 364;
cut(16).b = 1; cut(16).e = 106;
cut(17).b = 45; cut(17).e = 143;
cut(18).b = 58; cut(18).e = 136;
cut(19).b = 41; cut(19).e = 122;
cut(20).b = 93; cut(20).e = 163;
cut(21).b = 46; cut(21).e = 102;
cut(22).b = 4; cut(22).e = 54;

%% Apply correction
% RUN ONLY ONCE
for i=sel;
  fields = fieldnames(iop);
  for f=1:size(fields,1);
    %if ~strcmp(fields{f}, 'wl_bb9') && ~strcmp(fields{f}, 'wl_acs') && ~strcmp(fields{f}, 'scat_angle');
    if size(cut,2) >= i;
      iop(i).(fields{f}) = iop(i).(fields{f})(cut(i).b:cut(i).e, :);
    end;
    %end;
  end;
end;

clearvars('-except', 'acs', 'bb9', 'bb3', 'ctd', 'fchl', 'iop', 'meta', 'sel');

%% Save data L1a
%info='Data has been loaded from ascii files, c for acs is interpolated according to the wavelength of a, data has been merged and synchronize according to the CTD, data has been cut to have only downward profile, profile 6 has no c tubes on acs, caps are on bb9 and ws3s and descent rate vary';
save('cruise/iop_l1a', 'iop', 'meta');
for i=sel;
  % Adding wavelength for csv struct
  iop(i).wl_acs = acs(i).wl;
  iop(i).wl_bb9 = bb9(i).wl;
  iop(i).scat_angle = bb3(i).scat_angle;
end;
struct2csv(iop, 'cruise/iop_l1a.csv');

%% Load data
clearvars('-except', 'acs', 'bb9', 'bb3', 'ctd', 'fchl', 'iop', 'meta', 'sel');
load('cruise/iop_l1a');

%% 5. Plot raw data
fig(1, 'raw'); hold('all');
i=36;
ColorSet = varycolor(size(iop(i).p,1));
set(gca, 'ColorOrder', ColorSet); h=[];
for j=1:size(iop(i).p,1);
  h(j)=plot(meta.wl_acs, iop(i).c(j,:));
end;
xlabel('wavelength (nm)');
ylabel('c (m^{-1})');
title(sprintf('Profile %02d from IOP Package', i));
legend([h(1) h(floor(end/2)) h(end)], sprintf('%1.2f m', iop(i).p(1)),...
  sprintf('%2.2f m', iop(i).p(floor(end/2))),...
  sprintf('%2.2f m', iop(i).p(end)), 'Location', 'NorthWest');

% %% 
% set(gcf,'PaperUnits','centimeters','PaperPosition',[0 0 20 10]);
% print('-dpng', '-r300', 'cruise/img/bb9_darks');

%% 5. Level 1B: depth bin every meter
fields = fieldnames(iop(1));
iop_1b=struct();
for i=sel;
  p_p = 1:floor(max(iop(i).p));
  k=1;
  for p=p_p
    iop_1b(i).p(k,1) = p;
    foo = find(abs(iop(i).p - p) < 0.5);
    for f=1:size(fields,1);
      if ~strcmp(fields{f}, 'p') && ~strcmp(fields{f}, 'jd') && ~strcmp(fields{f}, 'wl_acs') && ~strcmp(fields{f}, 'wl_acs') && ~strcmp(fields(f), 'wl_bb9') && ~strcmp(fields{f}, 'scat_angle');
        if ~isempty(foo)
          if size(foo,1) > 1;
            iop_1b(i).(fields{f})(k,:) = nanmedian(iop(i).(fields{f})(foo, :));
          elseif size(foo,1) == 1
            iop_1b(i).(fields{f})(k,:) = iop(i).(fields{f})(foo, :);
          else
            iop_1b(i).(fields{f})(k,:) = nan(size(iop(i).(fields{f}),2),1);
          end;
        else
          iop_1b(i).(fields{f})(k,:) = nan(size(iop(i).(fields{f}),2),1);
        end;
      end;
    end;
    k = k + 1;
  end;
end;

%% Save Level 1 B
iop = iop_1b;
save('cruise/iop_l1b', 'iop', 'meta');
for i=sel;
  % Adding wavelength for csv struct
  iop(i).wl_acs = acs(i).wl;
  iop(i).wl_bb9 = bb9(i).wl;
  iop(i).scat_angle = bb3(i).scat_angle;
end;
struct2csv(iop, 'cruise/iop_l1b.csv');

%% Load Level 1 B
clearvars('-except', 'acs', 'bb9', 'bb3', 'ctd', 'fchl', 'iop', 'meta', 'sel');
load('cruise/iop_l1b');

%% 6. Corrections L2
%% 6.1 BB9
i_dark = 5;
% raw
fig(1,'raw'); 
plot(iop(i_dark).bb9(:,3),iop(i_dark).p);
% Scale factor
scale_factor_bb9 = [2.262 2.167 1.975 1.767 1.622 1.294 1.083 0.9036 0.6703] * 10^(-5);
% compute darks
for j=1:size(iop(i_dark).bb9,2)
  dark_bb9_med(j) = nanmedian(iop(i_dark).bb9(:,j));
  dark_bb9_prc(j,:) = prctile(iop(i_dark).bb9(:,j),[5 95]);
end;
% remove darks & aply scale factor
for i=sel;
  iop(i).bb9_d = iop(i).bb9 - dark_bb9_med(ones(size(iop(i).bb9,1),1),:);
  iop(i).bb9_sf = iop(i).bb9_d .* scale_factor_bb9(ones(size(iop(i).bb9,1),1),:);
end;
%%
fig(2, 'dark_p'); i=i_dark;
plot(iop(i).bb9, iop(i).p); axis('ij');
xlabel('b_b (counts)'); ylabel('pressure (dBar)');
title('BB9-132 dark vs pressure'); ylim([0 93]);
legend('412 nm', '440 nm', '488 nm', '510 nm', '532 nm', '595 nm', '650 nm', '715 nm', '880 nm', 'Location', 'SouthWest')
%%
fig(3, 'dark_t'); i=i_dark;
plot(iop(i).bb9, iop(i).t);
xlabel('b_b (counts)'); ylabel('temperature ({\circ}C)');
title('BB9-132 dark vs temperature'); ylim([5.8 16]);
legend('412 nm', '440 nm', '488 nm', '510 nm', '532 nm', '595 nm', '650 nm', '715 nm', '880 nm', 'Location', 'NorthWest')
%%
fig(4,'dark_prc');
errorbar(meta.wl_bb9,dark_bb9_med, dark_bb9_prc(:,1) - dark_bb9_prc(:,2) );
ylabel('b_b (counts)'); xlabel('wavelength (nm)');
title('BB9-132 darks');
%%
fig(3, 'sf');i=2;
plot(iop(i).bb9_sf, iop(i).p);
axis('ij');
xlabel('b_b (m^{-1})'); ylabel('Pressure (dBar)');
title('Monday - Station 1 - Ocean');
legend('412 nm', '440 nm', '488 nm', '510 nm', '532 nm', '595 nm', '650 nm', '715 nm', '880 nm', 'Location', 'NorthEast')
set(gca, 'FontSize', 18);

%% 
% set(gcf,'PaperUnits','centimeters','PaperPosition',[0 0 30 20]);
% print('-dpng', '-r300', 'cruise/img/ws3s_npq');

%% 6.2 WS3S (FCHL)
dark_ws3s = 0.078; % 0.051 according to manufacturer
                   % 0.078 for median(sort(iop(2).fchl)(1:10));
                   % 0.077 min iop(
scale_factor_ws3s = 14.8;
for i=sel;
  iop(i).fchl_d = iop(i).fchl - dark_ws3s;
  iop(i).fchl_sf = iop(i).fchl_d .* scale_factor_ws3s;
end;

fig(2, 'dark'); i=4;
plot(iop(i).fchl, iop(i).p, 'b', iop(i).fchl_d, iop(i).p, 'r'); axis('ij');
fig(3, 'sf'); hold('on');
i=3; plot(iop(i).fchl_sf, iop(i).p);
i=4; plot(iop(i).fchl_sf, iop(i).p);
i=9; plot(iop(i).fchl_sf, iop(i).p, '--');
i=10; plot(iop(i).fchl_sf, iop(i).p, '--');
axis('ij');
xlabel('Chlorophyll {\ita} fluorescence (mg.m^{-3})'); ylabel('Pressure (dBar)');
title('Monday - Station 1 and 2');
legend('station 1 cast 3', 'station 1 cast 4','station 2 cast 9','station 2 cast 10', 'Location', 'SouthEast');
set(gca, 'FontSize', 18);
%% 6.3 ecoVSF (BB3)
scale_factor_bb3 = [4.9932E-05 5.1319E-05 2.9128E-05];
% compute darks
for j=1:size(iop(i_dark).bb3,2)
  dark_bb3_med(j) = nanmedian(iop(i_dark).bb3(:,j));
  dark_bb3_prc(j,:) = prctile(iop(i_dark).bb3(:,j),[5 95]);
end;
% remove darks & aply scale factor
for i=sel;
  iop(i).bb3_d = iop(i).bb3 - dark_bb3_med(ones(size(iop(i).bb9,1),1),:);
  iop(i).bb3_sf = iop(i).bb3_d .* scale_factor_bb3(ones(size(iop(i).bb9,1),1),:);
end;
%%
fig(2, 'dark_p'); i=i_dark;
plot(iop(i).bb3, iop(i).p); axis('ij');
xlabel('b_b (counts)'); ylabel('pressure (dBar)');
title('ECO-VSF-004R dark vs pressure'); %ylim([0 93]);
%%
fig(3, 'dark_t'); i=i_dark;
plot(iop(i).bb3, iop(i).t);
xlabel('b_b (counts)'); ylabel('temperature ({\circ}C)');
title('ECO-VSF-004R dark vs temperature'); %ylim([5.8 16]);
%%
fig(4,'dark_prc');
errorbar(meta.scat_angle,dark_bb3_med, dark_bb3_prc(:,1) - dark_bb3_prc(:,2) );
ylabel('b_b (counts)'); xlabel('scattering angle (deg)');
title('ECO-VSF-004R darks');
%%
fig(3, 'sf');
plot(iop(i).bb9_sf, iop(i).p);

%% Keep only good variables in struct
for i=sel;
  iop_l2(i).p = iop(i).p;
  iop_l2(i).t = iop(i).t;
  iop_l2(i).s = iop(i).s;
  iop_l2(i).sig = iop(i).sig;
  iop_l2(i).a = iop(i).a;
  iop_l2(i).c = iop(i).c;
  iop_l2(i).bb9 = iop(i).bb9_sf;
  iop_l2(i).bb3 = iop(i).bb3_sf;
  iop_l2(i).fchl = iop(i).fchl_d;
end;

%% Save Level 2
iop = iop_l2;
save('cruise/iop_l2', 'iop', 'meta');
for i=sel;
  % Adding wavelength for csv struct
  iop(i).wl_acs = meta.wl_acs;
  iop(i).wl_bb9 = meta.wl_bb9;
  iop(i).scat_angle = meta.scat_angle;
end;
struct2csv(iop, 'cruise/iop_l2.csv');

%% Correction Level 3
load('cruise/iop_l2');
%% ACS
% Day 1, Station 1: cast 4 - cast 2
% Day 1, Station 2: cast 10 - cast 8
% Day 2, Station 3: cast 12 - cast 15
% profile 13 scrapped
% Day 2, Station 4: cast 17 or 18 - cast 19 or 22
% Day 3, Station 5: cast [38:43 - 36];
i_f = [2 8  12 18 36*ones(1,6)]; % filtered
i_t = [4 10 15 19 38:43]; % total
i_s = 1:10; % Stations id
for i=1:10;
  n = min(size(iop(i_f(i)).p,1), size(iop(i_t(i)).p,1));
  a = iop(i_t(i)).a(1:n,:) - iop(i_f(i)).a(1:n,:);
  c = iop(i_t(i)).c(1:n,:) - iop(i_f(i)).c(1:n,:);
  delta_t = iop(i_t(i)).t(1:n,:) - iop(i_f(i)).t(1:n,:); 
  delta_s = iop(i_t(i)).s(1:n,:) - iop(i_f(i)).s(1:n,:);
  [a_ts, c_ts] = ac_correct_ts(meta.wl_acs, meta.wl_acs, a, c, delta_t, delta_s);
  [a_sc] = ac_correct_scatter(meta.wl_acs, a_ts, c_ts);

  iop_l3(i_s(i)).p = iop(i_t(i)).p(1:n);
  iop_l3(i_s(i)).t = median([iop(i_t(i)).t(1:n), iop(i_f(i)).t(1:n)], 2);
  iop_l3(i_s(i)).s = median([iop(i_t(i)).s(1:n), iop(i_f(i)).s(1:n)], 2);
  iop_l3(i_s(i)).sig = median([iop(i_t(i)).sig(1:n), iop(i_f(i)).sig(1:n)], 2);
  iop_l3(i_s(i)).a = a_sc;
  iop_l3(i_s(i)).c = c_ts;
end;

%% Plot corrected ACs
fig(1, 'raw'); hold('all');
i_s=2;
ColorSet = varycolor(size(iop_l3(i_s).p,1));
set(gca, 'ColorOrder', ColorSet); h=[];
for j=1:size(iop_l3(i_s).c,1);
  %h(j)=plot(meta.wl_acs, iop_l3(i_s).c(j,:) - iop_l3(i_s).a(j,:));
  h(j)=plot(meta.wl_acs, iop_l3(i_s).a(j,:));
end;
xlabel('wavelength (nm)');
ylabel('b_p (m^{-1})');
title(sprintf('Station %d - ACS', i_s));
legend_info = {};
for i=1:10%[1 20:20:80 90];
  legend_info{end+1} = sprintf('%2.0f m', i);
end;
legend(h([1:10]),legend_info);
%% Correction ECO-BB9 and ECO_VSF Level 3
i_stn = [4, 10, 15, 19, 38:43];
% Remove SW effect
for i=1:size(i_stn,2);
  [iop_l3(i).bbp_BB9, iop_l3(i).beta_p_BB9] = compute_bbp(iop(i_stn(i)).bb9, iop(i_stn(i)).t, iop(i_stn(i)).s, meta.wl_bb9, 117);
  [iop_l3(i).bbp_VSF, iop_l3(i).beta_p_VSF] = compute_bbp(iop(i_stn(i)).bb3, iop(i_stn(i)).t, iop(i_stn(i)).s, 650, meta.scat_angle);
end;

%% WS3S corrections NPQ + scaling factor corrections
scale_factor_ws3s = 14.8;
for i=1:size(i_stn,2);
  iop_l3(i).fchl = iop(i_stn(i)).fchl * scale_factor_ws3s;
  iop_l3(i).mld = estimate_mld(iop(i_stn(i)).p, iop(i_stn(i)).sig, 'fixed_density', 0.005);
  iop_l3(i).fchl_xing = correct_npq(iop_l3(i).fchl, iop_l3(i).p, iop_l3(i).mld);
end;

%% Export Level 3
iop = iop_l3;
save('cruise/iop_l3', 'iop', 'meta');
for i=1:10;
  % Adding wavelength for csv struct
  iop(i).wl_acs = meta.wl_acs;
  iop(i).wl_bb9 = meta.wl_bb9;
  iop(i).scat_angle = meta.scat_angle;
end;
struct2csv(iop, 'cruise/iop_l3.csv');

