%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This program is designed to process all the sensors of the slowdrop
%package but the AC-S

%% Setup the Import Options
clear all
close all

opts = spreadsheetImportOptions("NumVariables", 9);

% Specify sheet and range
opts.Sheet = "Sheet1";
opts.DataRange = "A2:I48";

% Specify column names and types
opts.VariableNames = ["year", "mon", "day", "hr", "min", "run", "filt", "Epoch", "dark"];
opts.SelectedVariableNames = ["year", "mon", "day", "hr", "min", "run", "filt", "Epoch", "dark"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "string", "double", "double", "double"];
opts = setvaropts(opts, 6, "WhitespaceRule", "preserve");
opts = setvaropts(opts, 6, "EmptyFieldRule", "auto");

% Import the data
log_data = readtable("/Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/log_IOP_pack.xlsx", opts, "UseExcel", false);
clear opts

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

%% Clear temporary variables

for ii=[1:2,4:7,9:11,13,15:17,19:37,42,44,46] %3 & 8 12 14 18 43 45 47have no CTD files 38 39 40 41-no time stamp on CTD
    if log_data.Epoch(ii) == 1
        direct='/Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Epoch1/';
    elseif log_data.Epoch(ii) == 2
        direct='/Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Epoch2/';
    else
        direct='/Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Epoch3/';
    end
    cast_time=datenum(log_data.year(ii),log_data.mon(ii),log_data.day(ii),log_data.hr(ii),log_data.min(ii),0)
    ext=char(log_data.run(ii)),
    
    %read total profiles
    if ii<42 == (1 || 2)
        CTD_file=dir([direct,'run_24_CTD-ENGR.',ext]);
        bb9_file=dir([direct,'run_22_ECO.',ext]);
        CDOM_file=dir([direct,'run_23_ECO.',ext]);
        Anal_file=dir([direct,'run_29_ANALOG.',ext]);
    else
        CTD_file=dir([direct,'archive_24_CTD-ENGR.',ext]);
        bb9_file=dir([direct,'archive_22_ECO.',ext]);
        CDOM_file=dir([direct,'archive_23_ECO.',ext]);
        Anal_file=dir([direct,'archive_29_ANALOG.',ext]);
    end
    
    [time, P, T, S, C]=rd_SB37([direct,CTD_file.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));
    %read chl. fluorometer
    [time_F F]=rd_analog([direct,Anal_file.name]);
    
    %read CDOM Fluorometer
    [time_FDOM FDOM]=rd_Eco1ch([direct,CDOM_file.name]);
    
    %load BB9 data:
    if isempty(bb9_file)==0
        [time_BB9 wl BB9]=rd_BB9counts([direct,bb9_file.name]);
        time_bb9=time;
        
        %median filter BB9 with a 3-pt median filter (3s)
        
        for j=1:9
            BB9_new(:,j)=medfilt1(BB9(:,j),3);
        end
        clear BB9
        
        %determine begining and end time for merging
        I=find(time_ctd>=max([min(time_ctd); min(time_FDOM);min(time_F);min(time_BB9)]));
        I1=find(time_ctd<=min([max(time_ctd); max(time_FDOM);max(time_F);max(time_BB9)]));
        I2=find(time_ctd>=time_ctd(min(I)) & time_ctd<=time_ctd(max(I1)));
        time_merge=time_ctd(I2);
        dat1(:,1)=(time_merge-min(time_merge))*24*3600; %time in secs'
        dat1(:,[2:1:4])=[P(I2),T(I2),SS(I2)];
        dat1(:,5)=interp1(time_F,F,time_merge-12/24/3600,'linear');
        dat1(:,6)=interp1(time_FDOM,FDOM,time_merge,'linear');
        for j=1:9
            dat1(:,6+j)=interp1(time_BB9,BB9_new(:,j),time_merge+4/24/3600,'linear');
        end
        
        figure(1)
        plot((dat1(:,3)-min(dat1(:,3)))./(max(dat1(:,3))-min(dat1(:,3))),-dat1(:,2),'b')
        hold on
        plot((dat1(:,9)-min(dat1(:,9)))./(max(dat1(:,9))-min(dat1(:,9))),-dat1(:,2),'k')
        plot((dat1(:,5)-min(dat1(:,5)))./(max(dat1(:,5))-min(dat1(:,5))),-dat1(:,2),'g')
        plot((dat1(:,6)-min(dat1(:,6)))./(max(dat1(:,6))-min(dat1(:,6))),-dat1(:,2),'r')
        hold off
        pause
        
        [datT,prrT]=binpre(dat1',dat1(:,2),0.5,floor(max(dat1(:,2))),2,1);
%        [datT,prrT]=binpre(dat1',dat1(:,2),ceil(min(dat1(:,2))),floor(max(dat1(:,2))),1,1);
        datT=datT';
        
        wl_WL=[412,440,488,510,532,595,650,715,880];
        slopes_WL=[1.29e-5 1.415e-5 1.247e-5 1.008e-5 8.666e-6 5.333e-6 4.054e-6 3.24e-6 2.525e-6];
        darks_fac=[49 51 51 59 56 52 55 56 53];
        darks_field=[48.79	50.94	50.76	61.66	54.26	51.39	55.38	56.65 53];
        slopes_SS=[1.45e-5 1.45e-5 1.27e-5 9.74e-6 8.72e-6 5.11e-6 3.95e-6 3.16e-6 2.55e-6];
        wl_SS=[407,439,485,507,527,594,651,715,878];
        %apply calibrations
        for i=1:length(prrT)
            [betasw_WL(i,:),beta90sw,bsw]= betasw_ZHH2009(wl_WL,datT(i,3),124,datT(i,4));
            [betasw_SS(i,:),beta90sw,bsw]= betasw_ZHH2009(wl_SS,datT(i,3),124,datT(i,4));
        end
        datT(:,16:24)=((datT(:,7:15)-darks_field).*slopes_WL-betasw_WL)*2*pi*1.1;
        datT(:,7:15)=((datT(:,7:15)-darks_field).*slopes_SS-betasw_SS)*2*pi*1.1;
        figure(1)
        plot(wl_SS,datT(:,7:15))
        hold on
        figure(2)
        plot(wl_WL,datT(:,16:24))
        %minimization routine
        amp0=mean(datT(:,13));
        gamma0=1;
        wl=[407,439,485,527,594,651,715,878];
        for j=1:length(datT(:,13))
                 bbp=datT(j,[7:9,11:15]);
                [x1 Fval] = fminsearch(@least_squares_cp, [amp0 gamma0], opts, bbp, wl);
                gamma(j)=x1(2);
                bbp650(j)=x1(1);
        end
        figure(1)
        plot(wl,bbp650.*(wl'/650).^-gamma)
        hold off
        pause
        
        fileout=['Epoch',num2str(log_data.Epoch(ii)),'_',ext,'.mat'];
        
        save(fileout, 'datT', 'cast_time', 'gamma', 'prrT')
        clear gamma bbp650 BB9 FDOM F datT BB9_new time_merge dat1 betasw_WL betasw_SS
    else
        %determine begining and end time for merging
        I=find(time_ctd>=max([min(time_ctd); min(time_FDOM);min(time_F)]));
        I1=find(time_ctd<=min([max(time_ctd); max(time_FDOM);max(time_F)]));
        I2=find(time_ctd>=time_ctd(min(I)) & time_ctd<=time_ctd(max(I1)));
        time_merge=time_ctd(I2);
        dat1(:,1)=(time_merge-min(time_merge))*24*3600; %time in secs'
        dat1(:,[2:1:4])=[P(I2),T(I2),SS(I2)];
        dat1(:,5)=interp1(time_F,F,time_merge-12/24/3600,'linear');
        dat1(:,6)=interp1(time_FDOM,FDOM,time_merge,'linear');
        [datT,prrT]=binpre(dat1',dat1(:,2),0.5,floor(max(dat1(:,2))),2,1);
        datT=datT';
        
        fileout=['Epoch',num2str(log_data.Epoch(ii)),'_',ext,'.mat'];
        save(fileout, 'datT', 'cast_time', 'prrT')
        clear gamma bbp650 BB9 FDOM F datT time_merge dat1
    end
end
