clear all;
close all;
set_machine_et;

% Code to intake ODM level 1 data from ODM_save_all.m and calculate DSD parameters
% based on code from Kyla Drushka, Christian Klepp, Merhala Thurai, and adapted by EJT
% for SPURS-2 and PISTON deployments in summer 2020, updated Oct 2020
% CANTMAT takes drop number density N(D) mm^-1 m^-3 (see cant.txt files)
% gamma fitting code takes drop count data?
% ... (( drop number density N(D) mm^-1 m^-3 (see dsd.dat files) ))???

%%% NOTE: G-J Huang at CSU, creator of the original code, does not 
%    recommend interpolating the DSD to a regular size bin interval 
%    because it would introduce errors. 


% cruise = 'PISTON';
cruise = 'SPURS2';
version_num = 'v1';

path_data = ['/Users/eliz/Dropbox/DSD/' cruise];

indir = [path_data '/final/']; 
outdir = [path_data '/final/'];
plotdir = [path_data '/analysis/'];

% load latest DSD data file (might be updating if on the cruise currently)
these_files = dir([indir cruise '_dsd_level_1_*.mat']);
files = char(these_files(:).name);
[nfiles, nchars] = size(files);
load([indir files(end,:)]);
disp(['latest DSD file: ' files(end,:)]);

%==========================================================================
%   prepare data for processing
%==========================================================================

%%% info from DSD
nt = DSD.ntimes;        % number of times
nbins = DSD.nbins;      % number of bins
d = DSD.D;              % diameters of bins
Nd = DSD.Nd;            % drops per size bin per m3

%%% simple drop count per size bin
ndrops = DSD.nrain; 
ndrops_total = sum(sum(DSD.nrain));     % for fun... total number of all drops in all bins in all times
wh_no_drops = find(DSD.nsumrain == 0);  % no drops in this time step

%%% wind speed sensed by cup anemometer attached to disdrometer
% urel_anem = DSD.urel_anem;

start_dsd_processing = 1;
if start_dsd_processing == 1

%==========================================================================
%    gamma fitting
%==========================================================================

% Check each 1 min dsd record to make sure there are at least
% 4 consecutive size bins with non-zero values. If this condition is met, an entry
% at the end of the size bin dimension will have a "1" in it to
% indicate that it's a good row/column. This does NOT rewrite the dsd
% array. It just adds a dummy or ID variable at the end (k =  nbins+1) as a check.
n_check = zeros(1, nbins); % Nd(:,1:end);
for j=1:nt
  for k=1:nbins-3
    if(ndrops(j,k)>0 && ndrops(j,k+1)>0 && ndrops(j,k+2)>0 && ndrops(j,k+3)>0)
      n_check(j)=1;
      break;
    end
  end
end

n_check(j) = 0;

% if the 4 consecutive bin condition was not met, don't calculate the rest
% of the parameters
% In prior version, this nan'ed the row if condition wasn't met: ndrops(ii,1:nbins) = NaN;


%  Create dummy output arrays for gamma  distribution fitting
Dm_avg      = nan(1,nt);
W_avg       = nan(1,nt);   
N_avg       = nan(1,nt);
omega_0     = nan(1,nt);   
omega_1     = nan(1,nt);   
omega_2     = nan(1,nt);
omega_3     = nan(1,nt);   
omega_4     = nan(1,nt);   
omega_5     = nan(1,nt);
omega_6     = nan(1,nt);
D_mean      = nan(1,nt);   
ED2_mean    = nan(1,nt);  
var_mean    = nan(1,nt);
mu_sigma    = nan(1,nt);  
sigma       = nan(1,nt);     
mean_val    = nan(1,nt);
sigmamN     = zeros(1,nt);    
sigmamD     = zeros(1,nt);

