%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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
% specific to ACS-301
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=23.0;
load RR_CTD_data.mat

direct='/Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Dissolved_casts/'
ACS_files=dir([direct,'run21__E1*']);
CTD_files=dir([direct,'run24__E1*']);
FDOM_files=dir([direct,'run23__E1*']);
clear A
for ii=1:length(ACS_files)
    yr=ACS_files(ii).name(11:14);
    mo=ACS_files(ii).name(15:16);
    day=ACS_files(ii).name(17:18);
    mat_date_IOP(ii)=datenum(str2num(yr),str2num(mo),str2num(day));
    %start with dissolved profile:
    [time, P, T, S, C]=rd_SB37([direct,CTD_files(ii).name]);
    %find closest CTD cast
    I=find(abs(mat_date_IOP(ii)-mat_time)==min(abs(mat_date_IOP(ii)-mat_time)));
    CTD_Temp=Temp(I,:)
    CTD_Sal=Sal(I,:)
    CTD_depth=[1:1:110];
    
    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);
    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);
    % use a cruise long T-S relation to get salinity
      I=find(T>=6)
      SS(I)=32.3+(14-T(I))/32;
      I=find(T<6)
      SS(I)=32.55 +0.9*(6-T(I))
      SS=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_acs301([direct,ACS_files(ii).name]);
    time_acs=t(1:length(t))-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:164
        acs_new(:,j)=medfilt1(dat(:,j),5);
    end
    clear dat

    [time FDOM]=rd_Eco1ch([direct,FDOM_files(ii).name]);
    time_FDOM=time;
    FDOM=filtfilt([1 1 1 1 1 1]/6,1,(FDOM-80)/1000); %reduce noise by averaging
    
    %determine begining and end time for merging.
    I=find(time_ctd>=time_acs(1) & time_ctd>=time_FDOM(1));
    I1=find(time_ctd<=max(time_acs) & time_ctd<=max(time_FDOM));
    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:164
        dat1(:,4+j)=interp1(time_acs,acs_new(:,j),time_merge,'linear');
    end
   
    FDOM=interp1(time_FDOM,FDOM,time_merge,'linear');
    %bin data to 3m bins by averaging
    [datg,prrg]=binpre(dat1',dat1(:,2),ceil(min(dat1(:,2))),floor(max(dat1(:,2))),1,2);
    datg=datg';
    [FDOMg,prrg]=binpre(FDOM',dat1(:,2),ceil(min(dat1(:,2))),floor(max(dat1(:,2))),1,2);
   [Tg,prrg]=binpre(T',dat1(:,2),ceil(min(dat1(:,2))),floor(max(dat1(:,2))),1,2);
   [Sg,prrg]=binpre(SS',dat1(:,2),ceil(min(dat1(:,2))),floor(max(dat1(:,2))),1,2);
   
    %choose calibrations
    load /Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Calibrations/cal_301
    I=find(abs(cal_date_a-mat_date_IOP(ii))<=2);
    
    for i=1:length(I)
        for j=1:length(prrg)
            [datout_a(j,:), dS, dT]=TS_corr_ACS(datg(j,87:168)-cal_a(I(i),:), 'a', wla, 0, Tcal, 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
    II=find(sum_misfit==min(sum_misfit));
    cal_a_best=cal_a(I(II),:);
    disp(['best a-cal was: cal',num2str(I(II)),datestr(cal_date_a(I(II)))])
    
    [ag_cal_a, dS, dT]=TS_corr_ACS(datg(:,87:168)-cal_a_best, 'a', wla, datg(:,4), datg(:,3),   Tcal);
    
    clear sum_misfit

%     I=find(abs(cal_date_c-mat_date_IOP(ii))<=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(i),:), 'c', wlc, 0, Tcal, 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
%     II=find(sum_misfit==min(sum_misfit));
%     cal_c_best=cal_c(I(II),:);
%     disp(['best c-cal was: cal',num2str(I(II)),datestr(cal_date_c(I(II)))])    
%     [ag_cal_c, dS, dT]=TS_corr_ACS(datg(:,5:88)-cal_c_best, 'c', wlc, mean_Salinity, datg(:,2),   Tcal);
%     
    outfile=['dis',ACS_files(1).name(8:9),yr,mo,day,ACS_files(ii).name(24:26),'.mat']
    
    figure(1);
    %plot(wla,ag_cal_a, 'b',wlc,ag_cal_c, 'k')
    plot(wla,ag_cal_a, 'b')
    title(datestr(mat_date_IOP(ii)))
    xlabel('wavelength [nm]')
    %ylabel('a_g (b) and c_g (k) [m^-1');
    ylabel('a_g (b) [m^-1');
    export_fig([ACS_files(1).name(8:9),yr,mo,day,ACS_files(ii).name(24:26),'_a_g.png'])
   
%    figure(3);plot(ag_cal_a(:,15),-prrg,'b',ag_cal_c(:,15),-prrg,'k',FDOMg,-prrg,'--')
    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(24:26),'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', 'Tg', 'Sg');
    clear dat1 dat t acs_new ag_cal_c ag_cal_a datg prrg sum_misfit Tg Sg P SS
end

