clear all;
close all;
set_machine_et;
%%% Info: program to concatenate ODM .mat files from each day
%%%         run after ODM_save_PISTON.m

%%%     --- code for PISTON run by E. Thompson
%              Elizabeth Thompson
%              Applied Physics Lab, University of Washington
%              eliz@apl.uw.edu

%%%     --- code originally developed by K. Drushka
%              Kyla Drushka
%              Applied Physics Lab, University of Washington
%              kdrushka@apl.uw.edu

cruise = 'PISTON';
% cruise = 'SPURS2';
% cruise = 'TGT';

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

indir = [path_data '/processed/']; 
outdir = [path_data '/final/'];

version_num = 'v3';

%%%%% find files for each day

command = ['ls ' indir '*day*.mat > ' indir 'mat_files.txt'];
system(command);

files = char(textread([indir 'mat_files.txt'],'%s'));
[nfiles, nchars] = size(files);

process_files = 1;
if process_files == 1
    
    clear DSD;

for i = 1:nfiles
    
    filename = files(i,:);
    load(filename);

disp('**************************************************');
disp(['file # ' strnd(i,0) ' ' filename]);    
disp(['reading ' datestr(dsd.t(end),'YYYY mm - dd')]);

if i == 1
    first_time = dsd.t(1);
    
    % current readme of dsd files
%    ' t              = data and time, matlab format';
%    ' prate        = rain rate, mm/h (modified eq. 7);';...
%    ' voltage        = voltage, volts';...
%    ' urel_anem      = wind speed measured by cup anemometer attached to disdrometer, m/s';...
%    ' ndrops         = total number of drops in each bin, #';...
%    ' nsumrain              = total number of drops in all bins, #';...
%    ' psdrain      = same as n(d), drop count per volume, count per m^3 (Klepp eq. 2)
%    ' bin_nums       = number of bin, #';...
%    ' bin_sizes      = size of bin, mm';...
%    ' vfall          = estimated fall speed of drops per bin, m s^-1 (Gunn and Kinzer)';...
%    ' mass           = estimated mass of drops, kg (eq. 3)';...
    
    
DSD.README = [dsd.README; ...
    ' Dmax       = largest recorded drop';...
    ' Dmax_i     = index of largest recorded drop';...
    ' ntimes     = number of recorded times';...
    ' nbins      = number of bins';...
    ' dD         = size bin increments, mm';...
    ' refl       = reflectivity, dBZ';...
    ' reflfact   = reflectivity factor, mm^6/m^3';...
    ' lwc        = liquid water content, kg/m3';...
    ' r          = rain rate, mm/hr';...
    ' nd         = n(d) concentration, count per m^3 volume, same as psdrain';...
    ' Nd         = N(d) drop number density, count per mm size bin per m^3 volume';...
    ' odm_length = length of catchment area = 0.120 m = 120 mm';...
    ' odm_wdith  = width of catchment area = 0.022 m = 22 mm';...
    ' t_step     = time step of instrument = 60 sec'...
    ' P          = precip rate per time step, mm/hr';...
    ' Pd         = precip rate per bin per time step, mm/hr';...
    ' Pn         = precip rate per time step, mm/hr, using drop count and vfall instead of concentration';...
    ' Pdn        = precip rate per bin per time step, mm/hr, using drop count and vfall instead of concentration';...
    ' Prel       = precip rate per time step, mm/hr, using drop count, vfall, and rel wind instead of concentration';...
    ' Pdrel      = precip rate per bin per time step, mm/hr, using drop count, vfall, and rel wind instead of concentration';...
    ' W          = liquid water content, g/m3';...
    ' Wd         = liquid water content per bin, g/m3';...
    ' z          = reflectivity factor, mm^6/m^3';...
    ' zd         = reflectivity factor per bin, mm^6/m^3';...
    ' Z          = reflectivity, dBZ'... 
    ' Zd         = reflectivity per bin, dBZ'... 
    ];

    % all variables
    DSD.t               = dsd.t;            % data and time, matlab format
    DSD.praterain       = dsd.praterain;    % rain rate from read_odm Klepp for rain drops, mm/h, Eq 4
    DSD.pratesnow       = dsd.pratesnow;    % rain rate from read_odm Klepp for snow drops, mm/h, Eq 4
    DSD.prateall        = dsd.prateall;     % rain rate from read_odm Klepp for all drops, mm/h, Eq 4
    DSD.voltage         = dsd.voltage;      % voltage, volts
    DSD.urel_anem       = dsd.urel_anem;    % wind speed measured by disdrometer's attached cup anemometer, m/s
    DSD.nrain           = dsd.nrain;        % total number of rain drops in each bin, #
    DSD.nsnow           = dsd.nsnow;        % total number of snow drops in each bin, #
    DSD.nall            = dsd.nall;         % total number of all drops in each bin, #
    DSD.nsumrain        = dsd.nsumrain;     % total number of rain drops in all bins, #
    DSD.nsumsnow        = dsd.nsumsnow;     % total number of snow drops in all bins, #
    DSD.nsumall         = dsd.nsumall;      % total number of all drops in all bins, #
    DSD.psdrain         = dsd.psdrain;      % rain drop count per volume, count per m^3 (Klepp eq. 2)
    DSD.psdsnow         = dsd.psdsnow;      % snow drop count per volume, count per m^3 (Klepp eq. 2)
    DSD.psdall          = dsd.psdall;       % all drop count per volume, count per m^3 (Klepp eq. 2)
    DSD.bin_nums        = dsd.bin_nums;     % number of bin, #
    DSD.D               = dsd.D;            % sizes of bin, mm
    DSD.vfall           = dsd.vfall;        % estimated fall speed of drops per bin, m s^-1 (Klepp eq. 4)
    DSD.mass            = dsd.mass;         % estimated mass of drops, kg (Klepp eq. 3)
    % new var for this program:
    DSD.nbins           = length(DSD.D);    % num of bins