% Fill other arrays
for i=1:nt  
  if n_check(i)==0
    continue;
  end
  test=trapz(d,Nd(i,1:nbins).*d.^3);
  if test<=0.001
    continue;
  end
  Dm_avg(i)=trapz(d,Nd(i,1:nbins).*d.^4)/test;
     
  W_avg(i)=test*pi*.001/6;
     
  N_avg(i)=(W_avg(i)/(Dm_avg(i)^4))*81.487E03;
  omega_0(i)=trapz(d,Nd(i,1:nbins));
  omega_1(i)=trapz(d,Nd(i,1:nbins).*d);
  omega_2(i)=trapz(d,Nd(i,1:nbins).*d.^2);
  omega_3(i)=trapz(d,Nd(i,1:nbins).*d.^3);
  omega_4(i)=trapz(d,Nd(i,1:nbins).*d.^4);
  omega_5(i)=trapz(d,Nd(i,1:nbins).*d.^5);
  omega_6(i)=trapz(d,Nd(i,1:nbins).*d.^6);
  sigmamN(i)=trapz(d,Nd(i,1:nbins).*d.^3.*(d-Dm_avg(i)).^2);
  sigmamD(i)=(Dm_avg(i)^2)*omega_3(i);
   
% sigma is the "shape" of the dsd according to the ratio sigma/Dm
% or normalized std dev of the mass spectrum wrt to Dm
% the estimate of mu can be obtained as sigma/Dm=1/(4+mu)^.5
 
  sigma(i)=(sigmamN(i)/sigmamD(i))^0.5;
  if sigma(i)>0
    mu_sigma(i)=(1./(sigma(i)^2))-4;
  else
    mu_sigma(i)=NaN;
  end
 
  D_mean(i)=omega_1(i)/omega_0(i);
  ED2_mean(i)=omega_2(i)/omega_0(i);
  mean_val(i)=omega_2(i)/omega_0(i);

  var_mean(i)=ED2_mean(i)-D_mean(i)^2;
  if(var_mean(i)< .01)
    var_mean(i)=NaN;
  end

% Calculating Omega_n based on eq(21) of Sempere-Torres et al paper
% Note A_n in that equation is assumed to be unity.
end

% Fitting Omega_n=a(W**b) : See eq(22) of ST paper. Note W=psi there
%  choose W_avg threshold at 0.05
% First separate out steady rain from unsteady rain!!
% Take out this part

%%% test the code by outputting several values to the screen
% display('mean,std  of Dm_avg for R>=5 mm/h in convective rain');
% nanmean(Dm_avg)
% nanstd(Dm_avg)
% display('log (<Nw>) ');
% log10(nanmean(N_avg))
% display('mean and std of log(Nw)');
% nanmean(log10(N_avg))
% nanstd(log10(N_avg))
% display('mean and std of sigma');
% nanmean(sigma)
% nanstd(sigma)


% Normalize each dsd by N(D)/Nw and D/Dm
W_avg_tr=W_avg';
j_test=0;
d_norm = nan(nt,nbins);   Nd_norm = d_norm ;
for i = 1:nt
  if n_check(i)==0
    continue;
  end
  j_test=j_test+1;
  d_norm(i,1:nbins)=d(1:nbins)./Dm_avg(i);
  Nd_norm(i,1:nbins)=Nd(i,1:nbins)./N_avg(i);
end
display(['number of original DSD samples: ' strnd(length(Nd),0)]);
display(['number of points in normalized dsd in steady rain: ' num2str(j_test)]);


%xlabel('D/Dm');
%ylabel('N(D)/Nw');

mu=-3:0.5:15;
nummu=size(mu);
mu_init = 3;

