%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This program is designed to merge data collected by multiple instruments with the DH-4
%All data is merged to CTD data, the slowest sampling instrument.
%applies lags to instrument data
clear all
close all
opts = optimset('fminsearch');      
opts = optimset(opts,'MaxIter',20000000); 
opts = optimset(opts,'MaxFunEvals',20000);
opts = optimset(opts,'TolX',1e-8);
opts = optimset(opts,'TolFun',1e-8);


%instrument specific information
Tcal=17.8;

%day specific information
fnam_dis='013';
fnam_tot='009';
date=datenum(2018,08,25);
title(datestr(date))
direct='/Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Epoch2Day3/'

%start with dissolved profile:
[time, P, T, S, C]=rd_SB37([direct,'run_24_CTD-ENGR.',fnam_dis]);

figure(1)
plot(time,-P)
[X ,Y]=ginput(2)

%insures profile start at highest point.
II=find(P==min(P));
if (time(II)>X(1))
    X(1)=time(II);
end
%tar pressure to 1db at surface
dp=min(P)-1; 
I=find(time>=X(1) & time <=X(2));
time_ctd=time(I);
P=P(I)+dp;
T=T(I);
C=C(I);

%compute average salinity for salinity correction of a_g:
time_conductivity=time_ctd+16/(24*3600); %delay between pumped conductivity and temperature
CC=interp1(time_conductivity,C,time_ctd,'linear');
I=find(isfinite(CC));
SS=sw_salt(CC(I)/4.2914,T(I),P(I));
mean_Salinity=nanmean(SS);

%look at settling speed statistics
figure(1)
dt=median(diff(time_ctd))*24*3600;
hist(diff(P)/dt)
xlabel('settling speed of package')
disp('hit enter to continue')
pause

%load ac-s data
[t,dat,wla,wlc,itemp,xtemp]=rd_WAP_acs298([direct,'run_21_ACS.',fnam_dis]);
time_acs=t(1:length(dat))-6.5/24/3600; %assume a 6.5s delay for ac-9

%median filter the ac9 data with a 5 sample median filter  to remove spikes (1.25 sec,approximately matching CTD):
for j=1:168
 acs_new(:,j)=medfilt1(dat(:,j),5);
end
clear dat 

%determine begining and end time for merging.
I=find(time_ctd>=time_acs(1));
I1=find(time_ctd<=max(time_acs));
I2=find(time_ctd>=time_ctd(min(I)) & time_ctd<=time_ctd(max(I1)));
time_merge=time_ctd(I2); 

%create interpolated data matrix
dat1(:,1)=(time_merge-min(time_merge))*24*3600; %time in secs'
dat1(:,[2:3])=[P(I2),T(I2)];
for j=1:168
    dat1(:,4+j)=interp1(time_acs,acs_new(:,j),time_merge,'linear');
end

