%{
    ceilo_daily_plots_PISTON_MISOBOB.m

    incorporating parts from various DYNAMO/CALNEX scripts by LB & CF

    reads daily CL31 data files, makes plots, saves hourly and 10-min stats

    06/2018 BWB
    09/2019 EJT
%}

%% Initialize run parameters
close all;
fclose all;
clear
warning ('off','MATLAB:MKDIR:DirectoryExists');

% set month and day for start and end of run
% - because of way CL files are named...
st_month = 09;   st_day_of_month = 5;
end_month = 09;  end_day_of_month = 5;

yyyy = 2019;

% cruise ID and paths
cruise = 'PISTON_2019';    % cruise name
plot_str = 'PISTON-2019';  % cruise name, no '_' chrs
ship = 'Sally';                      % research vessel
prtit = true;                      % save plots to .png format

% system specific path defs
sysType = computer;
if strncmp(sysType,'MACI64',7)     % set Mac paths
     data_drive = '/Users/eliz/DATA/';
    path_prog = fullfile(data_drive,cruise,ship,'scientific_analysis','programs');
elseif strncmp(sysType,'PCWIN',7)  % set typical PSD DAC paths
    data_drive = 'D:\Data';
    path_prog = fullfile(data_drive,cruise,ship,'scientific_analysis','programs');
end

% matlab script path
restoredefaultpath
cd(fullfile(path_prog,'ceilo'));
addpath(genpath(fullfile(path_prog,'ceilo')));
rehash toolboxcache;

% set data and plot folders
path_raw_data = fullfile(data_drive,cruise,ship,'Ceilometer','Raw');
path_final_data = fullfile(data_drive,cruise,ship,'Ceilometer','Processed');
mkdir(fullfile(path_final_data,'cloudbase'));
mkdir(fullfile(path_final_data,'ten_min_stats'));
mkdir(fullfile(path_final_data,'hourly_stats'));

path_plots = fullfile(data_drive,cruise,ship,'Ceilometer','daily_plots');
mkdir(path_plots);
mkdir(fullfile(path_plots,'cloudbase'));
mkdir(fullfile(path_plots,'backscatter'));
mkdir(fullfile(path_plots,'layers'));
mkdir(fullfile(path_plots,'cloud_frac'));

% ET - make folder for backscatter output data
mkdir(fullfile(path_final_data,'backscatter'));

% initialize plot format
graphformat = '.png';  % select graphics format
graphdevice = '-dpng'; % select graphic device

%%  process daily files
stjd = md2yd_vn(yyyy, st_month, st_day_of_month);
endjd = md2yd_vn(yyyy, end_month, end_day_of_month);

for ddd = stjd:endjd
    close all;
    dstr = datestr(ddd + datenum(yyyy-1,12,31),1);  % for plots
    mm = datestr(ddd + datenum(yyyy-1,12,31),'mm');
    dd = datestr(ddd + datenum(yyyy-1,12,31),'dd');
    fdir = dir(fullfile(path_raw_data,['S9',mm,dd,'*.DAT'])); % EJT used to be 'A8...
    nfiles = length(fdir);
    
    % define default timestamp and empty result arrays
    % 15 sec standard resolution
    jd_ceil = ddd:15/86400:ddd+1-15/86400; jd_ceil = jd_ceil(:);
    cloudcode = NaN(size(jd_ceil)); cloud1 = cloudcode; cloud2 = cloudcode;
    cloud3 = cloudcode; ceil_range = NaN(770,1);
    ceil_b = NaN(770,length(jd_ceil));
    
    % load multiple files for ddd if required and copy into result arrays
    for ii = 1:nfiles
        fname = fdir(ii).name;
        dfl = fullfile(path_raw_data,fname);
        [jd,cldcd,cl1,cl2,cl3,range,bksktr] = load_CL31(ddd,dfl);
        if ii==1
            ceil_range = range';
        end
        for xx=1:length(jd)
            [~,zz] = min(abs(jd_ceil - jd(xx)));
            cloudcode(zz) = cldcd(xx);
            cloud1(zz) = cl1(xx);
            cloud2(zz) = cl2(xx);
            cloud3(zz) = cl3(xx);
            ceil_b(:,zz) = bksktr(:,xx);
        end
    end

    % save cloudbase data to output file
    dpath = fullfile(path_final_data,'cloudbase',[cruise,'_cloudbase_',sprintf('%03i',ddd),'.txt']);
    out = [jd_ceil,cloudcode,cloud1,cloud2,cloud3];
    save (dpath,'out','-ascii');

    % plot cloud layers detected timeseries
    il = cloudcode<=3;
    figure;
    plot(jd_ceil(il),cloudcode(il),'.');
    grid; ylim([0 3.5]); xlim([ddd,ddd+1]); yticks([0,1,2,3]);
    datetick('x','hh','keeplimits')
    xlabel([datestr(jd_ceil(1) + datenum(2019,0,0,0,0,0),'mm/dd') ' DOY ' strnd(jd_ceil(1),0)]);