% get new arrays ready
muvec_n         = nan(nt, 1);
Dofinal         = nan(nt, 1);
rainrate_avg_g  = nan(nt, 1);
f_converge      = nan(nt, 1);
Mu_init         = nan(nt,1);
d_norm_g        = nan(nt,1);
Nd_g            = nan(nt,1);
rep_r_g         = nan(nt,1);
temp_r_g        = nan(nt,1);
%  Using fminsearch to find the mu
for i = 1:nt
  if n_check(i)==0
    muvec_n(i) = NaN;
    Dofinal(i) = NaN;
    rainrate_avg_g(i) = NaN;
    f_converge(i) = NaN;
    Mu_init(i) = NaN;
    continue;
  end  
  D_norm = squeeze(d_norm(i,:));   N_norm = squeeze(Nd_norm(i,:));
  itr = 0;  flag = 0;   Mu = mu_init;
  while flag == 0 && itr<=3000
    itr = itr+600;
    options = optimset('MaxIter',itr);
    [my_mu,err,flag]=fminsearch('mu_cost',Mu,options,D_norm,N_norm);
    Mu = my_mu;
  end
  if flag == 1 && my_mu>=-3 && my_mu<=15
    muvec_n(i)=my_mu;
    f_converge(i) = 1;                   % Converge
    Mu_init(i) = mu_init;
  else
    itr = 3000;   err_now = err;
    if my_mu<=-3
      mu_now = -3;
    elseif my_mu>=15
      mu_now = 15;
    else
      mu_now = my_mu;
    end

    options = optimset('MaxIter',itr);
    if my_mu < 3 
      for m = 1:4
        Mu = mu_init+3*m;
        [my_mu,err,flag]=fminsearch('mu_cost',Mu,options,D_norm,N_norm);
        if flag == 1
          break;
        end
      end
    else
      for m = 1:2
        Mu = mu_init-3*m;
        [my_mu,err,flag]=fminsearch('mu_cost',Mu,options,D_norm,N_norm);
        if flag == 1
          break;
        end
      end
    end
    if flag == 1
      Mu_init(i) = Mu;
      if my_mu>=-3 && my_mu<=15
        muvec_n(i) = my_mu;
        f_converge(i) = 1;                   % Converge
      elseif my_mu<-3
        muvec_n(i) = -3;
        f_converge(i) = 2;                   % converge but force to mu limit
      else
        muvec_n(i) = 15;
        f_converge(i) = 2; 
      end
    else
      Mu_init(i) = NaN;
      f_converge(i) = 0;                     % not converge
      disp(['i = ' num2str(i) ' not converge']);
      if err>=err_now
          muvec_n(i)=mu_now;
      else
        if my_mu<=-3
          muvec_n(i) = -3;
        elseif my_mu>15
          muvec_n(i) = 15;
        else
          muvec_n(i) = my_mu;
        end
      end
    end
  end
  Dofinal(i)=Dm_avg(i)*(3.67+muvec_n(i))/(4.+muvec_n(i));

% Calculate the rainrate for my normalized gamma fitted dsd N(D)/Nw=f(mu)f(D/Dm)
% have to de-normalize before calculating R. 
%set upper limit of integration to be either 2.5*Dm or 8 mm which ever
%is less.
% un-normalize the same way as normalization ie re-multiply by Nw
% see variable temps(i) below 
  upper_lim = min([2.5*Dm_avg(i) 8]);

  dia = 0.35:.1:upper_lim;
  nbin2  = length(dia);
%  dia_act(i,1:nbin2)=dia(1:nbin2);
  d_norm_g(i,1:nbin2) = dia(1:nbin2)./Dm_avg(i);
  k_gam = 6/(4^4);
  f_mu = k_gam*((4+muvec_n(i)).^(4+muvec_n(i)));
  temps = N_avg(i)*(f_mu)/(gamma(muvec_n(i)+4));
  Nd_g(i,1:nbin2)=temps*(d_norm_g(i,1:nbin2).^(muvec_n(i))).*exp(-(4+muvec_n(i))*d_norm_g(i,1:nbin2));
  temp_r_g(i,1:nbin2)=(Nd_g(i,1:nbin2).*dia(1:nbin2).^3).*(9.65-10.3*exp(-.6.*dia(1:nbin2)));
  rainrate_avg_g(i)=trapz(dia(1:nbin2),temp_r_g(i,1:nbin2))*6*pi*1E-4;
end 

%Plot all expt normalized dsds according to N/N_w vs f(D/Dm) 
% and N/nc vs f(D/Dm)  Also plot the 
% various normalized gamma functions for all mu -values

% nanstd(muvec_n)

