%%%% Open Optics package process acs data

clear all
close all

mkdir Plots %%%%%% make folder called Plots that the plots go into.
depthplt = 18; %%% plot data from XX m 
printplots = 1;
exclude_upcast = 1;

%%% these are the wavelengths for acs 025. Serial #53000019
lambda_c = [401	405.1	409.3	413.5	417.8	422.8	427.4	431.8	436.3	440.9	445.6	450.9	455.5	460.1	464.8	469.9	474.8	480.3	485.3	490.1	494.8	499.6	504.7	509.7	515.1	520.6	525.4	530.4	535.3	540.6	545.4	550.5	555.8	561	565.9	570.8	575.6	579.8	584.4	588.8	593.7	597.7	602.3	607	612.3	617.5	622	626.5	631.4	636	640.7	645.3	650.1	654.7	659.3	664.1	668.5	672.9	677.3	681.7	685.9	689.8	693.6	697.6	701.7	705.4	708.9	712.6	716	720	723.2	726.2	729.7	732.7	735.9	738.9	741.5	744.2	746.8];
lambda_a = [402.1	406.8	410.5	415	419.6	424.4	429	433.8	438.2	442.3	447.8	452.9	457.6	462	466.9	471.7	476.9	482.4	487.3	491.8	496.7	501.5	506.9	512.1	517.5	522.7	527.9	532.9	537.5	542.7	547.8	552.9	558.1	563.3	568.1	572.9	577.3	582	586.4	591	595.7	600.3	604.8	608.8	613.9	618.5	623.4	628.4	632.9	637.3	642.1	647.1	651.6	656.5	661.1	665.2	670.1	674.4	678.8	683	687.1	691.5	695	699.1	703	706.7	710.5	714	717.3	720.9	724.1	727.8	730.8	734	736.9	739.7	742.4	745.6	748.1];
lambda = (lambda_a+lambda_c)./2;



% filenameCast = 'Cruise2_IOP_cast3_BB9_to_bb.csv';
% filenameCast = 'Cruise2_IOP_cast2_BB9_to_bb.csv';
 filenameCast = 'Cruise2_IOP_cast1_BB9_to_bb.csv';


filenameDarkFiltCast = 'Cruise2_IOP_dark2.csv';


%%%% Lab MQ correction. These files have the headers, wavelengths, time column and extra data (temp etc) on the
%%%% far right columns removed. 
filenameBlanka_noheader = 'acs025_Cruise2_CleanedMQ_a_228_nohead.txt';
filenameBlankc_noheader = 'acs025_Cruise2_CleanedMQ_c_227_nohead.txt';


CastImport = importdata(filenameCast);
DarkFiltCastImport = importdata(filenameDarkFiltCast);



%%% 1) eliminate duplicates at every 1 second interval in the raw data

[Castraw] = acs_elim_duplicates(CastImport.data);
[DarkFiltCastraw] = acs_elim_duplicates(DarkFiltCastImport.data);

[~,depthcut] = max(Castraw(:,5));
[~,depthcutDark] = max(DarkFiltCastraw(:,5));

isdown = 1:size(Castraw,1) < depthcut;
isdownDark = 1:size(DarkFiltCastraw,1) < depthcutDark;





%%% 2) Find the time lag of the flow through the package
goodCast = Castraw(:,5)>7;%isdown;
a740 = Castraw(goodCast,end-3);  %%%% wavelength = 740 nm
Temp = Castraw(goodCast,3); %% ctd temperature (C)
Sal= Castraw(goodCast,6);
Time = Castraw(goodCast,1); %% time from start (ms)
z = -Castraw(goodCast,5); %% depth from pressure (m)

a740Dark = DarkFiltCastraw(:,end-3);  %%%% wavelength = 740 nm
c740Dark = DarkFiltCastraw(:,end-3-length(lambda));  %%%% wavelength = 740 nm
TempDark = DarkFiltCastraw(:,3); %% ctd temperature (C)
TimeDark = DarkFiltCastraw(:,1); %% time from start (ms)
zDark = -DarkFiltCastraw(:,5); %% depth from pressure (m)


