%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This program is designed to merge BB9 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

%% Read data on all files available
opts = spreadsheetImportOptions("NumVariables", 8);

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

% Specify column names and types
opts.VariableNames = ["run", "date", "start_time", "end_time", "north_lat", "south_lat", "east_long", "west_long"];
opts.VariableTypes = ["string", "double", "double", "double", "double", "double", "double", "double"];

% Specify variable properties
opts = setvaropts(opts, "run", "WhitespaceRule", "preserve");
opts = setvaropts(opts, "run", "EmptyFieldRule", "auto");

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

wl_WL=[412,440,488,510,532,595,650,715,880];
wl_SS=[407,439,485,507,527,594,651,715,878];
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_fac=[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];
slopes_fau=[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=wl_SS; 
slopes=slopes_fau;
darks=darks_field;

 %Global parameters
opts = optimset('fminsearch');      
opts = optimset(opts,'MaxIter',2000000); 
opts = optimset(opts,'MaxFunEvals',200000);
opts = optimset(opts,'TolX',1e-8);
opts = optimset(opts,'TolFun',1e-8);


for i=22:25%1:length(BB{:,1})
    if i<=10
        direct='/Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Epoch1/';
    elseif i>10 & i<=21
        direct='/Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Epoch2/';
    else
        direct='/Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/Epoch3/';
    end
    
    date_str(i,:)=[num2str(BB{i,2}),'_'];
    if i>1 & date_str(i,:)==date_str(i-1,:)
        date_str(i,:)=[num2str(BB{i,2}),'b'];
    end
    if i>2 & date_str(i,1:8)==date_str(i-1,1:8) & date_str(i,1:8)==date_str(i-2,1:8)
        date_str(i,:)=[num2str(BB{i,2}),'c'];
    end
    %output file
    textfile=['EXPORTS-EXPORTSNP_ECO-BB9_process_',date_str(i,:),'_bbp_R1.sb'];
    docfile='Processing_ECOBB9_slowdrop_package.docx';
    
   for iwl=1:length(wl)
       wlarray{iwl}=['bbp',num2str(wl(iwl)),','];
       wlarray2{iwl}=['bbp',num2str(wl(iwl)),'_sd,'];
   end
   apunits=repmat('1/m,',1,2*length(wl));
   start_date=num2str(BB{i,2});
   end_date=num2str(BB{i,2});
   start_time=[num2str(floor(BB{i,3}*24)),':',num2str((BB{i,3}*24-floor(BB{i,3}*24))*60)];
   end_time=[num2str(floor(BB{i,4}*24)),':',num2str((BB{i,4}*24-floor(BB{i,4}*24))*60)];
   latN=[num2str(BB{i,5})];
   latS=[num2str(BB{i,6})];
   lonE=[num2str(BB{i,7})];
   lonW=[num2str(BB{i,8})];

    
    fnam=char(BB{i,1})

    %reading data in and getting their sampling time
    [time, P, T, S, C]=rd_SB37([direct,'run_24_CTD-ENGR.',fnam]);
    if i>=27 
        time=[0:1:length(time)-1]*0.0000157409995154012;
    end
    
    figure(1)
    plot(time,-P)
    [X ,Y]=ginput(2)

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

    [time wl BB9]=rd_BB9counts([direct,'run_22_ECO.',fnam]);
    time_bb9=time;
    if i>=27 
        time_bb9=[0:1:length(time)-1]*0.0000115748712640448;
    end

    %determine begining and end time for merging.
    I=find(time_ctd>=time_bb9(1) & time_ctd<=max(time_bb9));
    time_merge=time_ctd(I);

    %create interpolated data matrix
    dat1(:,1)=(time_merge-min(time_merge))*24*3600; %time in secs'
    dat1(:,[2:4])=[P(I),T(I),SS(I)];
    for j=1:9
        dat1(:,4+j)=interp1(time_bb9,BB9(:,j),time_merge,'linear');
    end
    %bin data to 2m bins
    depth_array=[0.5:2:110.5];
    for j=1:55
       I=find(dat1(:,2)>depth_array(j) & dat1(:,2)<=depth_array(j+1));
       if (isfinite(I))
          if length(I)>1
              dat2(j,1:13)=median(dat1(I,:));
              dat2(j,14:22)=std(dat1(I,5:13));
              N(j)=length(I);
          else 
              dat2(j,1:13)=dat1(I,:);
              N(j)=1;
              dat2(j,14:22)=[NaN NaN NaN NaN NaN NaN NaN NaN NaN];
          end
       else
          dat2(j,1:13)=ones(1,13)*NaN;
          dat2(j,14:22)=[NaN NaN NaN NaN NaN NaN NaN NaN NaN];
          N(j)=0;
       end
    end
    dat2(:,2)=[1.5:2:109.5];
    

    %apply calibrations
   for j=1:55
        [betasw(1,:),beta90sw,bsw]= betasw_ZHH2009(wl,dat2(j,3),124,dat2(j,4));
        dat2(j,5:13)=((dat2(j,5:13)-darks).*slopes-betasw)*2*pi*1.1;
        dat2(j,14:22)=dat2(j,14:22).*slopes*2*pi*1.1;
        amp0=dat2(j,11);
        gamma0=1;
        bbp=dat2(j,[5:7,9:12]);
        unc=dat2(j,[14:16,18:21]);
        if isfinite(bbp(2))==1
            [x1 Fval] = fminsearch(@least_squares_bbp, [amp0 gamma0], opts, bbp, unc, wl([1:3,5:8]));
            gamma(j)=x1(2);
            bbp650(j)=x1(1);
            bbp470=interp1(wl,dat2(j,5:13),470,'linear');
            C_phyto(j) = bbp470*12138-0.59;
        else
           gamma(j)=NaN;
            bbp650(j)=NaN;
            C_phyto(j) = NaN;
        end
   end
   COV=dat2(:,14:22)./dat2(:,5:13);
   for j=1:55
       I=find(COV(j,:)>0.5);
       dat2(j,5+I-1)=NaN;
       dat2(j,14+I-1)=NaN;
   end

   figure(2)
   subplot(2,2,1)
   plot(dat2(:,3),-dat2(:,2))
   ylabel('depth')
   xlabel('Temp')
   subplot(2,2,2)
   plot(dat2(:,4),-dat2(:,2))
   ylabel('depth')
   xlabel('Sal')
   subplot(2,2,3)
   plot(bbp650,-dat2(:,2))
   ylabel('depth')
   xlabel('bbp')
   subplot(2,2,4)
   plot(gamma,-dat2(:,2),'*')
   hold on
   plot(dat2(:,12)./dat2(:,13),-dat2(:,2),'-k')
   hold off
   ylabel('depth')
   xlabel('\gamma_{bbp}')
   export_fig(['bb9_fig',num2str(i),'a.png'])
   
   figure(3)
   subplot(2,1,1)
   plot(wl,dat2(:,5:13))
   ylabel('bbp(\lambda)')
   xlabel('wavelength')
   subplot(2,1,2)
   plot(wl,dat2(:,14:22)./dat2(:,5:13))
   ylabel('bbp_{COV}')
   xlabel('wavelength')
   export_fig(['bb9_fig',num2str(i),'b.png'])
   
   fileID=fopen(textfile,'w');

   %header:
   fprintf(fileID,'%-30s\n','/begin_header');
   fprintf(fileID,'%-30s\n','/investigators=Emmanuel_Boss,Nils_Haentjens');
   fprintf(fileID,'%-30s\n','/affiliations=UMaine');
   fprintf(fileID,'%-30s\n','/contact=emmanuel.boss@maine.edu');
   fprintf(fileID,'%-30s\n','/experiment=EXPORTS');
   fprintf(fileID,'%-30s\n','/cruise=EXPORTSNP');
   fprintf(fileID,'%-30s\n','/station=NA');
   fprintf(fileID,'%-30s\n',['/data_file_name=',textfile]);
   fprintf(fileID,'%-30s\n',['/documents=',docfile]);
   fprintf(fileID,'%-100s\n','/calibration_files=BB9-132_CharSheets_BOSS.pdf,CERT_BB9_132_181119actual.pdf,');
   fprintf(fileID,'%-30s\n','/data_type=cast');
   % fprintf(fileID,'%-30s\n','/data_status=final');
   fprintf(fileID,'%-30s\n',['/start_date=',start_date]);
   fprintf(fileID,'%-30s\n',['/end_date=',end_date]);
   fprintf(fileID,'%-30s\n',['/start_time=',start_time,':00[GMT]']);
   fprintf(fileID,'%-30s\n',['/end_time=',end_time,':00[GMT]']);
   fprintf(fileID,'%-30s\n',['/north_latitude=',latN,'[DEG]']);
   fprintf(fileID,'%-30s\n',['/south_latitude=',latS,'[DEG]']);
   fprintf(fileID,'%-30s\n',['/east_longitude=',lonE,'[DEG]']);
   fprintf(fileID,'%-30s\n',['/west_longitude=',lonW,'[DEG]']);
   fprintf(fileID,'%-30s\n','/water_depth=-999');
   fprintf(fileID,'%-30s\n','/delimiter=comma');
   fprintf(fileID,'%-30s\n','/missing=-999');
   fprintf(fileID,'%-30s\n','/instrument_model=ECOBB9');
   fprintf(fileID,'%-50s\n','/instrument_manufacturer=WETLabs');
   
   %spectral output
   fprintf(fileID,'%-30s\n','!COMMENTS');
   fprintf(fileID,'%-30s\n','!bbp values are based on binning 2m measurements');
   fprintf(fileID,'%-240s\n',['/fields=pressure,Wt,sal,',wlarray{:},wlarray2{:},'bincount,Phyto_Carbon,bbp_gamma']);
   fprintf(fileID,'%-150s\n',['/units=dbar,degreesC,PSU,',apunits,'none,ug/L,unitless']);
   fprintf(fileID,'%-15s\n','/end_header');
   
   
   %define the format for the data part formdate=repmat('%9s ',1,3);
   formrad=repmat('%8.7f,',1,(length(wl)*2));
   format=['%6.4f,','%6.4f,','%6.4f,',formrad,'%3.0f,','%6.4f,','%6.4f\n'];
      %delete lines with no data
   AA=dat2(:,2:22);
   B=AA*NaN;
   index=[];
   k=1;
   for l=1:55
       if isfinite(AA(l,6))
           index(k)=l;
           III=find(isnan(AA(l,:)));
           AA(l,III)=-999;
           k=k+1;
       end
   end
   
   
   III=isfinite(C_phyto)==0;
   C_phyto(III)=-999;
   III=isfinite(gamma)==0;
   gamma(III)=-999;
   III=isfinite(N)==0;
   N(III)=-999;

      
   data1=[AA(index,:),N(index)',C_phyto(index)',gamma(index)'];
   
   for ii = 1:length(index)
       fprintf(fileID,format,data1(ii,:));
   end
 
   fclose(fileID)
   clear dat2 dat1 gamma C_phyto N
end

 