% finally, calculate the last variables
Nw_D0 = nan(nt,1);
for i = 1:nt
   ff = (3.67*(muvec_n(i)+4)/(4*(3.67+muvec_n(i))))^4;
   Nw_D0(i) = N_avg(i)*ff;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%   SAVE THE DATA    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear DSDN;

%%% how many times did we get an Nw estimate?
wh_good = find(isfinite(Nw_D0) == 1);

% add on to DSD structure that we started with... except use the above
% wh_good ID to get everything on the same time base. rename the new
% structure as DSDN.

ff = fields(DSD);
for j = 1:length(ff)
    thisvar = DSD.(ff{j});
    [ti, tj] = size(thisvar);
    % for fields discretized in time...
    if ti > 100
        DSDN.(ff{j}) = thisvar(wh_good,:);
    % for all other metadata fields...
    else
        DSDN.(ff{j}) = thisvar;
    end
end

ID_gamma_good = zeros(1,nt);
ID_gamma_good(wh_good) = 1;

% add new variables to DSDN

DSDN.Nw         = Nw_D0(wh_good);     % normalized gamma distribution intercept parameter, mm-1 m-3
DSDN.D0         = Dofinal(wh_good);   % median volume diameter, mm
DSDN.mu         = muvec_n(wh_good);   % shape parameter, unitless
DSDN.ID_good_gamma = ID_gamma_good;   % ID of good retrievals in same length as original DSD time series
DSDN.good_gamma = wh_good;            % indicies of good retrievals from original DSD time series
DSDN.ntimes     = length(DSDN.t);     % redone after thresholding

DSDN.README_GAMMA = {'This structure adds to the original DSD structure: level 1 processing '; ...
                   'DSD were processed with code similar to Thompson et al. 2015 (JAS).'; ...
                   'The code estimates parameters of a normalized gamma distribution ';...
                   'for each raining time step that contains at least 4 consecutive bins of drops.'; ...
                   'The noramlized gamma DSD model is a fairly good assumption for short time steps.'; ...
                   'The three gamma DSD parameters are estimated from the samples (observations)'; ...
                   'using a cost function minimization optimation process. This determines which';...
                   'three parameters best represent the sample population, i.e. ';...
                   'they could theoretically recreate it.'; ...
                   'Nw is normalized gamma DSD intercept parameter in per mm per m^-3 units'; ...
                   'Nw is a proxy for the number concentration.'; ...
                   'D0 is the median drop diameter in mm'; ...
                   'mu is the size parameter, which describes the breadth and shape of the distribution'; ...
                   'good_gamma = indices of level 1 DSD time series where gamma parameter retrieval was possible '; ...
                   'ID_good_gamma = 1 or 0 for original DSD time series points when gamma retrieval was possible '; ...
                  };

%%% use similar file name, but indicate this is level 2... with derived DSD
%%% parameters and additional QC thresholding. 
first_time = DSDN.t(1);
last_time = DSDN.t(end);

time_stamp = [datestr(first_time,'YYYY_mmdd') '_to_' datestr(last_time,'YYYY_mmdd')];

filename = [outdir cruise '_dsd_level_2_' version_num '.mat'];    
save(filename, 'DSDN');

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


display(['number of Nw estimates: ' strnd(length(wh_good),0)]);

plot_simple_NwD0 = 1;

if plot_simple_NwD0 == 1

figure('position',[1,1,500, 500]); hold on;
plot(DSDN.D0, log10(DSDN.Nw),'ok','markersize',7);
ylim([2.5 5]);
xlim([0 2.5]);
plot(xlim,[3.85, 3.85],'--y','linewidth',3);
ylabel('log_{10}N_{W}');
xlabel('D_{0} [mm]');
box off;
%title('Test:triple the number count');
set(gca,'fontsize',19,'fontweight','bold','xtick',(0:0.5:2.5));
%title('ODM original dataset: 129 irregular size bins from 0.0124 - 22.2838 mm');
%legend(['Good Nw estimates: n = ' strnd(length(wh_good),0) ' out of ' strnd(nt,0) ' samples'],'location','best');
title('cruise');
annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{cruise},'FontSize',14,'FitBoxToText','off','LineStyle','none');
export_fig([plotdir cruise '_odm_output_' version_num ],'-png','-nocrop');