%%%%% find optimal time lag in the filtered Data through a lagged
%%%%% correlation analysis
[dtbestDark r2bestDark trangeDark r2rangeDark] = acs_findtimelag(a740Dark, TempDark, TimeDark);
[dtbestDarkc r2bestDarkc trangeDarkc r2rangeDarkc] = acs_findtimelag(c740Dark, TempDark, TimeDark);

acs_time_correctedDark = TimeDark - dtbestDark;
acs_time_correctedDarkc = TimeDark - dtbestDarkc;
%  shift time
a740Dark_shifted = interp1(acs_time_correctedDark,a740Dark,TimeDark);
c740Dark_shifted = interp1(acs_time_correctedDarkc,c740Dark,TimeDark);


a_Dark_shifted = interp1(acs_time_correctedDark,...
    DarkFiltCastraw(:,end-length(lambda)+1:end),DarkFiltCastraw(:,1));
c_Dark_shifted = interp1(acs_time_correctedDarkc,...
    DarkFiltCastraw(:,end-length(lambda)*2+1:end-length(lambda)),DarkFiltCastraw(:,1));

araw = Castraw(:,end-length(lambda)+1:end);
craw = Castraw(:,end-length(lambda)*2+1:end-length(lambda));

Temp = Castraw(:,3); %% ctd temperature (C)
Sal= Castraw(:,6);
Time = Castraw(:,1); %% time from start (ms)
z = -Castraw(:,5); %% depth from pressure (m)

TempDark = DarkFiltCastraw(:,3); %% ctd temperature (C)
SalDark= DarkFiltCastraw(:,6);
TimeDark = DarkFiltCastraw(:,1); %% time from start (ms)
zDark = -DarkFiltCastraw(:,5); %% depth from pressure (m)

[ablank cblank ablank_er cblank_er] = ...
acs_loadMQblanks(filenameBlanka_noheader, filenameBlankc_noheader);

%%%%%%


%%%%%%%%%%%%% MQ water temperature correction
Ta = 22.7;
Tc = 22.6;
[ablankT, cblankT] = acs_ts_correct_1(ablank,  cblank, 20.7, Ta, Tc, 0, lambda_a,lambda_c);



findz = find(Castraw(:,5)>depthplt); %%%% find depth (depthplt) in the unbinned files to plot
fz = findz(1);
figure (1);
subplot 211
hold on
plot(lambda_a, ablank, 'k', 'linewidth', 2)
plot(lambda_a, ablankT, 'r', 'linewidth', 2)
plot([400 750], [0 0], 'k--')
ylabel('a_{MQ} (m^{-1})')
set(gca, 'xticklabel', [])
title('Milli-Q blanks')
legend('Raw', 'TS-corr')
grid
subplot 212
hold on
plot(lambda_a, cblank, 'k', 'linewidth', 2)
plot(lambda_c, cblankT, 'r', 'linewidth', 2)
plot([400 750], [0 0], 'k--')
ylabel('c_{MQ} (m^{-1})')
xlabel('\lambda (nm)')
legend('Raw', 'TS-corr')
clean_up_figure_text
grid
if printplots == 1
orient tall
cd Plots
print('-dpdf','Fig1_MQ')
cd ..
end



%%% SAL AND TEMP CORRECTION DATA - do this initially to get a gross correction
%%% then after binning, do the residual correction.
ld = length(z);
for tt = 1:ld
    [arawT(tt,:), crawT(tt,:)] = acs_ts_correct_1(araw(tt,:),  craw(tt,:), ...
        21, Temp(tt), Temp(tt), Sal(tt),lambda_a,lambda_c);    
end


%%% SAL AND TEMP CORRECTION FILTERED DATA - do this initially to get a gross correction
%%% then after binning, do the residual correction.
ldz = length(zDark);
TshiftDarka = 0.2;
TshiftDarkc = -1;
for tt = 1:ldz
    [a_Dark_shiftedT(tt,:), c_Dark_shiftedT(tt,:)] = acs_ts_correct_1(a_Dark_shifted(tt,:),  ...
        c_Dark_shifted(tt,:), 21, TempDark(tt)-TshiftDarka, TempDark(tt)-TshiftDarkc, SalDark(tt),lambda_a,lambda_c);
end


% % %% Subtract Blank

a = arawT - ones(ld,1) * ablankT;
c = crawT - ones(ld,1) * cblankT;