else  % variables discretized in time
    DSD.t            = [DSD.t           ; dsd.t];            
    DSD.praterain    = [DSD.praterain   ; dsd.praterain];            
    DSD.pratesnow    = [DSD.pratesnow   ; dsd.pratesnow];            
    DSD.prateall     = [DSD.prateall    ; dsd.prateall];            
    DSD.voltage      = [DSD.voltage     ; dsd.voltage];     
    DSD.urel_anem    = [DSD.urel_anem   ; dsd.urel_anem];         
    DSD.nrain        = [DSD.nrain       ; dsd.nrain];       
    DSD.nsnow        = [DSD.nsnow       ; dsd.nsnow];       
    DSD.nall         = [DSD.nall        ; dsd.nall];       
    DSD.nsumrain     = [DSD.nsumrain    ; dsd.nsumrain];            
    DSD.nsumsnow     = [DSD.nsumsnow    ; dsd.nsumsnow];            
    DSD.nsumall      = [DSD.nsumall     ; dsd.nsumall];            
    DSD.psdrain      = [DSD.psdrain     ; dsd.psdrain]; 
    DSD.psdsnow      = [DSD.psdsnow     ; dsd.psdsnow]; 
    DSD.psdall       = [DSD.psdall      ; dsd.psdall]; 
end

    %%% new vars in this program
    DSD.ntimes  = length(DSD.t);


if i == nfiles
    last_time = dsd.t(end);
end

end  %%% loop for files


%% which var do you want to use for prate and integral rain parameters?
ndrops = DSD.nrain;
nsumdrops = DSD.nsumrain;


%%  rest of the stats
%%% figure out more statistics of DSD
%%% numbers
nbins=DSD.nbins;
D = DSD.D;
nt = DSD.ntimes;

DSD.Dmax    = nan(DSD.ntimes,1);
DSD.Dmax_i  = nan(DSD.ntimes,1);

for i = 1:DSD.ntimes
    %%% find the last bin that recorded a drop
    largest_drop_bin = find(ndrops(i,:) > 0, 1, 'last');
    if ~isempty(largest_drop_bin)
        DSD.Dmax(i) = DSD.D(largest_drop_bin);
        DSD.Dmax_i(i) = largest_drop_bin;
    end
end


%%%% derived variables:



%%% width of size bins in millimeters
% dD = diff(D); 
midpoints = D;
edges = nan(1,nbins+1);
for i = 2:nbins
    if i > 1 && i < nbins+1
        %%% assum edges are the midpoint between each D center
        edges(i) = sum([D(i) D(i-1)])/2;
    elseif i == nbins+1
    end
end
% now go back and do first and last edge: #1 and #nbins+1
edges(1)        = D(1) - (edges(2) - D(1));
edges(nbins+1)  = D(nbins) + (D(nbins) - edges(nbins)); 

% figure;
% plot(edges,zeros(nbins+1,1),'x');
% hold on;
% plot(edges(1), 0,'ok','markersize',10);
% plot(D(1), 0,'og','markersize',10);
% plot(edges(end), 0,'ok','markersize',10);
% plot(D(end), 0,'og','markersize',10);
% plot(D,zeros(nbins,1),'.');
% legend('edges','midpoint D');
% grid on;

dD = diff(edges);