%     xlabel('Decimal DOY (UTC)'); 
    ylabel('Number of Cloud Layers');
    title([plot_str,' ',dstr,', Ceilometer Cloud Base Layers Detected']);
    if prtit
        ppath = fullfile(path_plots,'layers',[cruise,'_cld_layers_',sprintf('%03i',ddd),graphformat]);
        print(graphdevice, ppath);
    end

    %%%%% save some of the ceilometer data so that I can plot it in another
    %%%%% program. 
    ceil.jd_ceil = jd_ceil;
    ceil.ceil_range = ceil_range;
    ceil.ceil_b = ceil_b;
    save_fname = [path_final_data '/backscatter/ceil_' datestr(jd_ceil(1) + datenum(2019,0,0,0,0,0),'yyyy_mmdd')];
    save(save_fname,'ceil');
    
    % plot backscatter timeseries image
    figure;
    imagesc(jd_ceil,ceil_range/1000,log10(max(0,ceil_b)));
    axis xy; ylim([0 7.5]); xlim([ddd,ddd+1]);
    datetick('x','hh','keeplimits')
    xlabel([datestr(jd_ceil(1) + datenum(2019,0,0,0,0,0),'mm/dd') ' DOY ' strnd(jd_ceil(1),0)]);
%     xlabel('Decimal DOY (UTC)'); 
    ylabel('Height (km)');
    h2 = colorbar; ylabel(h2,'1og_{10}(Backscatter), steradian');
    title([plot_str,' ',dstr,', Ceilometer Backscatter']);
    if prtit
        ppath = fullfile(path_plots,'backscatter',[cruise,'_backscatter_',sprintf('%03i',ddd),graphformat]);
        print(graphdevice, ppath);
    end

    % compute 10-min stats for ddd
    jd_bin = ddd:10/1440:ddd+1; jd_bin = jd_bin';   % 10-minute bin edges
    jd_10min = jd_bin(1:144);     % timestamp

    % empty result arrays
    Npnts = NaN(size(jd_10min)); cloud_cnt = Npnts; obscured_cnt = Npnts;
    part_obscured_cnt = Npnts; cloud_tot = Npnts; frac_cld = Npnts;
    frac_tot = Npnts; clear_cnt = Npnts; frac_clear = Npnts;
    one_layer_cnt = Npnts; multi_layer_cnt = Npnts;
    hgt_median = Npnts; hgtsm1 = Npnts; hgtsm2 = Npnts;

    for ii=1:144  % for 10 min stats
        idx = find(jd_ceil>=jd_bin(ii) & jd_ceil<jd_bin(ii+1));
        if ~all(isnan(cloudcode(idx)))
            Npnts(ii) = sum(isfinite(cloudcode(idx)));
            % number of cloud layer detects
            cloud_cnt(ii) = sum(cloudcode(idx)>0 & cloudcode(idx)<4);
            % number of totally obscured detects
            obscured_cnt(ii) = sum(cloudcode(idx)==4);
            % number of partly obscured but transparent detects
            part_obscured_cnt(ii) = sum(cloudcode(idx)==5);
            % total possible number of cloud detects
            cloud_tot(ii) = cloud_cnt(ii) + obscured_cnt(ii);
            % cloud fraction
            frac_cld(ii) = cloud_cnt(ii)/Npnts(ii); % of cloud detects
            frac_tot(ii) = cloud_tot(ii)/Npnts(ii); % of all possible
            % number of clear detects
            clear_cnt(ii) = sum(cloudcode(idx)==0);
            frac_clear(ii) = clear_cnt(ii)/Npnts(ii);
            % single / multi-layer detects
            one_layer_cnt(ii) = sum(cloudcode(idx)==1);
            multi_layer_cnt(ii) = sum(cloudcode(idx)==2 | cloudcode(idx)==3);
            if cloud_tot(ii)>0
                idx2 = (jd_ceil>=jd_bin(ii) & jd_ceil<jd_bin(ii+1) &...
                    cloudcode>0 & cloudcode<5);       % all possible cld
                hgt_median(ii) = nanmedian1(cloud1(idx2));  % median hgt, all possicle cld
                hgtsm1(ii) = prctilex(cloud1(idx2),16);     % one sigma lower limit for hgt
                hgtsm2(ii) = prctilex(cloud1(idx2),84);     % one sigma upper limit for hgt
            end
        end
    end

    % save 10-min stats to text file
    out = [jd_10min,cloud_cnt,obscured_cnt,part_obscured_cnt,cloud_tot,...
           frac_cld,frac_tot,clear_cnt,frac_clear,one_layer_cnt,multi_layer_cnt,...
           hgt_median,hgtsm1,hgtsm2];
    dfl = fullfile(path_final_data,'ten_min_stats',[cruise,'_ceilometer_10min_',sprintf('%03i',ddd),'.txt']);
    save(dfl,'out','-ascii');

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % process hourly stats
    jd_bin = ddd:1/24:ddd+1; jd_bin = jd_bin';   % hourly bin edges
    jd_hr = jd_bin(1:24);     % timestamp

    % empty result arrays
    Npnts = NaN(size(jd_hr)); cloud_cnt_hr = Npnts; obscured_cnt_hr = Npnts;
    part_obscured_cnt_hr = Npnts; cloud_tot_hr = Npnts; frac_cld_hr = Npnts;
    frac_tot_hr = Npnts; clear_cnt_hr = Npnts; frac_clear_hr = Npnts;
    one_layer_cnt_hr = Npnts; multi_layer_cnt_hr = Npnts;
    part_obs_cnt_hr = Npnts; hgt_median_hr = Npnts; hgtsm1_hr = Npnts;
    hgtsm2_hr = Npnts;

    for ii=1:24  % for hourly stats
        idx = find(jd_ceil>=jd_bin(ii) & jd_ceil<jd_bin(ii+1));
        if ~all(isnan(cloudcode(idx)))
            Npnts(ii) = sum(isfinite(cloudcode(idx)));
            % number of cloud layer detects
            cloud_cnt_hr(ii) = sum(cloudcode(idx)>0 & cloudcode(idx)<4);
            % number of totally obscured detects
            obscured_cnt_hr(ii) = sum(cloudcode(idx)==4);
            % number of partly obscured but transparent detects
            part_obscured_cnt_hr(ii) = sum(cloudcode(idx)==5);
            % total possible number of cloud detects
            cloud_tot_hr(ii) = cloud_cnt_hr(ii) + obscured_cnt_hr(ii);
            % cloud fraction
            frac_cld_hr(ii) = cloud_cnt_hr(ii)/Npnts(ii); % of cloud detects
            frac_tot_hr(ii) = cloud_tot_hr(ii)/Npnts(ii); % of all possible
            % number of clear detects
            clear_cnt_hr(ii) = sum(cloudcode(idx)==0);
            frac_clear_hr(ii) = clear_cnt_hr(ii)/Npnts(ii);
            % single / multi-layer detects
            one_layer_cnt_hr(ii) = sum(cloudcode(idx)==1);
            multi_layer_cnt_hr(ii) = sum(cloudcode(idx)==2 | cloudcode(idx)==3);
            if cloud_tot_hr(ii)>0
                idx2 = (jd_ceil>=jd_bin(ii) & jd_ceil<jd_bin(ii+1) &...
                          cloudcode>0 & cloudcode<5);           % all possible cld
                hgt_median_hr(ii) = nanmedian1(cloud1(idx2));  % median base ht, all possicle cld
                hgtsm1_hr(ii) = prctilex(cloud1(idx2),16);     % 1 sigma lower limit for ht
                hgtsm2_hr(ii) = prctilex(cloud1(idx2),84);     % 1 sigma upper limit for ht
            end
        end
    end

    % save hourly stats to text file
    out = [jd_hr,cloud_cnt_hr,obscured_cnt_hr,part_obscured_cnt_hr,cloud_tot_hr,...
           frac_cld_hr,frac_tot_hr,clear_cnt_hr,frac_clear_hr,one_layer_cnt_hr,...
           multi_layer_cnt_hr,hgt_median_hr,hgtsm1_hr,hgtsm2_hr];
    dfl = fullfile(path_final_data,'hourly_stats',[cruise,'_ceilometer_hr_',sprintf('%03i',ddd),'.txt']);
    save(dfl,'out','-ascii');

    % cloud fraction and cloudbase plots
    figure;
    plot(jd_hr,frac_cld_hr,'o'); grid; xlim([jd_hr(1),jd_hr(end)]);
    datetick('x','hh','keeplimits')
    xlabel([datestr(jd_bin(1) + datenum(2019,0,0,0,0,0),'mm/dd') ' DOY ' strnd(jd_bin(1),0)]);