end

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

plot_NwD0 = 1;
if plot_NwD0 == 1
    
figure('Position',[1,1,800, 800]); hold on;
plot(DSDN.D0, log10(DSDN.Nw),'o');
xlabel('D_0 [mm]');
ylabel('log_{10}N_{W}');
xlim([0 3]);
plot(xlim,[3.85, 3.85],'--','color',rgb('orange'),'linewidth',3);
box off;
grid on;
title(cruise);
legend(['Good Nw estimates: n = ' strnd(length(wh_good),0) ' out of ' strnd(nt,0) ' samples'],...
    'convective vs. stratiform threshold = 3.85 (Thompson et al. 2015)','location','best');
export_fig([plotdir cruise '_test_odm_output_original_' version_num],'-png','-nocrop');

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

figure('Position',[400,1,800, 800]);
scatter(DSDN.D0, log10(DSDN.Nw),100,DSDN.r,'o','filled');
hold on;
% wh_r10 = find(DSDN.rgauge > 10);
% plot(DSDN.D0(wh_r10), log10(DSDN.Nw(wh_r10)),'o','markeredgecolor','k','markersize',9);
c = colorbar;
cbpath = '/Users/eliz/Documents/MATLAB/colortables/';
load([cbpath 'nice_gfdl.rgb']);
colormap(nice_gfdl(10:end, :));
ylabel(c,'R [mm h^{-1}');
xlabel('D_0 [mm]');
ylabel('log_{10}N_{W}');
xlim([0 3]);
plot(xlim,[3.85, 3.85],'--','color',rgb('skyblue'),'linewidth',3);
box off;
grid on;
title(cruise);
legend(['Good Nw estimates: n = ' strnd(length(wh_good),0) ' out of ' strnd(nt,0) ' samples'],...
    'convective vs. stratiform threshold = 3.85 (Thompson et al. 2015)','location','best');
export_fig([plotdir cruise '_test_odm_output_R_' version_num],'-png','-nocrop');

%%%

% figure('Position',[400,1,800, 800]);
% scatter(DSDN.D0, log10(DSDN.Nw),100,DSDN.rgauge,'o','filled');
% hold on;
% wh_r10 = find(DSDN.rgauge > 10);
% % plot(DSDN.D0(wh_r10), log10(DSDN.Nw(wh_r10)),'o','markeredgecolor','k','markersize',9);
% c = colorbar;
% cbpath = '/Users/eliz/Documents/MATLAB/colortables/';
% load([cbpath 'nice_gfdl.rgb']);
% colormap(nice_gfdl(10:end, :));
% ylabel(c,'R_gauge [mm h^{-1}');
% xlabel('D_0 [mm]');
% ylabel('log_{10}N_{W}');
% xlim([0 3]);
% plot(xlim,[3.85, 3.85],'--','color',rgb('skyblue'),'linewidth',3);
% box off;
% grid on;
% title('ODM original dataset: 129 irregular size bins from 0.0124 - 22.2838 mm');
% legend(['Good Nw estimates: n = ' strnd(length(wh_good),0) ' out of ' strnd(nt,0) ' samples'],...
%     'convective vs. stratiform threshold = 3.85 (Thompson et al. 2015)','location','best');
% export_fig([plotdir cruise '_test_odm_output_Rgauge_' version_num],'-png','-nocrop');

%%%

figure('Position',[400,1,800, 800]);
scatter(DSDN.D0, log10(DSDN.Nw),100,DSDN.refl,'o','filled');
hold on;
% wh_r10 = find(DSDN.rgauge > 10);
% plot(DSDN.D0(wh_r10), log10(DSDN.Nw(wh_r10)),'o','markeredgecolor','k','markersize',9);
c = colorbar;
cbpath = '/Users/eliz/Documents/MATLAB/colortables/';
load([cbpath 'nice_gfdl.rgb']);
colormap(nice_gfdl(10:end, :));
ylabel(c,'Z_h dBZ');
xlabel('D_0 [mm]');
ylabel('log_{10}N_{W}');
xlim([0 3]);
plot(xlim,[3.85, 3.85],'--','color',rgb('skyblue'),'linewidth',3);
box off;
grid on;
title(cruise);
legend(['Good Nw estimates: n = ' strnd(length(wh_good),0) ' out of ' strnd(nt,0) ' samples'],...
    'convective vs. stratiform threshold = 3.85 (Thompson et al. 2015)','location','best');
