%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This program is designed to process the c-side of the ac-s used in
%EXPORTS01 to obtain the total attenuation coefficient.
%calibration is chosen based on providing the most realism of the spectrum.

clear all
close all

%Global parameters
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='024';
range0=[4:5];
range1=[1:168];
range2=[20:84];
%instrument specific information
Tcal=17.8;

%read total profiles
direct='/Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Total_casts/'
Tot_ACS_files=dir([direct,'run_E1_21*']);
Tot_CTD_files=dir([direct,'run_E1_24*']);
Tot_bb9_files=dir([direct,'run_E1_22*']);

for ii=range0 %length(Tot_ACS_files)
    yr=Tot_ACS_files(ii).name(15:18);
    mo=Tot_ACS_files(ii).name(19:20);
    day=Tot_ACS_files(ii).name(21:22);
    mat_date_IOP(ii)=datenum(str2num(yr),str2num(mo),str2num(day));
   [time, P, T, S, C]=rd_SB37([direct,Tot_CTD_files(ii).name]);
    
    %select range of data for the profile
    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)-8/24/3600;
    P=P(I)+dp;
    T=T(I);
    C=C(I);
    
    %compute average salinity based on a CTD regression to T:
      SS=zeros(1,length(T))';
      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));
     
    %load ac-s data
    [t,dat,wla,wlc,itemp,xtemp]=rd_WAP_acs024([direct,Tot_ACS_files(ii).name]);
     time_acs=t(1:length(t))-6.5/24/3600; %assume a 6.5s delay for ac-9
    %filter bad spectra
    %I=find((abs(dat(:,46)-dat(:,41))./((dat(:,46)+dat(:,41))))>0.04);
    I=find((abs(dat(:,46)-dat(:,41)))<0.017);
    time_acs=time_acs(I);
    dat=dat(I,:);
    
    if ii ~= 4
        [time wl BB9]=rd_BB9counts([direct,Tot_bb9_files(ii).name]);
        time_bb9=time;
    end
    
   
    %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
    %median filter BB9 with a 3-pt median filter (3s)
    if ii ~= 4
        for j=1:9
            BB9_new(:,j)=medfilt1(BB9(:,j),3);
        end
        clear BB9
    end


    %determine begining and end time for merging
    if ii ~= 4
        I=find(time_ctd>=time_acs(1) & time_ctd>=time_bb9(1));
        I1=find(time_ctd<=max(time_acs) & time_ctd<=max(time_bb9));
        I2=find(time_ctd>=time_ctd(min(I)) & time_ctd<=time_ctd(max(I1)));
        time_merge=time_ctd(I2);
    else
        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);
    end
    %create interpolated data matrix
    dat1(:,1)=(time_merge-min(time_merge))*24*3600; %time in secs'
    dat1(:,[2:1:4])=[P(I2),T(I2),SS(I2)];
    if ii ~= 4
        for j=1:9
            dat1(:,4+j)=interp1(time_bb9,BB9_new(:,j),time_merge,'linear');
        end
        for j=range1
            dat1(:,13+j)=interp1(time_acs,acs_new(:,j),time_merge,'linear');
        end
    else
        for j=range1
            dat1(:,4+j)=interp1(time_acs,acs_new(:,j),time_merge,'linear');
        end
    end        
    [datT,prrT]=binpre(dat1',dat1(:,2),ceil(min(dat1(:,2))),floor(max(dat1(:,2))),1,2);
    datT=datT';
    
    clear acs_new time dat1 P T SS BB9_new
    
    %choose calibrations
    load '/Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Calibrations/cal_024.mat'
    I=find(abs(cal_date_c-mat_date_IOP(ii))<=2);
       I=find(abs(cal_date_c-mat_date_IOP(ii))<=2);
    cal_c_best=median(cal_c(I,:));
    cal_c_unc=std(cal_c(I,:));
    if ii ~= 4
        for j=1:length(prrT)
            [c_pg2(j,:),c_pg(j,:), dS(j), dT(j)]=TS_corr_ACS(datT(j,14:97)-cal_c_best, 'c', wlc, datT(j,4), datT(j,3), Tcal);
            amp0=c_pg(j,70);
            gamma0=-1;
            [x1(j,:)] = fminsearch(@least_squares_cp, [amp0 gamma0], opts, c_pg(j,:), wlc);
        end
    else
        for j=1:length(prrT)
            [c_pg2(j,:),c_pg(j,:), dS(j), dT(j)]=TS_corr_ACS(datT(j,5:88)-cal_c_best, 'c', wlc, datT(j,4), datT(j,3), Tcal);
            amp0=c_pg(j,70);
            gamma0=-1;
            [x1(j,:)] = fminsearch(@least_squares_cp, [amp0 gamma0], opts, c_pg(j,:), wlc);
        end
    end

    figure(1)

    subplot(1,3,1)
    plot(dS, -prrT)
    subplot(1,3,2)
    plot(dT, -prrT)
    subplot(1,3,3)
    plot(x1(:,2), -prrT)
    figure(2)
    plot(wlc,c_pg)

    
    pause
    outfile=['Tot_cpg_',yr,mo,day,Tot_ACS_files(ii).name(24:26),'.mat']
    save(outfile, 'wlc', 'c_pg', 'prrT', 'datT', 'cal_c_best', 'cal_c_unc');
    clear dS dT x1 c_pg datT prrT c_pg2 sum_misfit misfit cal_c_best I II cal_c_unc
    fclose('all')
    close('all')
end