%     xlabel(['DOY (UTC)']); 
    ylabel('Cloud Fraction'); ylim([0,1]);
    title([plot_str,' ',dstr,', Ceilometer Cloud Fraction']);
    if prtit
        ppath = fullfile(path_plots,'cloud_frac',[cruise,'_cld_fraction_',sprintf('%03i',ddd),graphformat]);
        print(graphdevice, ppath);
    end

    ii = isfinite(hgt_median_hr);
    jj = isfinite(hgt_median);
    figure;
    subplot(2,1,1);
    errorbar(jd_hr(ii),hgt_median_hr(ii)/1000,(hgt_median_hr(ii)-hgtsm1_hr(ii))/1000,...
            (hgtsm2_hr(ii)-hgt_median_hr(ii))/1000,'o');
    ylabel(' Hourly median, km'); grid; xlim([jd_hr(1),jd_hr(end)]); ylim([0,7]);
    datetick('x','hh','keeplimits')
%     xlabel([datestr(jd_bin(1) + datenum(2019,0,0,0,0,0),'mm/dd') ' DOY ' strnd(jd_bin(1),0)]);
    title([plot_str,' ',dstr,', Cloudbase Ht, Median \pm 84/16th percentiles']);
    subplot(2,1,2);
    errorbar(jd_10min(jj),hgt_median(jj)/1000,(hgt_median(jj)-hgtsm1(jj))/1000,...
             (hgtsm2(jj)-hgt_median(jj))/1000,'o');
    ylabel('10 min median, km'); grid; xlim([jd_hr(1),jd_hr(end)]); ylim([0,7]);
    datetick('x','hh','keeplimits')
    xlabel([datestr(jd_bin(1) + datenum(2019,0,0,0,0,0),'mm/dd') ' DOY ' strnd(jd_bin(1),0)]);
    %     xlabel(['DOY (UTC)']);
    if prtit
        ppath = fullfile(path_plots,'cloudbase',[cruise,'_cld_base_',sprintf('%03i',ddd),graphformat]);
        print(graphdevice, ppath);
    end

%      close all
end

disp('End Ceilometer Processing')
