%{
    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
    11/2019 EJT
%}

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

%% cruise ID and paths
cruise = 'PISTON_2019';                 % Cruise name
ship = 'Sally_Ride';                    % Research vessel

cruise_title = strrep(cruise,'_',' ');  % Cruise name with space instead of underscore for plots
yyyy = cruise(end-3:end);               % year string
year = str2double(yyyy);                % year number

% set month and day for start and end of run
%%% PISTON 2019 9/5 - 9/25
st_month = 9;   st_day_of_month = 5;
end_month = 9;  end_day_of_month = 25;

%% plot and path things
prtit = true;    % save plots

% system specific path defs
machine=char(java.net.InetAddress.getLocalHost.getHostName); % doesn't work with floating unnamed IP nodenames!
if regexp(machine,'squall') % Simon's paths
        data_drive = '/Users/sdeszoek/Data/cruises/ATOMIC_2020/RHB/PSD';
        path_prog = fullfile(data_drive,cruise,ship,'scientific_analysis','programs');
else % Elizabeth paths
    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
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,'codes'));
mkdir(fullfile(path_plots,'cloud_frac'));

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

% 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
    
%% read data
    close all;
    dstr = datestr(ddd + datenum(yyyy,0,0),1);         % for plots
    fstr = datestr(ddd + datenum(yyyy,0,0),'mmdd');    % for files
    mm = datestr(ddd + datenum(yyyy,0,0),'mm');
    dd = datestr(ddd + datenum(yyyy,0,0),'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 = NaN(size(jd_ceil)); cloud2 = NaN(size(jd_ceil));
    cloud3 = NaN(size(jd_ceil)); 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

%     1    jd_ceil       decimal day of year at start of interval
%     2    cloudcode     cloud detection code, as follows:
%                             0 = clear sky
%                             1 = one cloud layer detected
%                             2 = two cloud layers detected
%                             3 = three cloud layers detected
%                             4 = totally obscured
%                             5 = partly obscured but determined transparent
%     3    cloud1        base height of lowest layer, m
%     4    cloud2        base height of next highest layer, m
%     5    cloud3        base height of highest layer, m


%% process 10-min stats
    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; cloud_cnt_any = Npnts;
    part_obscured_cnt = Npnts; cloud_tot = Npnts; frac_cld = Npnts; 
    frac_tot = Npnts; clear_cnt = Npnts; frac_clear = Npnts; frac_any = Npnts; 
    one_layer_cnt = Npnts; multi_layer_cnt = Npnts;
    hgt_50 = Npnts; hgt_16 = Npnts; hgt_84 = Npnts;
    hgt_5 = Npnts;    hgt_10 = Npnts;    hgt_15 = Npnts;     hgt_18 = 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); % cloud 1-3 present per time step
            frac_tot(ii) = cloud_tot(ii)/Npnts(ii); % cloud 1-4 present per time step
            % 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 1,2,3,4,5
                hgt_50(ii) = nanmedian1(cloud1(idx2));  % median hgt, all possicle cld
                hgt_16(ii) = prctilex(cloud1(idx2),16);     % one sigma lower limit for hgt ~ 16%
                hgt_84(ii) = prctilex(cloud1(idx2),84);     % one sigma upper limit for hgt ~ 84%
                hgt_5(ii) = prctilex(cloud1(idx2),5);     % 
                hgt_10(ii) = prctilex(cloud1(idx2),10);     % 
                hgt_15(ii) = prctilex(cloud1(idx2),15);     % 
                hgt_18(ii) = prctilex(cloud1(idx2),18);     % 
            end
        end
    end


%% 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_50_hr = Npnts; hgt_16_hr = Npnts;
    hgt_84_hr = Npnts; hgt_5_hr = Npnts;    hgt_10_hr = Npnts;    
    hgt_15_hr = Npnts;     hgt_18_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 1,2,3,4,5
                hgt_50_hr(ii)   = nanmedian1(cloud1(idx2));  % 50% = median
                hgt_16_hr(ii)   = prctilex(cloud1(idx2),16);     % 1 sigma lower limit for ht ~16%
                hgt_84_hr(ii)   = prctilex(cloud1(idx2),84);     % 1 sigma upper limit for ht ~84%
                hgt_5_hr(ii)    = prctilex(cloud1(idx2),5);     
                hgt_10_hr(ii)   = prctilex(cloud1(idx2),10);      
                hgt_15_hr(ii)   = prctilex(cloud1(idx2),15);      
                hgt_18_hr(ii)   = prctilex(cloud1(idx2),18);      
            end
        end
    end