%%% note, fortran code does this a bit differently:
%       dmm_low(nr(i))=EXP((float(nr(i))/94.)*ALOG(10.))-1.
%       dmm_upp(nr(i))=EXP(((float(nr(i))+1.)/94.)*ALOG(10.))-1.
%       dmm_interval(nr(i))=dmm_upp(nr(i))-dmm_low(nr(i))
%       dmm_mean(nr(i))=(dmm_low(nr(i))+dmm_upp(nr(i)))/2.
%       dsd_refl(nr(i))= n_a(nr(i))*(dmm_mean(nr(i))**6)
%       dsd_refl_sum=dsd_refl_sum+dsd_refl(nr(i))
%          if(dsd_refl_sum.ne.0.)dbz=10*(ALOG10(dsd_refl_sum))
%          if(dsd_refl_sum.eq.0.)dbz=-32.5  
%          if(rra.ne.0.)dbr=10*(ALOG10(rra))
%          if(rra.eq.0.)dbr=-32.5  
%       write(71,'(1x,i5,4i3,3(f10.4,x),f14.4,f6.0,f12.4,f5.2,i3.3,
%      *1x,f12.2,1x,f12.2,1x,f12.2)')
%      *iyr,imth,id,ihr,imin,ff,rra,rrds1,dd,summeanz,rrzeit
%      *,uref,ianzbin_r,dsd_refl_sum,dbz,dbr
     
%%%% calcuations
Nd = nan(nt, nbins);    % number density N(d): count per mm size bin per m3 volume... called psd number concentration in Klepp 2015
nd = nan(nt, nbins);    % concentration n(d): count per m3... called n(bin) = psd density in Klepp 2015
Pd = nan(nt, nbins);    % rain rate per bin mm/hr
P = nan(nt, 1);         % rain rate all bins mm/hr
Pdn = nan(nt, nbins);   % rain rate per bin mm/hr
Pn = nan(nt, 1);        % rain rate all bins mm/hr
Pdrel = nan(nt, nbins); % rain rate per bin mm/hr
Prel = nan(nt, 1);      % rain rate all bins mm/hr
Wd = nan(nt, nbins);    % liquid water content per bin g/m3
W = nan(nt, 1);         % liquid water content all bins g/m3
zd = nan(nt, nbins);    % reflecitivity factor per bin mm6/m3
z = nan(nt, 1);         % reflecitivity factor all bins mm6/m3

odm_length  = 0.120;    %% 120 mm so 0.120 m
odm_width   = 0.022;  	%% 22 mm so 0.022 m
t_step      = 60;       %% 60 second time integration

%%% weird parts of fortran code that have some strange geometry and embedded weird functions...
% dxu=evaldxu(rmmq(i),amm2(j+1),uradmm)
% ws=(uradmm+rmmq(i)-dxu-alt)/(uradmm+rmmq(i))
% alt=uradmm+rmmq(i)-dxu
% if(rnanz(j).gt.(float(n(i))*ws))then
%    rnanz(j)=rnanz(j)-(float(n(i))*2.*ws)

% calculate drop number density from drop counts, N(D) in per mm per m3.
% This is different than n(bin) from Eq 3 in Christian's paper.
% dD is the width of each size bin in mm. urel_anem and fall_v are both in
% m/s. ndrops is just a count. 
for i = 1:nt
    for j = 1:nbins
%         if j < nbins  %%% this was used back when I didn't know how to do
%         the bin widths properly.
            
            Nd(i,j) = ndrops(i,j) / ...
                (dD(j)*odm_width*odm_length*t_step* sqrt(DSD.urel_anem(i)^2 + DSD.vfall(j)^2));
            nd(i,j) = ndrops(i,j) / ...
                      (odm_width*odm_length*t_step* sqrt(DSD.urel_anem(i)^2 + DSD.vfall(j)^2));

%             M(j) = 4/3*pi*1000*(D(j)./2000).^3; mass per bin size in kg
            Pd(i,j) = 3600* nd(i, j).*DSD.vfall(j).*DSD.mass(j); % rain rate per bin per time step, mm/hr
            Pdn(i,j) = 3600* ndrops(i, j).*DSD.vfall(j).*DSD.mass(j); % rain rate per bin per time step, mm/hr
            Pdrel(i,j) = 3600* ndrops(i, j).*sqrt(DSD.vfall(j)^2 + DSD.urel_anem(i)^2).*DSD.mass(j); % rain rate per bin per time step, mm/hr
            Wd(i,j) = 1000* nd(i,j) .* DSD.mass(j); %%% lwc per bin per time step, g/m3
            zd(i,j) = nd(i,j) .* D(j)^6; %%% reflectivity factor per bin per time step, mm6/m3
            
