clear all;
close all;
set_machine_et;

%%% Info: program to read in ODM data, calculate rain rate, 
%%%         fix time on a few files, and output data files for sharing

%%%     --- 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 '/raw/']; 
outdir = [path_data '/processed/']; 

%%%%% find files for each day

command1 = ['ls ' indir 'ODM*rd* > ' indir 'rd_files.txt'];
system(command1);

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

process_files = 1;
if process_files == 1

%%% NOTE: will fail for empty files of 0 bytes. 
%    for i = 8:8
for i = 1:nfiles
    filename = files(i,:);
    
    s = dir(filename);
    if s.bytes == 0
        disp(['moving empty file to folder: ' s.name]);
        mv_empty_file_cmd = ['mv ' filename ' ' indir 'empty_files/'];
        system(mv_empty_file_cmd);
    else
%     disp(filename);
    clear dsd;
    dsd.README          = {[' Readme for ' cruise ' drop size distribution data'];...
                       '';...
                       ' Please acknowledge Elizabeth J. Thompson and Kyla Drushka';...
                       ' from the Applied Physics Laboratory at the University of Washington';...
                       ' elizabeth.thompson@noaa.gov, kdrushka@apl.uw.edu';...
                       '' ;...
                       ' Please also cite Klepp 2015 regarding instrument and calculations';...
                       ' see below for citation';...
                       '';...
                       ' The ODM-470 is an optical disdrometer developed by Christian Klepp';...
                       '     (Klepp 2015). It was specifically designed for ships in that it';...
                       '     rotates freely with the wind using the attached anemometer.';...
                       '     The anemometer wind speed is then used to calculate the drop density,';...
                       '     rain rate, etc... (see Klepp 2015).';...
                       '';...
                       ' NOTES:';...
                       '  - due to noise or imperfections, drops are often detected by unit ';...
                       '      even when it is definitely not raining. These do not produce a ';...
                       '      significant rain rate. Suggest thresholding data by a minimum ';...
                       '      number of total drops per time step, i.e. with the N variable.';...
                       '  - drops are counted each minute';...
                       '  - drop density is calculated as count per mm bin size per m^3 volume,';...
                       '      so the resulting units of drop_density are mm^-1 m^-3';...
                       '  - the first 14 size bins, below 0.44 mm, are ignored because drops';...
                       '      are too small to produces voltages higher than the noise floor';...
                       '  - 128 size bins total, which are logarithmically spaced';...
                       '  - some quantities were calculated with equations from Klepp 2015';...
                       '';...
                       ' For complete information about the ODM disdrometer and calculations,';...
                       ' see Klepp 2015: ';...
                       '';...
                       ' Klepp, C., 2015: The oceanic shipboard precipitation ';...
                       ' measurement network for surface validation- OceanRAIN. ';...
                       ' Atmos. Res., 163, 74?90.';...
                       ' http://dx.doi.org/10.1016/j.atmosres.2014.12.014';...
                       '-----------------------------------------------------------------------';...
                       '';...
                       ' matlab structure = dsd';...
                       '';...
                       ' see end of program ODM_save.m for explaination of contents:';...
                       'dsd.praterain        % rain rate from nrain, mm/h (Klepp eq. 7);';...
                        'dsd.pratesnow       % rain rate from nsnow, mm/h (Klepp eq. 7);';...
                        'dsd.prateall        % rain rate from nall (rain+snow), mm/h (Klepp eq. 7);';...
                        'dsd.voltage         % voltage, volts';...
                        'dsd.urel_anem       % rel wind speed measured by cup anemometer at ODM, m/s';...
                        'dsd.nrain           % total number of rain drops in each bin, #';...
                        'dsd.nsnow           % total number of snow drops in each bin, #';...
                        'dsd.nall            % total number of all (rain+snow) drops in each bin, #';...
                        'dsd.nsumrain        % total number of rain drops in all bins, #';...
                        'dsd.nsumsnow        % total number of snow drops in all bins, #';...
                        'dsd.nsumall         % total number of all (rain+snow) drops in all bins, #';...
                        'dsd.psdrain         % rain drop count per bin, count per m^3  (Klepp eq. 2)';...
                        'dsd.psdsnow         % snow drop count per bin, count per m^3 (Klepp eq. 2)';...
                        'dsd.psdall          % all drop (rain + snow) count per bin, count per m^3 (Klepp eq. 2)';...
                        'dsd.bin_nums        % number of bin, #';...
                        'dsd.D               % sizes of bin, mm';...
                        'dsd.vfall           % estimated fall speed of drops per bin, m s^-1 (Klepp eq. 4, but for mm))';...
                        'dsd.mass            % estimated mass of drops, kg (Klepp eq. 3, but for mm)';...
                       };
                   
%%%%% run the function that reads the raw files
out = read_odm_fix_DSD_et(filename);

% if out.sumall > out.sumrain
%     disp('sum(all) > sum(rain)... meaning drops were classified as snow');
%     disp(['sum(all) = ' sprintf('%i',out.sumall)]);
%     disp(['sum(snow) = ' sprintf('%i',out.sumrain)]);
%     disp(['sum(rain) = ' sprintf('%i',out.sumsnow)]);
% elseif out.sumsnow > 0out
%     disp('sum(snow) > 0 ... meaning drops were classified as snow');
% end

