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

%  ACS='301';
%  range0=[1:3];
%  range1=[1:164];
%  range2=[87:168];
%  range3=[5:86];
%  ACS='024'
%  range0=[4:5];
%  range1=[1:168];
%  range2=[89:172];
%  range3=[5:88];
  ACS='298' 
%  range0=[6:7]; %for E1
  %range0=[1:8];%for E2
range0=[1:2];%for E3
  range1=[1:168];
  range2=[89:172];
  range3=[5:88];
  Ext_tot=['014';'020'; '024'];%no timestamp
  Ext_dis=['012';'018'; '022'];%no timestamp
%    Ext_tot=['046';'003';'009'];% 009 with no timestamp
%    Ext_dis=['043';'001';'007'];%007 with no timestamp
% Ext_tot=['001';'005'; '009';'015';'021'; '028';'034';'041'];
% Ext_dis=['003';'007'; '013';'019';'025'; '031';'037';'039'];
%  Ext_tot=['031';'040'];
%  Ext_dis=['037';'046';];
% Ext_tot=['020';'025'];
% Ext_dis=['023';'027';];
% Ext_tot=['006';'010';'012'];
% Ext_dis=['004';'008';'017'];

%instrument specific information
Tcal=17.8;

%read total profiles
direct='/Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Total_casts/E3/'
% Tot_ACS_files=dir([direct,'run_E1_21*']);
% Tot_CTD_files=dir([direct,'run_E1_24*']);
% Dis_ACS_files=dir([direct,'dis_E1_21*']);
%Dis_CTD_files=dir([direct,'dis_E1_24*']);
% Tot_ACS_files=dir([direct,'run_E2_21*']);
% Tot_CTD_files=dir([direct,'run_E2_24*']);
% Dis_ACS_files=dir([direct,'dis_E2_21*']); 
% Dis_CTD_files=dir([direct,'dis_E2_24*']);
%Tot_ACS_files=dir([direct,'run_E3_21*']);
%Tot_CTD_files=dir([direct,'run_E3_24*']);
% Dis_ACS_files=dir([direct,'dis_E3_21*']);
% Dis_CTD_files=dir([direct,'dis_E3_24*']);


for ii=1:3%1:length(Tot_ACS_files)
    Tot_ACS_files(ii,:)=dir([direct,'arc_E3_21_ACS*.',Ext_tot(ii,:)]);
    Tot_CTD_files(ii,:)=dir([direct,'arc_E3_24_CTD*.',Ext_tot(ii,:)]);
    Dis_ACS_files(ii,:)=dir([direct,'dis_E3_21_ACS*.',Ext_dis(ii,:)]);
    Dis_CTD_files(ii,:)=dir([direct,'dis_E3_24_CTD*.',Ext_dis(ii,:)]);

    %start with total profile:
    [time, P1, T1, S1, C]=rd_SB37([direct,Tot_CTD_files(ii).name]);
    time=[0:1:length(time)-1]*0.0000157409995154012;
    time=time';
    
    figure(1)
    plot(time,-P1)
    [X ,Y]=ginput(2)
    
    %insures profile start at highest point.
    %tar pressure to 1db at surface
    dp=min(P1)-1;
    I=find(time>=X(1) & time <=X(2));
    time_ctd=time(I);
    P1=P1(I)+dp;
    T1=T1(I);
    C=C(I);
    
    %compute average salinity based on a CTD regression to T:
    SS1=zeros(1,length(T1))';
    I=find(T1>=6);
    SS1(I)=32.3+(14-T1(I))/32;
    I=find(T1<6);
    SS1(I)=32.55 +0.9*(6-T1(I));;
    
     
    %load ac-s data

%    [t,dat,wla,wlc,itemp,xtemp]=rd_WAP_acs301([direct,Tot_ACS_files(ii).name]);
    [t,dat,wla,wlc,itemp,xtemp]=rd_WAP_acs298([direct,Tot_ACS_files(ii).name]);
    t=[0:1:length(t)-1]*0.25/24/3600;
