%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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

dateday='20210505';
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=18.2;


direct=['/Users/emmanuel.boss/Desktop/EXPORTS02/SlowDROP/',dateday,'/'];
ACS_file='run_21_05_ACS.004';
CTD_file='run_23_05_CTD-ENGR.004';


yr=dateday(1:4);
mo=dateday(5:6);
day=dateday(7:8);
mat_date_IOP=datenum(str2num(yr),str2num(mo),str2num(day));
    %start with dissolved profile:
[time, P, T, C, S, V1, V2, D]=rd_SBE19([direct,CTD_file]);

figure(1)
plot(time,-P)
[X ,Y]=ginput(2)
    
    %insures profile start at highest point.
    %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);
S=S(I);
    
    
%load ac-s data
[t,dat,wla,wlc,itemp,xtemp]=rd_wap_acs111([direct,ACS_file]);
time_acs=t(1:length(t))-6.5/24/3600; %assume a 6.5s delay for AC-S due to filter

%median filter the ac9 data with a 5 sample median filter  to remove spikes (1.25 sec,approximately matching CTD):
for j=1:164
    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:5])=[P(I2),T(I2),S(I2),V1(I2)];
for j=1:164
    dat1(:,6+j)=interp1(time_acs,acs_new(:,j),time_merge,'linear');
end

%average decent rate:
disp(['average decent rate of package was ',num2str(100*mean(diff(dat1(:,2)))), 'cm/s'])

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

%choose calibrations
load /Users/emmanuel.boss/Desktop/EXPORTS02/SlowDROP/cals_SlowDROP/cal_111
I=find(abs(date_a-mat_date_IOP)<=2);

for i=1:length(I)
    for j=1:length(prrg)
        [datout_a(j,:), datout2_a(j,:), dS, dT]=TS_corr_ACS(datg(j,89:164)-cal_a(I(i),1:82-6), 'a', wla(1:82-6), Sg, Tg, Tcal);
        amp0=datout_a(j,10);
        s0=0.015;
        III=find(wla<560 & wla >430);
        %[x1] = fminsearch(@least_squares_ag, [amp0 s0], opts, datout_a(j,III), wla(III));
        %misfit(j)=sum(abs(datout_a(j,III)-x1(1)*exp(-x1(2)*(wla(III)-440))));
        misfit(j)=sum(abs(datout_a(j,III)-amp0*exp(-s0*(wla(III)-440))));
%         plot(wla(III),datout_a(j,III),'k',wla(III),x1(1)*exp(-x1(2)*(wla(III)-440)),':');
%         pause
    end
    sum_misfit(i)=sum(misfit);
end
II=find(sum_misfit==min(sum_misfit));
cal_a_best=cal_a(I(min(II)),:);
disp(['best a-cal was: cal',num2str(I(min(II))),datestr(date_a(I(min(II))))])
    
[ag_cal_a, ag_cal2_a, dS, dT]=TS_corr_ACS(datg(:,89:164-6)-cal_a_best(1:82-6-6), 'a',  wla(1:82-6-6), Sg, Tg,   Tcal);
wla=wla(1:82-6-6);
    
clear sum_misfit;

I=find(abs(date_c-mat_date_IOP)<=2);
    
for i=1:length(I)
    for j=1:length(prrg)
        [datout_c(j,:), datout2_c(j,:), dS, dT]=TS_corr_ACS(datg(j,7:88-6)-cal_c(I(i),7:88-6), 'c', wlc(1:82-6), 0, Tcal, Tcal);
        amp0=datout_c(j,10);
        s0=0.015;
        III=find(wlc<560 & wlc >430);
%         [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))));
        misfit(j)=sum(abs(datout_a(j,III)-amp0*exp(-s0*(wla(III)-440))));

    end
    sum_misfit(i)=sum(misfit);
end
II=find(sum_misfit==min(sum_misfit));
cal_c_best=cal_c(I(min(II)),:);
disp(['best c-cal was: cal',num2str(I(min(II))),datestr(date_c(I(min(II))))])    
[cg_cal_c, cg_cal2_c, adS, dT]=TS_corr_ACS(datg(:,7:88-12)-cal_c_best(1:82-6-6), 'c', wlc(1:82-6-6), Sg, Tg,   Tcal);
wlc=wlc(1:82-6-6);

outfile=['dis',yr,mo,day,'.mat']
    
figure(1);
plot(wla,ag_cal2_a, 'b',wlc,cg_cal2_c, 'k')
title(datestr(mat_date_IOP))
xlabel('wavelength [nm]')
ylabel('a_g (b) and c_g (k) [m^-1');
%ylabel('a_g (b) [m^-1');
export_fig(['dis',yr,mo,day,'_ac_g.png'])
   
% figure(2);
%  title(datestr(mat_date_IOP(ii)))
%  plot(ag_cal_a(:,15),-prrg,'b',FDOMg,-prrg,'--')
%  export_fig([ACS_files(1).name(8:9),yr,mo,day,ACS_files(ii).name(20:22),'CDOM_a_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
%     
%     II=input('which do you prefer ag (1) cg (2) or their average (3)?   ')
%     if II==1
%          save(outfile, 'wla', 'ag_cal_a', 'prrg');;
%     elseif II==2
%          save(outfile, 'wlc', 'ag_cal_c', 'prrg');;
%     elseif II==3
%         ag_cal_c=interp1(wlc,ag_cal_c',wla,'linear','extrap');
%         ag_average=(ag_cal_c'+ag_cal_a)/2;
%         save(outfile, 'wla', 'ag_average', 'prrg');;
%     end    
    save(outfile, 'wla', 'ag_cal_a', 'prrg', 'wlc', 'cg_cal_c');