ag = a_Dark_shiftedT - ones(ldz,1) * ablankT;
cg = c_Dark_shiftedT - ones(ldz,1) * cblankT;


figure (2);
subplot 211
hold on
plot(lambda_a, araw(fz,:), 'k', 'linewidth', 2)
plot(lambda_a, arawT(fz,:), 'b', 'linewidth', 2)
plot(lambda_a, a(fz,:), 'r', 'linewidth', 2)
plot([400 750], [0 0], 'k--')
ylabel('a (m^{-1})')
set(gca, 'xticklabel', [])
title('Unfiltered ACS data')
legend('Raw', 'TS-corr','blank-corr')
grid
subplot 212
hold on
plot(lambda_c, craw(fz,:), 'k', 'linewidth', 2)
plot(lambda_c, crawT(fz,:), 'b', 'linewidth', 2)
plot(lambda_c, c(fz,:), 'r', 'linewidth', 2)
plot([400 750], [0 0], 'k--')
ylabel('c (m^{-1})')
xlabel('\lambda (nm)')
legend('Raw', 'TS-corr','blank-corr')
clean_up_figure_text
grid
if printplots == 1
orient tall
cd Plots
print('-dpdf','Fig2_Unfiltered')
cd ..
end

figure (3);
subplot 211
hold on
plot(lambda_a, a_Dark_shifted(fz,:), 'k', 'linewidth', 2)
plot(lambda_a, a_Dark_shiftedT(fz,:), 'b', 'linewidth', 2)
plot(lambda_a, ag(fz,:), 'r', 'linewidth', 2)
plot([400 750], [0 0], 'k--')
ylabel('a_{g} (m^{-1})')
set(gca, 'xticklabel', [])
title('Filtered ACS data')
legend('Raw', 'TS-corr', 'blank-corr')
grid
subplot 212
hold on
plot(lambda_c, c_Dark_shifted(fz,:), 'k', 'linewidth', 2)
plot(lambda_c, c_Dark_shiftedT(fz,:), 'b', 'linewidth', 2)
plot(lambda_c, cg(fz,:), 'r', 'linewidth', 2)
plot([400 750], [0 0], 'k--')
ylabel('c_{g} (m^{-1})')
xlabel('\lambda (nm)')
legend('Raw', 'TS-corr', 'blank-corr')
clean_up_figure_text
grid
if printplots == 1
orient tall
cd Plots
print('-dpdf','Fig3_Filtered')
cd ..
end

