% =============================================
% M-file to process LISST-100X data "by hand" at sea instead of using LISST.exe
%
% INPUT:  profile data (*.dat), calibration factors
% OUTPUT: mat-file with fully calibrated VD [uL/L] data, binned to 2m
%       One mat-file is created for each HOT cruise
%
% M-files needed:
%       bin.m, updown.m from Chris Wingard
%       getscat.m, invert.p, and vdcorr.m from Sequoia
%
% A. White, UH, Aug 2020
% fill out parameters beginning on line 20 each cruise
% data/psd slope and zsc jpg are saved
% =============================================
    clear all
    close all
    % begin in the mfiles dir
    inversion_type= 1; %('Which inversion would you like to use?  1=sphere, 2=random, 3=n 1.17 ... ');
    mdir = (cd);  % dir for folder containing mfiles
   
%% change this section for each at sea processing
    dias = load('Dias32_b.asc');
    cd ../LISST_1421_calfiles        % change to folder containing cruise data
    fzsc_name = ([cd '/factory_zsc_1421_mar2019.asc']);
    fzsc = load('factory_zsc_1421_mar2019.asc');
    ringareafile = ([cd '\Ringarea_1421.asc']); %ring area same for 1285 and 1421
    cd ../
    mkdir('out')
    outfilename = ([cd '\out\rdata.log']);
    deepzscname = ([cd '\out\deepzsc.asc']);
    crnum = 321;
    yr = 2020;
    VCC = 31936; % from InstrumentData.txt, changes with factory calibration
    flref = 773;  % from element 36 in factory_zsc file
    class_vol = (1/6)*pi*(dias.^3); % um^3, equivalent 4/3PIr^3 and 1/6PId^3
    class_vol = class_vol*1e-9; class_vol=class_vol';    % ul, 10^9 um3 in a ul 

% get diameters 
    x = 1.25; %lower limit for first size range for type B, spherical particles
    rho=200^(1/32);
    bins(:,1) = x*rho.^([0:31]); %lower limit
    bins(:,2) = x*rho.^([1:32]); % upper limit
    bins(:,3) = sqrt(bins(:,1).*bins(:,2)); %mid-point
    bins(:,4) = bins(:,2) - bins(:,1); %width
    binwidth = bins(:,2) - bins(:,1); 
    binwidth =binwidth';
   
%% load zsc - assumes there is one zsc file
  zscfiles = dir('h*cal*.asc');
  zsc = load(zscfiles(1).name);
  figure(1)
    semilogx(dias,zsc(1:32,:))
    set(gca,'XScale','log');  grid;
    hold on
    axis tight
    ylabel('Factory Z-scat (raw counts)');
    xlabel('Particle diameter [\mum]'); 
    semilogx(dias,fzsc(1:32),'k--');
  
%% load dat files - should be only one dat file if not concatenate
  datfiles = dir('L*.DAT');
% if there were 1 file
%   [scat,tau,zsc,data_out,cscat] = getscat(datfiles(1).name,fzsc_name,1,ringareafile);
  
  % in the case of HOT321 the file split into 4 - merge them
    data_out = [];
    for i=1:length(datfiles)
        dfile = [cd '/' datfiles(i).name];
        cd mfiles
        [scat,tau,zsc,data,cscat] = getscat(dfile,fzsc_name,1,ringareafile);
        cd ../
        data_out = vertcat(data_out,data);
    end
    
    save(outfilename,'data_out','-ascii');
 
    depth= data_out(:,37)./100;
    day  = floor(data_out(:,39)/100); hour = (data_out(:,39)/100-day)*100;
    sdate = datenum(yr,1,1)+ day + hour/24;
   
    figure(2)
       plot(depth,'k.'); hold on; xlabel('sample #'); ylabel('depth')
       deepz = find(depth>= max(depth)-5); n=[1:length(depth)];
       plot(n(deepz),depth(deepz),'r.')
       deep_zsc = nanmean(data_out(deepz,:),1); %bottom 2m of each cast
       save(deepzscname,'deep_zsc','-ascii');
       figure(1)
       semilogx(dias,deep_zsc(1:32),'r-');
       legend('factory','DI','Deep')
       cd out
       fn = ['HOT-LISST_ZSC_' num2str(crnum) '.jpg'];
       print(gcf,'-djpeg','-r300',fn);
       cd ..\
       
% invert the raw data using the deep zsc
      [scat,tau,zsc,data_out,cscat] = getscat(outfilename,deepzscname,1,ringareafile); %all data corrected for deepest 2m
      [vd, dias] = invert(cscat,2,0,0,0,0,1,15); % spherical inversion % 50 iterations default
      lref = data_out(:,36); % get laser ref for VD correction
      cd(mdir)
      vdc = vdcorr(vd,VCC,flref,lref); % correct for concentration and laser power
  