export_fig([plotdir cruise '_test_odm_output_dbz_' version_num],'-png','-nocrop');

%%%

figure('Position',[400,1,800, 800]);
scatter(DSDN.D0, log10(DSDN.Nw),100,DSDN.lwc,'o','filled');
hold on;
% wh_r10 = find(DSDN.rgauge > 10);
% plot(DSDN.D0(wh_r10), log10(DSDN.Nw(wh_r10)),'o','markeredgecolor','k','markersize',9);
c = colorbar;
cbpath = '/Users/eliz/Documents/MATLAB/colortables/';
load([cbpath 'nice_gfdl.rgb']);
colormap(nice_gfdl(10:end, :));
ylabel(c,'LWC g m^{-3}');
xlabel('D_0 [mm]');
ylabel('log_{10}N_{W}');
xlim([0 3]);
plot(xlim,[3.85, 3.85],'--','color',rgb('skyblue'),'linewidth',3);
box off;
grid on;
title(cruise);
legend(['Good Nw estimates: n = ' strnd(length(wh_good),0) ' out of ' strnd(nt,0) ' samples'],...
    'convective vs. stratiform threshold = 3.85 (Thompson et al. 2015)','location','best');
export_fig([plotdir cruise '_test_odm_output_lwc_' version_num],'-png','-nocrop');


end

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

%%% problems right now: (first time through)



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