%% plots 

    plot_it = 0;
    if plot_it == 1

        % plot cloud codes detected timeseries
    il = cloudcode<=5; % so... all codes
    figure;
    plot(jd_ceil(il),cloudcode(il),'.');
    grid; ylim([0 5.5]); xlim([ddd,ddd+1]); yticks([0,1,2,3,4,5]);
    datetick('x','hh','keeplimits')
    xlabel([datestr(jd_ceil(1) + datenum(yyyy,0,0,0,0,0),'mm/dd') ' DOY ' strnd(jd_ceil(1),0)]);
%     xlabel('Decimal DOY (UTC)'); 
    ylabel('Cloud Code');
    text(jd_ceil(1), 0.2, 'clear','fontweight','bold','fontsize',14);
    text(jd_ceil(1), 1.2, '1 layer','fontweight','bold','fontsize',14);
    text(jd_ceil(1), 2.2, '2 layers','fontweight','bold','fontsize',14);
    text(jd_ceil(1), 3.2, '3 layers','fontweight','bold','fontsize',14);
    text(jd_ceil(1), 4.2, 'obscured','fontweight','bold','fontsize',14);
    text(jd_ceil(1), 5.2, 'transparent partly obscured','fontweight','bold','fontsize',14);
    title([cruise_title,' ',dstr,', Ceilometer Cloud Codes']);
    if prtit
        ppath = fullfile(path_plots,'codes',[cruise '_cld_codes_' fstr '_DOY_' sprintf('%03i',ddd) graphformat]);
        print(graphdevice, ppath);
    end

    
    % plot backscatter timeseries image to 7.5 km and 3 km
    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(yyyy,0,0,0,0,0),'mm/dd') ' DOY ' strnd(jd_ceil(1),0)]);
    ylabel('Height (km)');
    h2 = colorbar; ylabel(h2,'1og_{10}(backscatter), steradian^{-1} km^{-1}');
    title([cruise_title,' ',dstr,', Ceilometer backscatter']);
    grid on;
    set(gca,'layer','top');
    if prtit
        ppath = fullfile(path_plots,'backscatter',[cruise '_backscatter_' fstr '_DOY_' sprintf('%03i',ddd) graphformat]);
        print(graphdevice, ppath);
    end
        %%% same but only up to 3 km height
    ylim([0 3]);
    grid on;
    set(gca,'layer','top');
    if prtit
        ppath = fullfile(path_plots,'backscatter',[cruise '_backscatter_3km_' fstr '_DOY_' sprintf('%03i',ddd) graphformat]);
        print(graphdevice, ppath);
    end    
    
    % cloud fraction and cloudbase plots
    figure; hold on;
    
    plot(jd_hr+(1/48),frac_cld_hr,'-om'); 
    plot(jd_10min,frac_tot,'.'); grid; 
    xlim([jd_hr(1),jd_hr(end)]);
    legend('1-hour','10-min','location','best');
    datetick('x','hh','keeplimits')
    xlabel([datestr(jd_bin(1) + datenum(yyyy,0,0,0,0,0),'mm/dd') ' DOY ' strnd(jd_bin(1),0)]);