%bin data to 3m bins by averaging
[datg,prrg]=binpre(dat1',dat1(:,2),ceil(min(dat1(:,2))),floor(max(dat1(:,2))),3,2);
datg=datg';

%choose calibrations
load /Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Calibrations/cal_298
I=find(abs(cal_date_a-date)<=2);

for i=1:length(I)
    for j=1:length(prrg)
      [datout_a(j,:), dS, dT]=TS_corr_ACS(datg(j,89:172)-cal_a(i,:), 'a', wla, mean_Salinity, datg(j,2), Tcal);
      amp0=datout_a(j,10);
      s0=0.015;
      [x1] = fminsearch(@least_squares_ag, [amp0 s0], opts, datout_a(j,:), wla);
      misfit(j)=sum(abs(datout_a(j,:)-x1(1)*exp(-x1(2)*(wla-440))));
    end
   sum_misfit(i)=sum(misfit);
end
I=find(sum_misfit==min(sum_misfit));
cal_a_best=cal_a(I,:);
disp(['best a-cal was: cal',I,cal_date_a(I)])

[ag_cal_a, dS, dT]=TS_corr_ACS(datg(:,89:172)-cal_a_best, 'a', wla, mean_Salinity, datg(:,2),   Tcal);


I=find(abs(cal_date_c-date)<=2);

for i=1:length(I)
    for j=1:length(prrg)
      [datout_c(j,:), dS, dT]=TS_corr_ACS(datg(j,5:88)-cal_c(i,:), 'c', wlc, mean_Salinity, datg(j,2), Tcal);
      amp0=datout_a(j,10);
      s0=0.015;
      [x1] = fminsearch(@least_squares_ag, [amp0 s0], opts, datout_c(j,:), wlc);
      misfit(j)=sum(abs(datout_c(j,:)-x1(1)*exp(-x1(2)*(wlc-440))));
    end
   sum_misfit(i)=sum(misfit);
end
I=find(sum_misfit==min(sum_misfit));
cal_c_best=cal_c(I,:);
disp(['best c-cal was: cal',I,cal_date_c(I)])

[ag_cal_c, dS, dT]=TS_corr_ACS(datg(:,5:88)-cal_c_best, 'c', wlc, mean_Salinity, datg(:,2),   Tcal);


figure;
plot(wla,ag_cal_a, 'b',wlc,ag_cal_c, 'k')
title(datestr(date))
xlabel('wavelength [nm]')
ylabel('a_g (b) and c_g (k) [m^-1])
export_fig([fnam,'_a_g_c_g.png'])

figure;
plot(wla,cal_a,'b-',wla,cal_a_best,'k-')
figure;
plot(wlc,cal_c,'b-',wlc,cal_c_best,'k-')

pause

clear dat1 dat t acs_new

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%analyze total data


%reading data in and getting their sampling time
[time, P, T, S, C]=rd_SB37([direct,'run_24_CTD-ENGR.',fnam_tot]);
figure(3)
plot(time,-P)
[X ,Y]=ginput(2)

%insures profile start at highest point.
II=find(P==min(P));
if (time(II)>X(1))
    X(1)=time(II);
end
%tar pressure to 1db at surface
dp=min(P)-1; 
I=find(time>=X(1) & time <=X(2));
time_ctd=time(I);
P=P(I)+dp;
T=T(I);
C=C(I);


[t,dat,wla,wlc,itemp,xtemp]=rd_WAP_acs298([direct,'run_21_ACS.',fnam_tot]);
time_acs=t(1:length(dat))-6.5/24/3600; %assume a 6.5s delay for ac-9

%median filter the ac9 data with a 5 sample median  filter (1.25 sec):
for j=1:168
 acs_new(:,j)=medfilt1(dat(:,j),5);
end

%determine begining and end time for merging.
I=find(time_ctd>=time_acs(1));
I1=find(time_ctd<=max(time_acs));
I2=find(time_ctd>=time_ctd(min(I)) & time_ctd<=time_ctd(max(I1)));
time_merge=time_ctd(I2); 

%create interpolated data matrix
dat1(:,1)=(time_merge-min(time_merge))*24*3600; %time in secs'
dat1(:,[2:3])=[P(I2),T(I2)];
for j=1:168
    dat1(:,4+j)=interp1(time_acs,acs_new(:,j),time_merge,'linear');
end


%bin data to 1m bins
[datpg,prrpg]=binpre(dat1',dat1(:,2),ceil(min(dat1(:,2))),floor(max(dat1(:,2))),3,2);
datpg=datpg';

%apply calibrations
load /Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Calibrations/cal_298_20180820
datpg(:,5:172)=datpg(:,5:172)-[calc cala];

figure(3)
subplot(1,2,1)
plot(datg(:,[10:10:60]),-datg(:,2),'-k')
hold on
plot(datpg(:,[10:10:60]),-datpg(:,2),'-b')
hold off
title('c_g, c_{pg}')
subplot(1,2,2)
plot(datg(:,[90:10:150]),-datg(:,2),'-k')
hold on
plot(datpg(:,[90:10:150]),-datpg(:,2),'-b')
hold off
title('a_g, a_{pg}')

[C,IA,IB] = intersect(prrg,prrpg)
ap=datpg(IB,89:172)-datg(IA,89:172);
cp=datpg(IB,5:88)-datg(IA,5:88);
figure(4)
plot(wla,ap)
figure(5)
plot(wla,cp)