% Separate up/down cast 
      figure
      plot(depth,'k.'); hold on; 
      n = input('How many casts are there?\n')
      title('Select start and end of consecutive up/down casts')
      [n1,n2]=ginput(n*2+1); n1=floor(n1); 
            
       % Separate up/down and bin
       % Bin applies 3 stdev filter to the up/down data 
       % bin_filt(data,bincol,start,stop,binsize,bintype)
        m=1;
            for j=1:2:(n*2)
              vd_down(m).depth = depth(n1(j):n1(j+1));
              vd_down(m).vd    = vdc(n1(j):n1(j+1),:);
              vd_up(m).depth = depth(n1(j+1):n1(j+2));
              vd_up(m).vd    = vdc(n1(j+1):n1(j+2),:);
              a = bin_filt([vd_down(m).depth, vd_down(m).vd],1,0,220,2,'mean');
              b = bin_filt([vd_up(m).depth, vd_up(m).vd],1,0,220,2,'mean');
              vd_down(m).vd2db = a(:,2:end); vd_up(m).vd2db = b(:,2:end);
              vd_down(m).binz = a(:,1); vd_up(m).binz = b(:,1);
              
              %  particle VD normalized to bin width = uL/L/um    
              vd_down(m).VDn = vd_down(m).vd2db./binwidth; 
              vd_up(m).VDn = vd_up(m).vd2db./binwidth; 
              % #/L        
              vd_down(m).PSD =vd_down(m).vd2db./repmat(class_vol,[size(vd_down(m).vd2db,1),1]); % #/L
              vd_up(m).PSD =vd_up(m).vd2db./repmat(class_vol,[size(vd_up(m).vd2db,1),1]); % #/L
              % #/L/um 
               vd_down(m).PSDn =  vd_down(m).PSD./binwidth;% #/L/um
               vd_up(m).PSDn =  vd_up(m).PSD./binwidth;% #/L/um
              
                figure(4)
                loglog(dias,vd_down(m).PSDn(1:75,:));xlim([ 0 100]); hold on
                loglog(dias,vd_up(m).PSDn(1:75,:));xlim([ 0 100]); hold on
                m=m+1; 
            end 
      figure(4)
       cd ..\out
       fn = ['HOT-LISST_PSD_' num2str(crnum) '.jpg'];
       xlabel('diameter'); ylabel('psd/bin width,#/L/um')
       print(gcf,'-djpeg','-r300',fn);
       cd ..\
            
%% class sums
dbz=[0:2:220];
for c=1:size(vd_down,2)
 vd_down(c).sum_sm  = nansum(vd_down(c).VDn(:,1:3),2);   % total particle volume, 1.25-2 um
 vd_down(c).sum_mid = nansum(vd_down(c).VDn(:,4:17),2);   % total particle volume, 1.25-2 um
 vd_down(c).sum_lg  = nansum(vd_down(c).VDn(:,18:27),2);   % total particle volume, 1.25-2 um
 vd_down(c).sum     = nansum(vd_down(c).VDn(:,1:27),2);   % total particle volume, 1.25-2 um

 vd_up(c).sum_sm  = nansum(vd_up(c).VDn(:,1:3),2);   % total particle volume, 1.25-2 um
 vd_up(c).sum_mid = nansum(vd_up(c).VDn(:,4:17),2);   % total particle volume, 1.25-2 um
 vd_up(c).sum_lg  = nansum(vd_up(c).VDn(:,18:27),2);   % total particle volume, 1.25-2 um
 vd_up(c).sum     = nansum(vd_up(c).VDn(:,1:27),2);   % total particle volume, 1.25-2 um
 
 vd_up(c).nL_sum   = vd_up(c).PSD(:,1:32); % #/L
 vd_down(c).nL_sum = vd_down(c).PSD(:,1:32); % #/L

 figure (5); 
    subplot(1,3,1)
        plot(vd_down(c).sum_sm,dbz); hold on; revy; xlabel('Small VDc')
     subplot(1,3,2);title('Downcast 2db VDc')
        plot(vd_down(c).sum_mid,dbz); hold on; revy;xlabel('Middle VDc')
     subplot(1,3,3)
        plot(vd_down(c).sum_lg,dbz); hold on; revy; xlabel('Large VDc')
 figure(6); 
    subplot(1,3,1)
        plot(vd_up(c).sum_sm,dbz); hold on; revy; xlabel('Small VDc')
     subplot(1,3,2);title('Upcast 2db VDc')
        plot(vd_up(c).sum_mid,dbz); hold on;revy;xlabel('Middle VDc')
     subplot(1,3,3)
        plot(vd_up(c).sum_lg,dbz); hold on;revy;xlabel('Large VDc')
    
end
            cd out
            outfile = ['HOT-LISST_PROC_' num2str(crnum) '_cast-01.mat'];
            save(outfile,'dbz', 'dias','class_vol','depth','vdc','dias','vd_up','vd_down','sdate')
            cd ..\