%%% problem now: (second time through with N(D)

%%% log10(Nw) is low in convection, but not by an order of magnitude

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

end  %%% loop to do Nw D0 processing

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% psd_conc = odm(:,22:149); %%% Number Concentration per bin
%%%% ndrops = odm_num(:,22:149);   %%% Drop count per bin
%%%% d = 128 size bins

nbins = length(d);

N_mode = nan(1,nbins);
N_mean = nan(1,nbins);
N_med = nan(1,nbins);
N_5 = nan(1,nbins);
N_25 = nan(1,nbins);
N_75 = nan(1,nbins);
N_95 = nan(1,nbins);
N_CI_w = nan(1,nbins);
N_count = nan(1,nbins);

%%% stats of drop counts
for i = 1:nbins
   N_mode(i) = mode(ndrops(:,i));
   N_mean(i) = nanmean(ndrops(:,i));
   N_med(i)  = nanmedian(ndrops(:,i));
   N_5(i) = prctile(ndrops(:,i),5);
   N_25(i) = prctile(ndrops(:,i),25);
   N_75(i) = prctile(ndrops(:,i),75);
   N_95(i) = prctile(ndrops(:,i),95);
   N_CI_w(i) = conf_int(ndrops(:,i));
   N_count(i) = sum(ndrops(:,i));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% plot 2D histogram of N(D)

plot_N_of_D = 1;
if plot_N_of_D == 1
    
figure('Position',[100,200,500,500]);
size_bins_ex = repmat(d(1:70),nt,1);
ndrops_ex = ndrops(:,1:70);

N = ndhist(size_bins_ex, ndrops_ex); 
set(gcf,'Colormap',flipud(pink));
c = colorbar; ca = caxis;
set(get(c,'ylabel'),'string','Freq. [%]','HorizontalAlignment','Right');
set(c,'Location','North');
ylabel('N(D)');
xlabel('D [mm]');
set(gca,'yscale','log');
%ylim([0,1.5]);
xlim([0.25,4]);
title(cruise)
%hold on; plot([0 100], [0.25 0.25], '--k','lineWidth',2);
%title('(a)');
%set(gca, 'xtick',dbzticks,'Xticklabel',dbzticklabels,'ytick',zdrticks,'yticklabel',zdrticklabels);
set(gca,'TickDir','out','TickLength',[0.02 0.03]);
%set(gca,'FontSize',13);
%%%%%%%%%
export_fig([plotdir cruise '_N_of_D_2Dhist_' version_num],'-png','-nocrop');

end

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

%%% Try calculating mean N(D) for myself... 

% N_mean = nan(1,nbins);
% N_med = nan(1,nbins);
% N_5 = nan(1,nbins);
% N_25 = nan(1,nbins);
% N_75 = nan(1,nbins);
% N_95 = nan(1,nbins);
% N_CI_w = nan(1,nbins);

plot_mean_N_of_D = 1;
if plot_mean_N_of_D == 1

figure('position',[100,400,500,500]);
thismap = rgb('darkorchid');

boundedline(d,N_med,N_CI_w,'-',...
    'alpha','nan','remove','cmap',thismap);
title(cruise);
ylabel('N(D)');
xlabel('D [mm]');
grid on; xlim([0.25 5]);
set(gca,'yscale','log');
set(gca,'TickDir','out','TickLength',[0.01 0.01]);

%=================
export_fig([plotdir cruise '_N_of_D_median_' version_num],'-png','-nocrop');
%=================

end

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

plot_all_N_of_D = 1;
if plot_all_N_of_D == 1

figure('position',[100,400,500,500]);
thismap = rgb('darkorchid');

plot([ 1 1], [1 1],'o','color',rgb('white'));
hold on;

for i = 1:nt
 pp = plot(d,ndrops(i,:),'.','color',thismap);
end
%m1 = plot(d, N_med,'-','color',rgb('magenta'),'linewidth',2);
%m2 = plot(d, N_mode,'-','color',rgb('lime'),'linewidth',2);
m3 = plot(d, N_mean,'-','color',rgb('black'),'linewidth',2);

legend([pp m3],'all points','mean');

title([cruise ' Drop Size Distributions']);

ylabel('N(D)');
xlabel('D [mm]');
grid on; xlim([0.25 7]);ylim([1 1000]);
set(gca,'yscale','log');
set(gca,'TickDir','out','TickLength',[0.01 0.01]);

%=================
export_fig([plotdir cruise '_N_of_D_all_' version_num],'-png','-nocrop');
%=================

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
drop_hist = histcounts(ndrops,d);

plot_count_N_of_D = 1;
if plot_count_N_of_D == 1

figure('position',[100,400,500,500]);

thismap = rgb('darkorchid');

first_good = find(N_count > 0,1,'first');
last_good = find(N_count > 0 ,1,'last');
N_count_f = N_count;
this_filtered_part = movmean(N_count(first_good:last_good),5,'omitnan');
N_count_f(first_good:last_good) = this_filtered_part;
N_count_f(last_good) = 1;

% plot(d,N_count,'-','color',rgb('red'),'LineWidth',2);
% hold on;
plot(d,N_count_f,'-','color',thismap,'LineWidth',2);

%plot(d,N_count,'-','color',thismap,'LineWidth',2);
title([cruise ' Drop Size Distributions']);

ylabel('total drop count per bin [#]');
xlabel('D [mm]');
grid on; xlim([0.25 8]);
set(gca,'yscale','log');
set(gca,'TickDir','out','TickLength',[0.01 0.01]);

% d1 = designfilt('lowpassiir','FilterOrder',12, ...
%     'HalfPowerFrequency',0.15,'DesignMethod','butter');
% N_count_filt = filtfilt(d1,N_count);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% subplot(1,2,2);
% plot(d,N_count/ndrops_total,'-','color',thismap,'LineWidth',2);
% title('SPURS-2 2016 Drop Size Distributions');
% 
% ylabel('Normalized Drop Count [%]');
% xlabel('D [mm]');
% grid on; xlim([0.25 8]);
% set(gca,'yscale','log');
% set(gca,'TickDir','out','TickLength',[0.01 0.01]);max(

%=================
export_fig([plotdir cruise '_N_of_D_count_' version_num],'-png','-nocrop');
%=================

end