%    [t,dat,wla,wlc,itemp,xtemp]=rd_WAP_acs024([direct,Tot_ACS_files(ii).name]);
    a_nir=medfilt1(dat(:,164),20);
    %figuring out the delay:
    for i=1:200
           delta_t(i)=0.2*(i-100)/24/3600;
           NIR=interp1(t-delta_t(i),a_nir,time_ctd,'linear');
           I=find(isfinite(NIR) & T1 >11 );
           P=polyfit(NIR(I),T1(I),1);
           res(i)=sum(abs(T1(I)-polyval(P,NIR(I))));
           %plot(T1-polyval(P,NIR),-P1)
    end
    I=find(res==min(res));
    disp(['delay particulate =  ',num2str(delta_t(I)*3600*24)])

    
    time_acs=t-delta_t(I); 
    
    %median filter the ac9 data with a 5 sample median filter  to remove spikes (1.25 sec,approximately matching CTD):
    for j=range1
        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:1:4])=[P1(I2),T1(I2),SS1(I2)];
    
    
    for j=range1
        dat1(:,4+j)=interp1(time_acs,acs_new(:,j),time_merge,'linear');
    end
   
    [datT, err, prrT]=binpre(dat1',dat1(:,2),ceil(min(dat1(:,2))),floor(max(dat1(:,2))),1,2);
    datT=datT';
    err=err';
    
        %filter out spectra with jumps in their middle

    flag=[];
    for i=1:length(prrT)
        if max(diff(datT(i,30:50)))<0.002; %good spectra
          flag=[flag i];
        end
    end
    if length(flag)<length(prrT)
        disp('lost spectra');
        disp(length(prrT)-length(flag));
        datT=datT(flag,:);
        err=err(flag,:);
        prrT=prrT(flag);
    end
    
    clear acs_new time dat1 P1 T1 SS1 res delta_t NIR t C flag

    %start with dissolved profile:
    [time, P2, T2, S2, C]=rd_SB37([direct,Dis_CTD_files(ii).name]);
    time=[0:1:length(time)-1]*0.0000157409995154012;
    time=time';

    
    figure(1)
    plot(time,-P2)
    [X ,Y]=ginput(2)
    
    %insures profile start at highest point.
    %tar pressure to 1db at surface
    dp=min(P2)-1;
    I=find(time>=X(1) & time <=X(2));
    time_ctd=time(I);
    P2=P2(I)+dp;
    T2=T2(I);
    C=C(I);
    
    %compute average salinity based on a CTD regression to T:
    SS2=zeros(1,length(T2))';
    I=find(T2>=6);
    SS2(I)=32.3+(14-T2(2))/32;
    I=find(T2<6);
    SS2(I)=32.55 +0.9*(6-T2(I));;
    
    %align depth between two casts based on temperature.
  
    for i=1:1:60
        delta(i)=i*0.1-3;
        T22=interp1(P2+delta(i),T2,P2,'linear');
        I=find(isfinite(T22));
        diffT(i)=sum(T22(I)-T2(I));
    end
    I=find(abs(diffT)==min(abs(diffT)));
    P2=P2+delta(I);
      
    %load ac-s data
%   [t,dat,wla,wlc,itemp,xtemp]=rd_WAP_acs301([direct,Dis_ACS_files(ii).name]);
%   t=[0:1:length(t)-1]*0.25/24/3600;
%
    [t,dat,wla,wlc,itemp,xtemp]=rd_WAP_acs298([direct,Dis_ACS_files(ii).name]);
    t=[0:1:length(t)-1]*0.25/24/3600;
%    [t,dat,wla,wlc,itemp,xtemp]=rd_WAP_acs024([direct,Dis_ACS_files(ii).name]);
    a_nir=medfilt1(dat(:,164),20);
    %figuring out the delay:
    for i=1:400
           delta_t(i)=0.2*(i-100)/24/3600;
           NIR=interp1(t-delta_t(i),a_nir,time_ctd,'linear');
           I=find(isfinite(NIR) & T2 >11 & P2 >20 & P2<40);
           P=polyfit(NIR(I),T2(I),1);
           res(i)=sum(abs(T2(I)-polyval(P,NIR(I))));
           %plot(T1-polyval(P,NIR),-P1)
    end
    I=find(res==min(res));

    I=find(res==min(res));
    disp(['delay dissolved =  ',num2str(delta_t(I)*3600*24)])

    time_acs=t(1:length(t))-delta_t(I); 
    
    %median filter the ac9 data with a 5 sample median filter  to remove spikes (1.25 sec,approximately matching CTD):
    for j=range1
        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
    dat2(:,1)=(time_merge-min(time_merge))*24*3600; %time in secs'
    dat2(:,[2:4])=[P2(I2),T2(I2),SS2(I2)];
    for j=range1
        dat2(:,4+j)=interp1(time_acs,acs_new(:,j),time_merge,'linear');
    end
   
    [datg,err2,prrg]=binpre(dat2',dat2(:,2),ceil(min(dat2(:,2))),floor(max(dat2(:,2))),1,1); %1meter median bins
    datg=datg';
    err2=err2';
    
    %filter out spectra with jumps in their middle
    k=0
    flag=[];
    for i=1:length(prrg)
        if max(diff(datg(i,30:50)))<0.002 ; %good spectra
          flag=[flag i];
        end
    end
    if length(flag)<length(prrg)
        disp('lost spectra');
        disp(length(prrg)-length(flag));
        datg=datg(flag,:);
        err2=err2(flag,:);
        prrg=prrg(flag);
    end
            
    
    [C,IA,IB] = intersect(prrT,prrg);
    P=prrT(IA);
    T=datT(IA,3);
    S=datT(IA,4);
    dat(1:length(C),1:4)=datT(IA,1:4);
    dat(:,5:4+max(range1))=datT(IA,5:4+max(range1))-datg(IB,5:4+max(range1));
    unc_cp=sqrt(err(IA,range3).^2+err2(IB,range3).^2);
    unc_ap=sqrt(err(IA,range2).^2+err2(IB,range2).^2);
    ap_nc=dat(:,range2);
    cp=dat(:,range3);
    for i=1:length(P)
         if (max(diff(ap_nc(i,30:50)))>0.002 | max(unc_ap(i,20:80)>0.005)); %bad spectra
             ap_nc(i,:)=NaN;
             unc_ap(i,:)=NaN;
         end
         if (max(diff(cp(i,30:50)))>0.002 | max(unc_cp(i,20:80)>0.005)); %bad spectra
             cp(i,:)=NaN;
             unc_cp(i,:)=NaN;
         end
    end

    ap_nc2=interp1(wla,ap_nc',wlc,'linear')'; %interpolate a's into c's wavelenths.
    [ap_new,cp_new,wl] = ACS_ResidTempScatCorr(ap_nc2, cp, wlc);   
    ap=ap_new;
    cp=cp_new;
    unc_ap=interp1(wla,unc_ap',wl,'linear'); 
    unc_cp=interp1(wlc,unc_cp',wl,'linear'); 
    
    figure(2),plot(wl,cp_new);
    ylabel('c_p(\lambda)[m^-1]')
    xlabel('wavelength [nm]')
%    export_fig(['cp_',Tot_ACS_files(1).name(5:6),'_',num2str(ii),'.png'])
    export_fig(['cp_',Tot_ACS_files(ii).name(5:6),'_',num2str(ii+3),'.png'])
    figure(3),plot(wl,ap_new);    
    ylabel('a_p (\lambda)[m^-1]')
    xlabel('wavelength [nm]') 
%    export_fig(['ap_',Tot_ACS_files(1).name(5:6),'_',num2str(ii),'.png'])
    export_fig(['ap_',Tot_ACS_files(ii).name(5:6),'_',num2str(ii+3),'.png'])
    %figure(4),plot(P,cp(:,50),'k',P,ap_nc2(:,15),'b');
   


%    outfile=['part',Tot_ACS_files(1).name(5:6),'_',num2str(ii),'.mat']
    outfile=['part',Tot_ACS_files(ii).name(5:6),'_',num2str(ii+3),'.mat']
    save(outfile, 'P', 'T', 'S' ,'wl', 'ap', 'cp', 'unc_ap', 'unc_cp');
    clear acs_new time dat2 P1 T1 SS1 P2 T2 SS2 dat1 dat2 datg datT prrg prrT cp ap res NIR delta_t unc flag
 
end