%%% parameterized fall velocity [m/s] per drop diameter [mm] from Gunn & Kinzer 1949
% fall_v / out.vfall = 0.7738
% out.vfall = 0.7738 * fall_v
%%% this seems wrong. Redo-ing it now.
% vfall = 9.65-10.3*exp(-.6.*out.D_p);

%%% this reference voltage of bad data should be checked for each time
%%% series, just in case. It also maybe shouldn't be used at all...
% bad_rain = find(out.voltage < 1);
% out.prate(bad_rain) = 0.0;

disp(['reading ' datestr(out.t(1),'YYYY mm - dd')]);

if strcmp(cruise,'PISTON') == 1
    %for files before 2018 09-18 05:16:54, subtract one hour to time
    hr1 = datenum(2018,9,18,1,0,0) - datenum(2018,9,18,0,0,0);
    min1 = datenum(2018,9,18,0,1,0) - datenum(2018,9,18,0,0,0);
    bad_time_1   = datenum(2018, 9, 18,  5, 17,  0);
    bad_time_2       = datenum(2018, 9, 28,  6,  0,  0);
    bad_time_3      = datenum(2018, 9, 29, 19, 16, 11);
    bad_time_4        = datenum(2018, 9, 30, 0, 0, 0);

    %if i <= 9

%%% I set the computer with the wrong UTC time when installed. 
%%%      It was originally too fast by one hour. Subtract 1 hour.
    out.tnew = out.t;
    out.tnew(out.t < bad_time_1) = out.t(out.t < bad_time_1)-hr1;

    %%% accidentally turned off the computer and lost a minute. Add 1 minute.
    out.tnew(out.t > bad_time_1 & out.t <= bad_time_2) = ...
        out.t(out.t > bad_time_1 & out.t <= bad_time_2)+min1;

    %%% accidentally turned off the computer again and lost another minute. Add 2 minutes.
    out.tnew(out.t > bad_time_2 & out.t < bad_time_3) = ...
        out.t(out.t > bad_time_2 & out.t < bad_time_3)+(2*min1);

    %%% Fixing everything beyond this point. Subtract 3 total minutes.
    out.tnew(out.t > bad_time_4) = ...
        out.t(out.t > bad_time_4)-(3*min1);


    %end

    % %first_duplicated_time = datenum(2018,9,18,5, 16, 54);
    % if out.t < last_bad_time 
    %     out.tnew  = out.t - hr1;
    % end

    % duplicated_data = find(out.t >= last_bad_time & out.t <= last_bad_time+hr1);

    % %%% then get rid of duplicated data
    % out.prate(duplicated_data) = 0.0;
    % out.voltage(duplicated_data) = nan;
    % out.nrain(duplicated_data,1:128) = zeros(1,128);
    % out.psd(duplicated_data,1:128) = zeros(1,128);
    % out.sumrain(duplicated_data) = 0.0;
    % out.urel_anem(duplicated_data) = 


    %%%%% save data for each day: output from each minute
    % if i <= 9
    
    dsd.t           = out.tnew;            % data and time, matlab format
else
    dsd.t = out.t;
end
% else
%     dsd.t           = out.t;
% end

dsd.praterain       = out.praterain;    % rain rate from nrain, mm/h (Klepp eq. 7);
dsd.pratesnow       = out.pratesnow;    % rain rate from nsnow, mm/h (Klepp eq. 7);
dsd.prateall        = out.prateall;     % rain rate from nall (rain+snow), mm/h (Klepp eq. 7);
dsd.voltage         = out.voltage;      % voltage, volts
dsd.urel_anem       = out.urel_anem;    % rel wind speed measured by cup anemometer at ODM, m/s
dsd.nrain           = out.nrain;        % total number of rain drops in each bin, #
dsd.nsnow           = out.nsnow;        % total number of snow drops in each bin, #
dsd.nall            = out.nall;         % total number of all (rain+snow) drops in each bin, #
dsd.nsumrain        = out.nsumrain;     % total number of rain drops in all bins, #
dsd.nsumsnow        = out.nsumsnow;     % total number of snow drops in all bins, #
dsd.nsumall         = out.nsumall;      % total number of all (rain+snow) drops in all bins, #
dsd.psdrain         = out.psdrain;      % rain drop count per bin, count per m^3  (Klepp eq. 2)
dsd.psdsnow         = out.psdsnow;      % snow drop count per bin, count per m^3 (Klepp eq. 2)
dsd.psdall          = out.psdall;       % all drop (rain + snow) count per bin, count per m^3 (Klepp eq. 2)
dsd.bin_nums        = out.binnum;       % number of bin, #
dsd.D               = out.D_p;          % sizes of bin, mm
dsd.vfall           = out.vfall;        % estimated fall speed of drops per bin, m s^-1 (Klepp eq. 4, but for mm))
dsd.mass            = out.Mparticle;    % estimated mass of drops, kg (Klepp eq. 3, but for mm)

%%% don't save with the dsd.t(1) because it is actually from the previous
%%% day for the first 9 files!!
filename = [outdir cruise '_dsd_day_' datestr(out.t(1),'YYYY_mmdd') '.mat'];
disp(['saving ' datestr(out.t(1),'YYYY mm - dd')]);

save(filename, 'dsd');

    end %%% if file is not empty
end %%% for each file

end %%% if processing files