%             Wd(i,j) = 1000* ndrops(i,j) .* DSD.mass(j) / ... 
%                       (odm_width * odm_length * t_step * sqrt(DSD.urel_anem(i)^2 + DSD.vfall(j)^2));  %%% lwc per bin per time step, kg/m3
%             zd(i,j) = ndrops(i,j) .* D(j).^6 / ... 
%                       (odm_width * odm_length * t_step * sqrt(DSD.urel_anem(i)^2 + DSD.vfall(j)^2));  %%% lwc per bin per time step, kg/m3

%         else
%             Nd(i,j) = 0;
%             nd(i,j) = 0;
%             Pd(i,j) = 0;
%             Pdn(i,j) = 0;
%             Pdrel(i,j) = 0;
%             Wd(i,j) = 0;
%             zd(i,j) = 0;
%         end
    end
    P(i) = nansum(Pd(i,:));  % rain rate in mm/hr of all bins per time step, using n(d) concentration as specified in Klepp 2015
    Pn(i) = nansum(Pdn(i,:));  % rain rate in mm/hr of all bins per time step, using drop count instead of concentration
    Prel(i) = nansum(Pdrel(i,:));  % rain rate in mm/hr of all bins per time step, using rel velocity and vfall
    W(i) = nansum(Wd(i,:));  % liquid water content in g/m3 of all bins per time step
    z(i) = nansum(zd(i,:));        % reflectivity factor in mm6/m3 of all bins per time step
    
end

% Mparticle=4/3*pi*1000*(D./2000).^3;     % mass of drops, in kg (Klepp eq. 3 - modified! should be 2000 not 200 in the denominator, to convert from diam[mm] to radius[m])
% Prate = nansum(3600 * nrain.*vfall.*Mparticle); % rain rate in mm/hr (Klepp eq. 6)
% Prate = nansum(3600 * ndrops.*vfall.*Mparticle); % rain rate in mm/hr (Klepp eq. 6)

%%% reset nans to zeros *** this was super important for Nw Do processing
Nd(find(isnan(Nd) == 1)) = 0;
nd(find(isnan(nd) == 1)) = 0;

DSD.edges = edges;
DSD.dD = dD;
DSD.Nd = Nd;
DSD.nd = nd;
DSD.Pd = Pd;
DSD.P = P;
DSD.Pdn = Pdn;
DSD.Pn = Pn;
DSD.Pdrel = Pdrel;
DSD.Prel = Prel;
DSD.Wd = Wd;
DSD.W = W;
DSD.zd = zd;
DSD.z = z;
DSD.Zd = 10*log10(zd);
DSD.Z = 10*log10(z);
DSD.odm_length = odm_length; % = 0.120;    %% 120 mm so 0.120 m
DSD.odm_width  = odm_width; % = 0.022;  	%% 22 mm so 0.022 m
DSD.t_step     = t_step; % = 60;       %% 60 second time integration

% the optical disdrometer catchment area is not constant because it
% actually takes relative wind speed and residence time into account
%%% NOTE: average area seems to be nearly 2.4 x 10^6 mm^3, or 0.0024 m^3.
%%% Small! 

% ? = dsd.ndrops./dsd.Nd; I don't konw what this calculates

% 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

wh_no_drops = find(nsumdrops == 0);
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;

% make a rain estimate if there were sufficient data for the row
r           = nan(nt,1);   % rain rate in mm/hr
reflfact    = nan(nt,1);   % z in mm6 m-3: radar reflectivity factor
lwc         = nan(nt,1);   % lwc in g m-3: liquid water content

for i=1:nt
  if(n_check(i)==1)   
    r(i)        = 6 * pi * 1e-4 *   trapz(D, ndrops(i,1:nbins) .* D.^3 .* DSD.vfall); %%% the drops didn't just fall through, as this equation suggests. 
    reflfact(i) =                   trapz(D, ndrops(i,1:nbins) .* D.^6 );
    lwc(i)      = pi/6 * 1e-3 *     trapz(D, ndrops(i,1:nbins) .* D.^3 );
  end
end

% Z in dBZ: radar reflectivity
refl = 10*log10(reflfact);

DSD.refl = refl;
DSD.reflfact = reflfact;
DSD.lwc = lwc;
DSD.r = r;

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

time_stamp = [datestr(first_time,'YYYY_mmdd') '_to_' datestr(last_time,'YYYY_mmdd')];
filename = [outdir cruise '_dsd_level_1_all_data_' version_num '.mat'];

save(filename, 'DSD');


end %%% if processing files

these = find(DSD.t <= datenum(2018,9,20,0,0,0) & ...
    DSD.t >= datenum(2018,9,17,0,0,0));