%     xlabel(['DOY (UTC)']); 
    ylabel('Cloud Fraction'); ylim([0,1]);
    title([cruise_title,' ',dstr,', Ceilometer Cloud Fraction']);
    if prtit
        ppath = fullfile(path_plots,'cloud_frac',[cruise '_cld_fraction_' fstr '_DOY_' sprintf('%03i',ddd) graphformat]);
        print(graphdevice, ppath);
    end

    ii = isfinite(hgt_50_hr);
    jj = isfinite(hgt_50);
    figure;
    subplot(2,1,1);
    errorbar(jd_hr(ii),hgt_50_hr(ii)/1000,(hgt_50_hr(ii)-hgt_16_hr(ii))/1000,...
            (hgt_84_hr(ii)-hgt_50_hr(ii))/1000,'o');
    ylabel('km'); grid; xlim([jd_hr(1),jd_hr(end)]); ylim([0,7]);
    datetick('x','hh','keeplimits')
    text(jd_hr(2), 0.95*max(ylim),'1-hour','fontweight','bold','fontsize',20);
%     xlabel([datestr(jd_bin(1) + datenum(yyyy,0,0,0,0,0),'mm/dd') ' DOY ' strnd(jd_bin(1),0)]);
    title([cruise_title,' ',dstr,', cloud base heights: 16%, 50%, 84% values of 1st cloud layer']);
    subplot(2,1,2);
    errorbar(jd_10min(jj),hgt_50(jj)/1000,(hgt_50(jj)-hgt_16(jj))/1000,...
             (hgt_84(jj)-hgt_50(jj))/1000,'o');
    ylabel('km'); grid; xlim([jd_hr(1),jd_hr(end)]); ylim([0,7]);
    datetick('x','hh','keeplimits')
    text(jd_hr(2), 0.95*max(ylim),'10-minute','fontweight','bold','fontsize',20);
    xlabel([datestr(jd_bin(1) + datenum(yyyy,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_' fstr '_DOY_' sprintf('%03i',ddd) graphformat]);
        print(graphdevice, ppath);
    end

    
    % plot combo image
    figure;
    imagesc(jd_ceil,ceil_range/1000,log10(max(0,ceil_b)));
    hold on;
    b = plot(jd_hr+(1/48), hgt_5_hr/1000,'om'); % plot centered on half hours
    a = plot(jd_10min, hgt_5/1000,'--w'); % 
%     b = plot(jd_hr, hgtsm10_hr/1000,'--w'); % 
%     c = plot(jd_hr, hgtsm15_hr/1000,'--w'); % 
%     d = plot(jd_hr, hgtsm18_hr/1000,'-w'); % 
%     a = plot(jd_hr, hgtsm1_hr/1000,'-w'); % 16th percentile
%     b = plot(jd_hr, hgt_median_hr/1000,'-w'); % 50th percentile
%     c = plot(jd_hr, hgtsm2_hr/1000,'-k'); % 84th percentile
    axis xy; ylim([0 3]); xlim([ddd,ddd+1]);
    datetick('x','hh','keeplimits')
    xlabel([datestr(jd_ceil(1) + datenum(yyyy,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^-1 km^-1');
%     leg = legend([a, b, c, d], '5 %','10 %', '15 %','18 %','location','north');
    leg = legend([a, b], '5% 10-min','5% 1-hour','location','north');
    title(leg, 'cloud base height');
    title([cruise_title,' ',dstr,', Ceilometer backscatter']);
    grid on;
    set(gca,'layer','top');
    if prtit
        ppath = fullfile(path_plots,'backscatter',[cruise '_backscatter_combo_3km_' fstr '_DOY_' sprintf('%03i',ddd) graphformat]);
        print(graphdevice, ppath);
    end
    end  %%% if plot_it == 1
    
%      close all

%% legacy text files

    % 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_50_hr,hgt_16_hr,hgt_84_hr];
    dfl = fullfile(path_final_data,'hourly_stats',[cruise '_ceilometer_hr_' fstr '_DOY_' sprintf('%03i',ddd) '.txt']);
    save(dfl,'out','-ascii');

    % 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_50,hgt_16,hgt_84];
    dfl = fullfile(path_final_data,'ten_min_stats',[cruise '_ceilometer_10min_' fstr '_DOY_' sprintf('%03i',ddd) '.txt']);
    save(dfl,'out','-ascii');

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

%% stats to save for matlab statistics files
    clear c;

    % year day
    c.jd = jd_ceil;
    c.jd_10min = jd_10min;
    c.jd_hr = jd_hr;
    % matlab datetime
    c.t           = jd_ceil + datenum(2020,0,0,0,0,0);
    c.t_10min     = jd_10min + datenum(2020,0,0,0,0,0);
    c.t_hr        = jd_hr + datenum(2020,0,0,0,0,0);
    % time since 2020 1/1 00:00
    start_day =  datenum(2020,1,1,0,0,0);
    c.te          = (c.t - start_day) *24*3600;
    c.te_10min    = (c.t_10min - start_day) *24*3600;
    c.te_hr       = (c.t_hr - start_day) *24*3600;

    % variables to save
    c.range = ceil_range;
    c.ceil_b = ceil_b;
    c.cloud_code = cloudcode; % cloud detection code; note 5 is never detected
    c.cloud_base_first = cloud1;  % base height of first cloud layer detected
    c.cloud_base_second = cloud2; % base height of second cloud layer detected
    c.cloud_base_third = cloud3; % base height of third cloud layer detected

    % cloud fraction in terms of time
    c.frac_cloud      = frac_tot;       % percentage of time period when classification of 1,2,3,4,5 happened
    c.frac_cloud_hr   = frac_tot_hr;    % percentage of time period when classification of 1,2,3,4,5 happened
    c.frac_clear      = frac_clear;     % percentage of time period when classification of 0 happened
    c.frac_clear_hr   = frac_clear_hr;  % percentage of time period when classification of 0 happened

    % counts of code detections per time step.... ? is it just counting
    % detections or adding the numbers?
    c.count_cloud_total                       = cloud_tot;        % any of 1,2,3,4,5
    c.count_cloud_total_hr                    = cloud_tot_hr;     % any of 1,2,3,4,5
    c.count_cloud_layers_1_to_3               = cloud_cnt;        % any of 1,2,3
    c.count_cloud_layers_1_to_3_hr            = cloud_cnt_hr;     % any of 1,2,3
    c.count_cloud_obscured                    = obscured_cnt;     % any of 4
    c.count_cloud_obscured_hr                 = obscured_cnt_hr;  % any of 4
    c.count_cloud_part_obsured_transparent    = part_obscured_cnt;    % any of 5 %%% this never gets classified?
    c.count_cloud_part_obsured_transparent_hr = part_obscured_cnt_hr; % any of 5 %%% this never gets classified? 
    c.count_clear                             = clear_cnt;        % 0
    c.count_clear_hr                          = clear_cnt_hr;     % 0

    % cloud base height = 5th pecentile of lowest layer detected, first_cloud_layer, cloud1 in 15-sec file
    c.cloud_base        = hgt_5/1000;   % km 
    c.cloud_base_hr     = hgt_5_hr/1000; % km
    
    % percentile levels of the height of the clowest cloud layer detected,
    % cloud1 in 15 sec file, for cloud types 1-5 (all types). The different
    % percentile levels tell you something about how variable the clouds
    % are, how sheared they are, etc. We choose a "best" percentile to be
    % our estimate of cloud base height. In ATOMIC 2020, it was 5th
    % percentile as tested with plots in this program. Investigate to see
    % what threshold helps keep the line from jumping up and down a lot. 
    c.prcnt_cld_5         = hgt_5/1000; % km
    c.prcnt_cld_5_hr      = hgt_5_hr/1000; % km
    c.prcnt_cld_10        = hgt_10/1000; % km
    c.prcnt_cld_10_hr     = hgt_10_hr/1000; % km    
    c.prcnt_cld_16        = hgt_16/1000; % km
    c.prcnt_cld_16_hr     = hgt_16_hr/1000; % km
    c.prcnt_cld_50        = hgt_50/1000; % km
    c.prcnt_cld_50_hr     = hgt_50_hr/1000; % km
    c.prcnt_cld_84        = hgt_84/1000; % km
    c.prcnt_cld_84_hr     = hgt_84_hr/1000; % km

    
    save_fname = [path_final_data '/all/' cruise '_ceilometer_' datestr(c.t_hr(1),'mmdd') '_DOY_' sprintf('%03i',ddd)];
    save(save_fname,'c');
    disp(['saving file: ' save_fname]);

end % for each day

disp('End Ceilometer Processing')