%%%%%%%%%%%%%%%%%%%%%% interpolate onto lambda (as a and c channels have
%%%%%%%%%%%%%%%%%%%%%% slightly different wavelengths

% %% uncomment if we want to interpolate onto the new wavelengths
ai = interp1(lambda_a,a', lambda,'linear', 'extrap')';
ci = interp1(lambda_c, c', lambda,'linear','extrap')';
agi = interp1(lambda_a,ag', lambda,'linear','extrap')';
cgi = interp1(lambda_c, cg', lambda,'linear','extrap')';




% %% Bin 

%%%% average into XXm bins.
XX = 1;
zbin = -(1:XX:80);
zh = -(XX/2:XX:80+XX/2);

%%%% zbin2 = actualy mean depth of each binned measurement.
%%%% zbin = uniformly spaced grid of depths. This assumes that the data are
%%%% normally distributed around the bin center.

for i = 1:length(zbin);
    findz = find(z<zh(i) & z>zh(i+1) & isdown');
    if length(findz)>0
        abin(i,:) = nanmedian(ai(findz,:));
        cbin(i,:) = nanmedian(ci(findz,:)); 
        abinstd(i,:) = nanstd(ai(findz,:));
        cbinstd(i,:) = nanstd(ci(findz,:)); 
        nbin(i) = length(findz);
        zbin2(i,:) = nanmedian(z(findz));
       
    else
        abin(i,:) = nanmedian(ai(findz,:))+NaN;
        cbin(i,:) = nanmedian(ci(findz,:))+NaN;
        abinstd(i,:) = nanstd(ai(findz,:))+NaN;
        cbinstd(i,:) = nanstd(ci(findz,:))+NaN;
        nbin(i) = 0;
        zbin2(i,:) = nanmean(z(findz))+NaN;
       
    end
    
    findz = find(zDark<zh(i) & zDark>zh(i+1) & isdownDark');
    if length(findz)>0
        agbin(i,:) = nanmedian(agi(findz,:));
        cgbin(i,:) = nanmedian(cgi(findz,:));
        agbinstd(i,:) = nanstd(agi(findz,:));
        cgbinstd(i,:) = nanstd(cgi(findz,:));
        ngbin(i) = length(findz);
    else
        agbin(i,:) = nanmedian(agi(findz,:))+NaN;
        cgbin(i,:) = nanmedian(cgi(findz,:))+NaN;    
        agbinstd(i,:) = nanstd(agi(findz,:))+NaN;
        cgbinstd(i,:) = nanstd(cgi(findz,:))+NaN;    
        ngbin(i) = 0;
    end
end




bbin = cbin - abin;
bbinstd = sqrt(cbinstd.^2 + abinstd.^2);

apbin = abin - agbin;   
apbinstd = sqrt(abinstd.^2 + agbinstd.^2);

cpbin = apbin + bbin;

ap_uncorr = apbin;
a_uncorr = abin;
ag_uncorr = agbin;
cp = cpbin;
bp_uncorr = cpbin - apbin;

wavel = lambda;


%%%%% Now do a residual correction for temperature. 
[ap_TSalScatCorr, deltaT] = ACS_ResidTempScatCorr(ap_uncorr, cpbin, wavel);
apbin = ap_TSalScatCorr;

figure (7)
plot(deltaT, zbin)
findnan = find(isnan(deltaT) == 1);
zbin(findnan)=[];

[fitdT] = polyfit(zbin', deltaT, 1);

%%% there is an offset in the ag in the infrared. We find the median in
%%% this range and subtract
i = find(lambda > 705);
agbin_new = agbin - nanmedian(agbin(:,i)')'*ones(1,79);

%%% calculate a from ap and ag
abin = apbin + agbin_new;

bbin = cbin-abin;
bpbin = cpbin-apbin;



figure (4);
subplot 211
hold on
plot(lambda, ap_uncorr(depthplt,:), 'k', 'linewidth', 2)
plot(lambda, ap_TSalScatCorr(depthplt,:), 'r', 'linewidth', 2)
plot([400 750], [0 0], 'k--')
ylabel('a_p (m^{-1})')
set(gca, 'xticklabel', [])
title('Residual and Scatter Correct ACS data')
legend('TS-corr', 'Res/Scatt-corr')
grid
subplot 212
hold on
plot(lambda, bp_uncorr(depthplt,:), 'k', 'linewidth', 2)
plot(lambda, bpbin(depthplt,:), 'r', 'linewidth', 2)
plot([400 750], [0 0], 'k--')
ylabel('c_{p} (m^{-1})')
xlabel('\lambda (nm)')
legend('TS-corr', 'Res/Scatt-corr')
clean_up_figure_text
grid
if printplots == 1
orient tall
cd Plots
print('-dpdf','Fig4_ResScatCorr')
cd ..
end



% 
% %%% print some examples of T correcting cg
% 
% 
% break
% for ff = 1:4
% TshiftDarkc = -1;
% for tt = 1:ldz
%     [~, c_Dark_shiftedT(tt,:)] = acs_ts_correct_1(a_Dark_shifted(tt,:),  ...
%         c_Dark_shifted(tt,:), 21, TempDark(tt)-TshiftDarka, TempDark(tt)-TshiftDarkc, SalDark(tt),lambda_a,lambda_c);
% end
% 
% 
% figure (5);
% hold on
% plot(lambda, ap_uncorr(depthplt,:), 'k', 'linewidth', 2)
% plot(lambda, ap_TSalScatCorr(depthplt,:), 'r', 'linewidth', 2)
% plot([400 750], [0 0], 'k--')
% ylabel('a_p (m^{-1})')
% set(gca, 'xticklabel', [])
% title('Residual and Scatter Correct ACS data')
% legend('TS-corr', 'Res/Scatt-corr')
% grid
% if printplots == 1
% cd Plots
% print('-dpdf','Fig4_ResScatCorr')
% cd ..
% end
% 
% end
