

function fname = Make_LISST_Binned_from_WAP(lisst_file, outdir, exclude, plotflag, saveflag)

%         lisst_file = 'Mindoro_extracted\MindoroTransit1_08_042_22_ASCII.020'
%         outdir = 'Mindoro1\'
%         exclude = [];
%         plotflag = 'plot'
%         saveflag = 'nosave'
% 

    % *** DEPENDS ON DH4 CONFIG ***
    portstr = '22_ASCII'

    [lisst_time, lisst_msec, lisst_dat40] = Read_LISST_pctime(lisst_file);    
    
    % If there was no valid pctime, then create acs_time as vector of
    % datenums where start of file is zero
    if isempty(lisst_time)
        lisst_time = datenum(0,0,0,0,0,lisst_msec/1000);
    end
    
    hf1 = 0;
    tr = [3 10 20 30];
    if strcmpi(plotflag,'plot')
        hf1 = figure('Name',[lisst_file ' (Raw)'])
        subplot(211), hold on
            title('Raw Ring Scatter')
            plot(lisst_time,lisst_dat40(:,tr(1)), 'b-', 'LineWidth',0.5)
            plot(lisst_time,lisst_dat40(:,tr(2)), 'g-', 'LineWidth',0.5)
            plot(lisst_time,lisst_dat40(:,tr(3)), 'r-', 'LineWidth',0.5)
            plot(lisst_time,lisst_dat40(:,tr(4)), 'c-', 'LineWidth',0.5)
            datetick('keeplimits'), grid on
            ylabel('Counts')
        subplot(212), hold on
            title('Transmitted and Reference Power')
            plot(lisst_time,lisst_dat40(:,36), 'k-', 'LineWidth',0.5);
            plot(lisst_time,lisst_dat40(:,33), 'r-', 'LineWidth',0.5);
            datetick('keeplimits'), grid on
            ylabel('Counts')
    end    
    
    
    % Remove filt/unfilt transitions
    tmp = datevec(lisst_time);
    m = tmp(:,5);                
    lisst_dat40(m<2|(m>=10&m<12),:) = nan;
    lisst_dat40(m<2|(m>=10&m<12),:) = nan;     
    
    
    % exclude is Nx2 matrix of datenums, each row represents a beginning
    % and end of range to exclude
    if ~isempty(exclude)
        msk = repmat(logical(0),size(lisst_time));
        for kr = 1:size(exclude,1)
           disp(['    Masking region from ' datestr(exclude(kr,1),0) ' to ' datestr(exclude(kr,2),0)]) 
           tmp = lisst_time>=exclude(kr,1) & lisst_time<=exclude(kr,2);
           if isempty(find(tmp))
               disp('    ** No data found for this region')
           else
               msk = msk | tmp;  % Build mask from rows of exclude matrix
           end       
        end
        % NaN out masked data
        lisst_dat40(msk,:) = nan;
    end
    

    
    time0 = datevec(lisst_time(1));
    time0(6) = 0;
    time1 = datevec(lisst_time(end));
    time1(6) = 0;
    bin_delt = datenum(0,0,0,0,1,0)
    bin_time = (datenum(time0):bin_delt:datenum(time1))';
             
    dat40_bin.time = bin_time;
    dat40_bin.ns = nan(size(bin_time,1),1);
    dat40_bin.med = nan(size(bin_time,1),40);
    dat40_bin.mean = nan(size(bin_time,1),40);
    dat40_bin.std = nan(size(bin_time,1),40);
    dat40_bin.p5 = nan(size(bin_time,1),40);
    dat40_bin.p95 = nan(size(bin_time,1),40);
    dat40_bin.p16 = nan(size(bin_time,1),40);
    dat40_bin.p84 = nan(size(bin_time,1),40);
        
    %keyboard
    tic
    for kt = 1:size(bin_time,1)

        if mod(kt,fix(size(bin_time,1)/10))==0
            disp([num2str(kt/fix(size(bin_time,1)/10)*10) '% complete (' num2str(toc) ' sec)'])
        end

        idx_bin = find((lisst_time>=bin_time(kt))&(lisst_time<bin_time(kt)+bin_delt));
        if isempty(idx_bin)
            disp(['No data found for bin at ' datestr(bin_time(kt))])
            continue
        end
        
        dat40_bin.ns(kt) = length(idx_bin);
        dat40_bin.mean(kt,:) = nanmean(lisst_dat40(idx_bin,:));
        dat40_bin.std(kt,:) = nanstd(lisst_dat40(idx_bin,:));
        dat40_bin.med(kt,:) = nanmedian(lisst_dat40(idx_bin,:));
        dat40_bin.p5(kt,:) = prctile(lisst_dat40(idx_bin,:),5);
        dat40_bin.p95(kt,:) = prctile(lisst_dat40(idx_bin,:),95);
        dat40_bin.p16(kt,:) = prctile(lisst_dat40(idx_bin,:),16);
        dat40_bin.p84(kt,:) = prctile(lisst_dat40(idx_bin,:),84);
        
    end

    disp(['100% complete'])
    toc


    lisst_scat.time = dat40_bin.time;
    lisst_scat.ns = dat40_bin.ns;
    lisst_scat.mean = dat40_bin.mean(:,1:32);
    lisst_scat.std = dat40_bin.std(:,1:32);
    lisst_scat.med = dat40_bin.med(:,1:32);
    lisst_scat.p5 = dat40_bin.p5(:,1:32);
    lisst_scat.p95 = dat40_bin.p95(:,1:32);
    lisst_scat.p16 = dat40_bin.p16(:,1:32);
    lisst_scat.p84 = dat40_bin.p84(:,1:32);
        
    lisst_tr.time = dat40_bin.time;
    lisst_tr.ns = dat40_bin.ns;
    lisst_tr.mean = dat40_bin.mean(:,33);
    lisst_tr.std = dat40_bin.std(:,33);
    lisst_tr.med = dat40_bin.med(:,33);
    lisst_tr.p5 = dat40_bin.p5(:,33);
    lisst_tr.p95 = dat40_bin.p95(:,33);
    lisst_tr.p16 = dat40_bin.p16(:,33);
    lisst_tr.p84 = dat40_bin.p84(:,33);
        
    lisst_ref.time = dat40_bin.time;
    lisst_ref.ns = dat40_bin.ns;
    lisst_ref.mean = dat40_bin.mean(:,36);
    lisst_ref.std = dat40_bin.std(:,36);
    lisst_ref.med = dat40_bin.med(:,36);
    lisst_ref.p5 = dat40_bin.p5(:,36);
    lisst_ref.p95 = dat40_bin.p95(:,36);
    lisst_ref.p16 = dat40_bin.p16(:,36);
    lisst_ref.p84 = dat40_bin.p84(:,36);
        
    
    
    
    
    
    
    
