

function fname = Make_ACS_Binned_from_WAP(acs_file, outdir, exclude, plotflag, saveflag)

    
    [acs_time, acs_msec, raw_c, raw_a, acs_cwl, acs_awl] = Read_acs_pctime(acs_file);

    if isempty(acs_msec)
        fname = ''
        return
    end
    
    % *** Need to deal with acs_time empty
    % If there was no valid pctime, then create acs_time as vector of
    % datenums where start of file is zero
    if isempty(acs_time)
        acs_time = datenum(0,0,0,0,0,acs_msec/1000);
    end
    
    hf1 = 0;
    awl = 10;
    cwl = 30;
    if strcmpi(plotflag,'plot')
        hf1 = figure('Name',[acs_file ' (Raw)'])
        subplot(211), hold on
            title('Raw Absorption')
            plot(acs_time,raw_a(:,awl), 'k-', 'LineWidth',0.5)
            datetick('keeplimits'), grid on
            ylabel(['Measured a(' num2str(acs_awl(awl)) ') [1/m]'])
        subplot(212), hold on
            title('Raw Attenuation')
            plot(acs_time,raw_c(:,cwl), 'k-', 'LineWidth',0.5);
            datetick('keeplimits'), grid on
            ylabel(['Measured c(' num2str(acs_cwl(cwl)) ') [1/m]'])
    end    
    
    
    % Remove filt/unfilt transitions
    tmp = datevec(acs_time);
    m = tmp(:,5);                
    raw_a(m<2|(m>=10&m<12),:) = nan;
    raw_c(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(acs_time));
        for kr = 1:size(exclude,1)
           disp(['    Masking region from ' datestr(exclude(kr,1),0) ' to ' datestr(exclude(kr,2),0)]) 
           tmp = acs_time>=exclude(kr,1) & acs_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
        raw_a(msk,:) = nan;
        raw_c(msk,:) = nan;    
    end
    

    
    % Kludgie way to deal with the filter-wheel discontinuity
    del_a = nan(size(acs_time));
    del_c = nan(size(acs_time));
    for k = 1:length(acs_time)
        if ~any(isnan(raw_a(k,30:35)))
            del_a(k) = interp1(acs_awl(30:35),raw_a(k,30:35),acs_awl(36),'linear','extrap') - raw_a(k,36);
        end
        if ~any(isnan(raw_c(k,30:35)))
            del_c(k) = interp1(acs_cwl(30:35),raw_c(k,30:35),acs_cwl(36),'linear','extrap') - raw_c(k,36);
        end
    end
    del_a = del_a - nanmedian(del_a);
    del_c = del_c - nanmedian(del_c);
    
    raw_a(find(abs(del_a)>0.002),:) = nan;
    raw_c(find(abs(del_c)>0.002),:) = nan;    
    
    if strcmpi(plotflag,'plot')
        figure(hf1)
        subplot(211), hold on
            plot(acs_time,raw_a(:,awl), 'm.')
        subplot(212), hold on
            plot(acs_time,raw_c(:,cwl), 'g.');
    end
    
    
    time0 = datevec(acs_time(1));
    time0(6) = 0;
    time1 = datevec(acs_time(end));
    time1(6) = 0;
    bin_delt = datenum(0,0,0,0,1,0)
    bin_time = (datenum(time0):bin_delt:datenum(time1))';
         
    ndata1 = nan(size(bin_time));
    ndata2a = nan(size(bin_time));
    ndata2c = nan(size(bin_time));
    ndata3a = nan(size(bin_time));
    ndata3c = nan(size(bin_time));
       
    acs_a.time = bin_time;
    acs_c.time = bin_time;
    acs_a.wavel = acs_awl;
    acs_c.wavel = acs_cwl;
    
    acs_a.med = nan(size(bin_time,1),size(raw_a,2));
    acs_a.mean = nan(size(bin_time,1),size(raw_a,2));
    acs_a.std = nan(size(bin_time,1),size(raw_a,2));
    acs_c.med = nan(size(bin_time,1),size(raw_c,2));
    acs_c.mean = nan(size(bin_time,1),size(raw_c,2));
    acs_c.std = nan(size(bin_time,1),size(raw_c,2));

    %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((acs_time>=bin_time(kt))&(acs_time<bin_time(kt)+bin_delt));
        if isempty(idx_bin)
            disp(['No data found for bin at ' datestr(bin_time(kt))])
            continue
        end
        
        ndata1(kt) = length(idx_bin);
        ndata2a(kt) = length(find(all(isfinite(raw_a(idx_bin,:)),2)));
        ndata2c(kt) = length(find(all(isfinite(raw_c(idx_bin,:)),2)));
        
        prcalo = prctile(raw_a(idx_bin,:),5);
        prcahi = prctile(raw_a(idx_bin,:),95);
        prcclo = prctile(raw_c(idx_bin,:),5);
        prcchi = prctile(raw_c(idx_bin,:),95);
        threshahi = prcahi+1.5*(prcahi-prcalo);
        threshalow = prcalo-1.5*(prcahi-prcalo);
        threshchi = prcchi+1.5*(prcchi-prcclo);
        threshclow = prcclo-1.5*(prcchi-prcclo);
        nd = length(idx_bin);
        outla = any(raw_a(idx_bin,:)<repmat(threshalow,nd,1)|raw_a(idx_bin,:)>repmat(threshahi,nd,1), 2);
        outlc = any(raw_c(idx_bin,:)<repmat(threshclow,nd,1)|raw_c(idx_bin,:)>repmat(threshchi,nd,1), 2);

        if strcmpi(plotflag,'plot')
            figure(hf1),
            subplot(211)
                plot([bin_time(kt) bin_time(kt)+bin_delt],  [threshalow(awl) threshalow(awl)], 'r--')
                plot([bin_time(kt) bin_time(kt)+bin_delt],  [threshahi(awl) threshahi(awl)], 'r--')
                plot(acs_time(idx_bin(outla)),raw_a(idx_bin(outla),awl), 'ro')
            subplot(212)
                plot([bin_time(kt) bin_time(kt)+bin_delt],  [threshclow(cwl) threshclow(cwl)], 'r--')
                plot([bin_time(kt) bin_time(kt)+bin_delt],  [threshchi(cwl) threshchi(cwl)], 'r--')
                plot(acs_time(idx_bin(outlc)),raw_c(idx_bin(outlc),cwl), 'ro')
        end
        
        raw_a(idx_bin(outla),:) = nan;
        raw_c(idx_bin(outlc),:) = nan;

        ndata3a(kt) = length(find(all(isfinite(raw_a(idx_bin,:)),2)));
        ndata3c(kt) = length(find(all(isfinite(raw_c(idx_bin,:)),2)));
        
        acs_a.mean(kt,:) = nanmean(raw_a(idx_bin,:));
        acs_a.std(kt,:) = nanstd(raw_a(idx_bin,:));
        acs_a.med(kt,:) = nanmedian(raw_a(idx_bin,:));
        acs_c.mean(kt,:) = nanmean(raw_c(idx_bin,:));
        acs_c.std(kt,:) = nanstd(raw_c(idx_bin,:));
        acs_c.med(kt,:) = nanmedian(raw_c(idx_bin,:));

    end

    disp(['100% complete'])
    toc

    
    acs_a.ns = ndata1;
    acs_a.nsclean = ndata3a;
    acs_c.ns = ndata1;
    acs_c.nsclean = ndata3c;

    
    if strcmpi(plotflag,'plot')
        figure(hf1)
        subplot(211)
            plot(bin_time+datenum(0,0,0,0,0,30), acs_a.med(:,awl), 'bo-', 'LineWidth',3.0)
            plot([bin_time+datenum(0,0,0,0,0,30) bin_time+datenum(0,0,0,0,0,30)]', ...
                [acs_a.med(:,awl)+acs_a.std(:,awl) acs_a.med(:,awl)-acs_a.std(:,awl)]', ...
                'b.-','LineWidth',2.0)
        subplot(212)
            plot(bin_time+datenum(0,0,0,0,0,30), acs_c.med(:,cwl), 'bo-', 'LineWidth',3.0)
            plot([bin_time+datenum(0,0,0,0,0,30) bin_time+datenum(0,0,0,0,0,30)]', ...
                [acs_c.med(:,cwl)+acs_c.std(:,cwl) acs_c.med(:,cwl)-acs_c.std(:,cwl)]', ...
                'b.-','LineWidth',2.0)

    end
    
    
    

    % Dissect input filename
    tmp = [strfind(acs_file,'/') strfind(acs_file,'\')];
    fname = acs_file((tmp(end)+1):end);
    fname = strrep(fname, '.','_');
    fname = strrep(fname, '21_ACS','ACS');
    fname = [outdir fname '.mat'];
     
    mkdir(outdir)
    save(fname, 'acs_a','acs_c')


return

        



function [pctime, msec,c,a,cwl,awl] = Read_acs_pctime(fname)

    disp(['    Reading ' fname ' using Read_acs_w_pctime...']);

	[dat] = textread(fname,'','delimiter','\t', 'headerlines',97);

	if ~isempty(dat)
		msec = dat(:,1);
		c = dat(:,2:84);
		a = dat(:,85:167);
	else
		[msec,a,c] = deal([])
	end
		
	cwl = [401.0	405.1	410.1	414.4	419.2	423.8	428.6	433.2	437.9	442.5	447.4	452.0	456.7	461.5	466.2	471.0	475.9	480.8	485.5	490.0	494.6	499.3	504.0	509.0	513.7	518.5	523.4	528.2	532.7	537.3	542.1	546.9	551.0	555.8	560.1	565.2	569.6	574.1	578.3	582.7	587.1	591.8	596.2	600.8	605.6	610.3	615.2	619.9	624.4	628.9	633.5	638.0	642.8	647.4	652.2	656.8	661.3	666.0	670.5	675.0	679.3	683.5	687.7	691.8	696.0	700.0	703.9	707.6	711.3	715.0	718.5	722.0	725.2	728.9	732.1	735.3	738.5	741.6	744.8	747.3	750.0	752.5	755.2];
	awl = [400.3	405.0	409.4	413.7	418.3	423.1	428.1	432.7	437.1	441.9	446.7	451.3	456.2	461.1	465.7	470.6	475.5	480.4	485.2	489.9	494.2	498.9	503.8	508.5	513.3	518.2	523.0	527.9	532.0	536.8	541.5	546.3	550.9	555.3	559.8	562.5	566.9	571.2	575.4	579.8	584.1	588.5	593.0	597.7	602.1	606.8	611.7	616.0	620.9	625.4	630.0	634.7	639.3	643.8	648.4	653.2	657.8	662.4	666.9	671.1	675.7	680.1	684.3	688.5	692.4	696.6	700.4	704.4	708.1	711.6	715.3	719.0	722.4	725.6	729.1	732.3	735.4	738.6	741.5	744.6	747.5	750.0	752.4];		
	
    
    % 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, '_ACS.');
    tmp = fname((idx-2):(idx+3));
    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