%             hf1 = figure('Name',[lisst_file ' (Raw)'])
%         subplot(211), hold on
%             title('Raw Ring Scatter')
%             plot(lisst_time,lisst_dat40(:,tr), '-', 'LineWidth',0.5)
%             datetick('keeplimits'), grid on
%             ylabel('Counts')
%         subplot(212), hold on
%             title('Transmitted and Reference Power')
%             plot(lisst_time,lisst_dat40(:,36), 'k-', 'LineWidth',0.5);
%             plot(lisst_time,lisst_dat40(:,33), 'r-', 'LineWidth',0.5);
%             datetick('keeplimits'), grid on
%             ylabel('Counts')

    if strcmpi(plotflag,'plot')
        hf2 = figure('Name',[lisst_file ' (Binned)'])
        subplot(211), hold on
        title('Raw Ring Scatter')        
        plot(bin_time+datenum(0,0,0,0,0,30), lisst_scat.med(:,tr(1)), 'bo-', 'LineWidth',3.0)
        plot([bin_time+datenum(0,0,0,0,0,30) bin_time+datenum(0,0,0,0,0,30)]', ...
            [lisst_scat.p5(:,tr(1)) lisst_scat.p95(:,tr(1))]', ...
            'b.-','LineWidth',2.0)
        plot(bin_time+datenum(0,0,0,0,0,30), lisst_scat.med(:,tr(2)), 'go-', 'LineWidth',3.0)
        plot([bin_time+datenum(0,0,0,0,0,30) bin_time+datenum(0,0,0,0,0,30)]', ...
            [lisst_scat.p5(:,tr(2)) lisst_scat.p95(:,tr(2))]', ...
            'g.-','LineWidth',2.0)
        plot(bin_time+datenum(0,0,0,0,0,30), lisst_scat.med(:,tr(3)), 'ro-', 'LineWidth',3.0)
        plot([bin_time+datenum(0,0,0,0,0,30) bin_time+datenum(0,0,0,0,0,30)]', ...
            [lisst_scat.p5(:,tr(3)) lisst_scat.p95(:,tr(3))]', ...
            'r.-','LineWidth',2.0)
        plot(bin_time+datenum(0,0,0,0,0,30), lisst_scat.med(:,tr(4)), 'co-', 'LineWidth',3.0)
        plot([bin_time+datenum(0,0,0,0,0,30) bin_time+datenum(0,0,0,0,0,30)]', ...
            [lisst_scat.p5(:,tr(4)) lisst_scat.p95(:,tr(4))]', ...
            'c.-','LineWidth',2.0)
        
        set(gca,'XLim',[min(lisst_time) max(lisst_time)])
        datetick('keeplimits'), grid on
        ylabel('Counts')
        
        subplot(212), hold on
            title('Transmitted and Reference Power')

        plot(bin_time+datenum(0,0,0,0,0,30), lisst_ref.med, 'ko-', 'LineWidth',3.0)
        plot([bin_time+datenum(0,0,0,0,0,30) bin_time+datenum(0,0,0,0,0,30)]', ...
            [lisst_ref.p5 lisst_ref.p95]', ...
            'k.-','LineWidth',2.0)
            
        plot(bin_time+datenum(0,0,0,0,0,30), lisst_tr.med, 'ro-', 'LineWidth',3.0)
        plot([bin_time+datenum(0,0,0,0,0,30) bin_time+datenum(0,0,0,0,0,30)]', ...
            [lisst_tr.p5 lisst_tr.p95]', ...
            'r.-','LineWidth',2.0)

        
            
    end
    
    
    

    % Dissect input filename
    tmp = [strfind(lisst_file,'/') strfind(lisst_file,'\')];
    fname = lisst_file((tmp(end)+1):end);
    fname = strrep(fname, '.','_');
    fname = strrep(fname, portstr,'LISST');
    fname = [outdir fname '.mat'];
    
    if ~exist(outdir, 'dir'), mkdir(outdir), end;
    save(fname, 'lisst_scat', 'lisst_tr', 'lisst_ref')


return




function [pctime, msec, dat40] = Read_LISST_pctime(fname)

    disp(['  Reading ' fname ' using rd_lisst...']);
    [msec, dat] = textread(fname,'%u\t%s', 'delimiter','\n', 'whitespace',' ');
        
    % find all the starts of records {
    dat = strtrim(dat);
    packetidx = find(strcmp('{', dat))+1;
    if isempty(packetidx)
       pctime = [];
       msec = [];
       dat40 = [];
       return        
    end
    % deal with the possibility that the last packetidx doesnt have a whole
    % 40 rows of data
    if length(dat)-packetidx(end)<40
        packetidx(end) = [];
    end
    
    % predeclare size of data matrix
    dat40 = repmat(-888,length(packetidx),40);
    
    for kp = 1:length(packetidx)
        datrow = str2num(char(dat(packetidx(kp):(packetidx(kp)+39))))';
        if isempty(datrow) | size(datrow,2)~=40
            datrow = repmat(-888,1,40);
        end
        dat40(kp,:) = datrow;
    end

    msec = msec(packetidx);

    % Read PC clock ASCII data
    % PC time file looks like ...
    % 9000	08/20/06	23:00:10.200
    % 10000	08/20/06	23:00:11.200
    % 11000	08/20/06	23:00:12.200
    % 12000	08/20/06	23:00:13.200
    
    % Look for the extracted pc time file at wap id 19
    idx = findstr(fname, '_ASCII.');
    tmp = fname((idx-2):(idx+5));
    fname2 = strrep(fname, tmp, '19_ASCII');
    disp(['    Looking for PC time file ' fname2])
    
    if exist(fname2)
        disp('      Found. Reading.')
        
        [clockmsec, MM,DD,YY,hh,mm,sssss] = textread(fname2,'%u %d/%d/%d %d:%d:%f');
        clocktime = datenum(2000+YY,MM,DD,hh,mm,sssss); % Assumes year is >2000
        disp(['      PC times from '  datestr(clocktime(1)) ' to '  datestr(clocktime(end))])
                
        % Special case where clockmsec has non-distinct msec, interp1
        % doesn't handle this case. I remove these elements and let interp1
        % interpolate over these regions, so timing might be off there
        idx = find(diff(clockmsec)==0);
        clockmsec(idx) = [];
        clocktime(idx) = [];

        pctime = interp1(clockmsec,clocktime, msec, 'linear','extrap');
       
    else
        warning('PC time file NOT FOUND')

        disp('    Not found, returning empty pctime')
        pctime = [];
        
    end
    
    
return


        
