clear all;
close all;
set_machine_et;

% volume and cross sectional area
% K15 says cross sectional area is 26.4 cm2, volume is 45.6 cm3
% l = 0.120; % meters... length/height of cylinder
% w = 0.022; % meters... diameter of cylinder
% h = 0.011; % meters... assume half of the diameter
% 
% v_cyl = pi*l*((0.5*w)^2);
% v_rect = l*w*h;
% a = l*w;
% v_rect_cm = v_rect*10^6;
% v_cyl_cm = v_cyl*10^6;
% a_cm = a*10^4;


%%% the fortran code uses a 0.01 threshold of prate. 
%         if(rru.lt.0.01)rru=0.0
%         if(rra.lt.0.01)rra=0.0
%         if(rra.gt.99999.)rra=0.0


% Code to intake ODM data from my code vs. Kyla/CHristian's code
% for SPURS-2 and PISTON deployments

% cruise = 'TGT';         cruise_title = 'TGT 2017 - 2018';   version_num = 'v3';
cruise = 'PISTON';    cruise_title = 'PISTON 2018';       version_num = 'v3';
% cruise = 'SPURS2';    cruise_title = 'SPURS2 2016-2017';  version_num = 'v3';

path_data = ['/Users/eliz/Dropbox/DSD/' cruise '/'];
syndir = '/Users/eliz/Dropbox/DSD/synthesis/';
indir = [path_data 'final/']; 
outdir = [path_data 'final/'];
plotdir = [path_data 'analysis/'];
metdir = [path_data 'ship_met/'];

%%% ET dsd data
% load([indir cruise '_dsd_level_2_v1.mat']); % Nw, D0, mu, and everything on that time base of "good" retrievals
load([indir cruise '_dsd_level_1_all_data_' version_num '.mat']);  % everything else
et_xlim = [min(DSD.t) max(DSD.t)];


%%% files to compare to
if strcmp(cruise,'TGT') == 1

    kfile_2 = '/Users/eliz/Dropbox/DSD/TGT/processed_Kyla_Apr2020/odm_dsd_tgt_mod.mat';
        if isfile(kfile_2) == 0
            k_p18 = load('/Users/eliz/Dropbox/DSD/TGT/processed_Kyla_Apr2020/odm_dsd_piston.mat');
            k = k_p18.dsd;
        else
            load(kfile_2)    
        end
%     load([metdir 'PISTON_met_1min.mat']); %%% PISTON 2018 1-min met data
    R_m = nan(3600,1);
    t_m = nan(3600,1);
    yd_m = nan(3600,1);
    went_bad_on_TGT = datenum(2018, 5, 15, 0, 0, 0);
    test_case = 42211;
elseif strcmp(cruise,'SPURS2') == 1
    kfile_2 = '/Users/eliz/Dropbox/DSD/SPURS2/processed_Kyla_Apr2020/odm_dsd_spurs2_mod.mat';
        if isfile(kfile_2) == 0
            k_s16 = load('/Users/eliz/Dropbox/DSD/SPURS2/processed_Kyla_Apr2020/odm_dsd_spurs2_2016.mat');
            k_s17 = load('/Users/eliz/Dropbox/DSD/SPURS2/processed_Kyla_Apr2020/odm_dsd_spurs2_2017.mat');
            k = cat_2mats(k_s16.dsd, k_s17.dsd);
            k.README = [k_s16.dsd.README ; k_s17.dsd.README];
        else
            load(kfile_2)    
        end
        k.w = k.t*nan;
    m_s16 = load([metdir 'SPURS2_2016_OneMinute_Met_Rev3']); %%% SPURS-2 2016 met data
    m_s17 = load([metdir 'SPURS2_2017_OneMinute_Met_Rev3']); %%% SPURS-2 2017 met data
    R_m = [m_s16.Prate  m_s17.Prate]';
    yd_m = [m_s16.Yday  m_s17.Yday]';
    t_m = [m_s16.Yday + datenum(2016,0,0,0,0,0)  m_s17.Yday + datenum(2017,0,0,0,0,0)]';
    test_case = 40593; 
elseif strcmp(cruise,'PISTON') == 1
    kfile_2 = '/Users/eliz/Dropbox/DSD/PISTON/processed_Kyla_Apr2020/odm_dsd_piston_mod.mat';
        if isfile(kfile_2) == 0
            k_p18 = load('/Users/eliz/Dropbox/DSD/PISTON/processed_Kyla_Apr2020/odm_dsd_piston.mat');
            k = k_p18.dsd;
        else
            load(kfile_2)    
        end
    load([metdir 'PISTON_2018_flux_met_sea_1min_r1.mat']); %%% PISTON 2018 1-min met data
    R_m = b.prate;
    t_m = b.t;
    test_case = 9237;
end

if strcmp(cruise,'PISTON') ~=1
    R_k         = k.p;
    reflfact_k  = k.reflectivity;
    t_k         = k.t;
    psd_k       = nan(3600, 1);
    refl_k      = 10*log10(reflfact_k); 
else  %%% if cruise == PISTON, dsd from Kyla are missing
    R_k         = nan(3600,1);
    t_k         = nan(3600,1);
    reflfact_k  = nan(3600,1);
    refl_k      = nan(3600,1);
    psd_k       = nan(3600,1);
    k.ndrops    = nan(3600,1);  
end
    nk = length(t_k);

% Kyla's data:  
% 
%   dsd = struct with fields:
% 
%                t: [8565×1 double]
%             wspd: [8565×1 double]  %% from anemometer
%         sumdrops: [8565×1 double]
%           ndrops: [8565×128 double]
%              psd: [8565×128 double]
%              D_p: [1×128 double]
%        dropclass: [1×128 double]
%                p: [8565×1 double]
%                w: [8565×1 double]   ... fall velocity?
%             Vref: [8565×1 double]
%     reflectivity: [8565×1 double]


%%% if this hasn't already bee done before.... 

if isfield(k, 'ndrops') == 1 && isfield(k, 'Dmax') == 0
    disp('making Dmax and nsumrain fields');
    nk = length(t_k);
    Dmax_k = nan(nk,1);
    Dmax_i_k = nan(nk,1);
    nsumrain_k = nan(nk,1);
    
    %%% do not do calculations if kyla's dataset is missing and therfore nans
    if isnan(t_k(300)) == 0
        disp('time data are not nans, so go ahead');
        % total number of drops in each time step
        nsumrain_k = sum(k.ndrops,2);

        for i = 1:length(t_k)
            %%% find the last bin that recorded a drop
            largest_drop_bin = find(k.ndrops(i,:) > 0, 1, 'last');
            if ~isempty(largest_drop_bin)
                Dmax_k(i) = DSD.D(largest_drop_bin);
                Dmax_i_k(i) = largest_drop_bin;
            end
        end

        % save new fields to k and save k as new file
        disp('saving new kyla file');
        k.Dmax = Dmax_k;
        k.Dmax_i = Dmax_i_k;
        k.nsumrain = nsumrain_k;
        save(kfile_2,'k');
    end
    
else
    disp('dmax already exists, so did not recalculate');
    Dmax_k = k.Dmax;
    Dmax_i_k = k.Dmax_i;
    nsumrain_k = k.nsumrain;
end


figure;
subplot(1,2,1); hold on;
plot(t_k, Dmax_k,'xb');
plot(DSD.t, DSD.Dmax,'oc');
legend('Kyla','ET','location','best');
ylabel('D_{MAX}, mm');
grid on;
title(cruise_title);
% ylim([0 ]);
xlim(et_xlim);
datetick('x','mm/dd','keeplimits');
if strcmp(cruise,'TGT') == 1
    plot([went_bad_on_TGT went_bad_on_TGT], ylim,'--r','linewidth',2);
    legend('Kyla','ET','ODM went bad May 15, 2018','location','best');
end

subplot(1,2,2); hold on;
plot(t_k, nsumrain_k,'xb');
plot(DSD.t, DSD.nsumrain,'oc');
legend('Kyla','ET','location','best');
ylabel('total drop count');
grid on;
title(cruise_title);
% ylim([0 250]);
xlim(et_xlim);
datetick('x','mm/dd','keeplimits');
if strcmp(cruise,'TGT') == 1
    plot([went_bad_on_TGT went_bad_on_TGT], ylim,'--r','linewidth',2);
    legend('Kyla','ET','ODM went bad May 15, 2018','location','best');
end

export_fig([plotdir cruise '_test_odm_Dmax_ndrops_' version_num],'-png','-nocrop','-r300');

%% make comparison plots and stats:

% psd, Nd, nd, drop_density_Klepp

figure;
subplot(1,2,1); hold on;
plot(DSD.D(1:40), DSD.nd(test_case, 1:40),'oc'); 
plot(DSD.D(1:40), DSD.psdrain(test_case, 1:40),'v','color',rgb('purple')); 
plot(DSD.D(1:40), DSD.psdsnow(test_case, 1:40),'v','color',rgb('orange')); 
plot(DSD.D(1:40), DSD.nrain(test_case, 1:40),'x','color',rgb('purple'));
plot(DSD.D(1:40), DSD.nsnow(test_case, 1:40),'x','color',rgb('orange'));
legend('ET n(d) rain','Klepp psdrain','Klepp psdsnow','Klepp # nrain','Klepp # nsnow');
xlabel('D (mm)');
ylabel('count, # ..... or ....  n(d), # m^{-3}');
grid on;
title(cruise_title);

subplot(1,2,2); hold on;
plot(DSD.D(1:40), DSD.Nd(test_case, 1:40),'oc'); 
legend('ET');
xlabel('D (mm)');
ylabel('N(D), # mm^{-1} m^{-3}');
grid on;
title(cruise_title);

export_fig([plotdir cruise '_test_odm_nd_Nd_' version_num],'-png','-nocrop','-r300');


% R_k, R_m, r_Klepp, r
figure; hold on;
plot(t_m, R_m, '-', 'color',rgb('gray'),'linewidth',2);
plot(t_k, R_k,'xb','linewidth',4);
legend('ship optical gauge','Kyla DSD','location','best');
grid on;
ylim([0 250]);
xlim(et_xlim);
datetick('x','mm/dd','keeplimits');
title(cruise_title);
ylabel('P, mm hr^{-1}');
if strcmp(cruise,'TGT') == 1
    plot([went_bad_on_TGT went_bad_on_TGT], ylim,'--r','linewidth',2);
    legend('ship optical gauge','Kyla DSD','ODM went bad May 15, 2018','location','best');
end
export_fig([plotdir cruise '_test_odm_R_kyla_' version_num],'-png','-nocrop','-r300');

figure; hold on;
plot(DSD.t, DSD.pratesnow,'v','color',rgb('orange'),'linewidth',4);
plot(t_m, R_m, '-', 'color',rgb('gray'),'linewidth',2);
plot(DSD.t, DSD.praterain,'v','color',rgb('purple'),'linewidth',4);
legend('Klepp prate snow','ship optical gauge','Klepp prate rain','location','best');
grid on;
ylim([0 250]);
xlim(et_xlim);
datetick('x','mm/dd','keeplimits');
title(cruise_title);
ylabel('P, mm hr^{-1}');
if strcmp(cruise,'TGT') == 1
    plot([went_bad_on_TGT went_bad_on_TGT], ylim,'--r','linewidth',2);
    legend('Klepp prate snow','ship optical gauge','Klepp prate rain','ODM went bad May 15, 2018','location','best');
end
export_fig([plotdir cruise '_test_odm_R_klepp_' version_num],'-png','-nocrop','-r300');

figure; hold on;
plot(t_m, R_m, '-', 'color',rgb('gray'),'linewidth',2);
plot(DSD.t, DSD.r,'oc','linewidth',4);
plot(DSD.t, DSD.P,'o','color',rgb('green'),'linewidth',4);
plot(DSD.t, DSD.Prel,'o','color',rgb('blue'),'linewidth',4);
plot(DSD.t, DSD.Pn,'o','color',rgb('purple'),'linewidth',4);
legend('ship optical gauge','ET integral','ET n(d) - BEST','ET urel&count','ET vfall&count','location','best');
grid on;
ylim([0 250]);
xlim(et_xlim);
datetick('x','mm/dd','keeplimits');
title(cruise_title);
ylabel('P, mm hr^{-1}');
if strcmp(cruise,'TGT') == 1
    plot([went_bad_on_TGT went_bad_on_TGT], ylim,'--r','linewidth',2);
    legend('ship optical gauge','ET integral','ET n(d) - BEST','ET urel&count','ET vfall&count','ODM went bad May 15, 2018','location','best');
end
export_fig([plotdir cruise '_test_odm_R_et_versions_' version_num],'-png','-nocrop','-r300');

figure; hold on;
plot(t_m, R_m, '-', 'color',rgb('gray'),'linewidth',2);
plot(DSD.t, DSD.P,'oc','linewidth',4);
legend('ship optical gauge','ET n(d)','location','best');
grid on;
ylim([0 250]);
xlim(et_xlim);
datetick('x','mm/dd','keeplimits');
title(cruise_title);
ylabel('P, mm hr^{-1}');
if strcmp(cruise,'TGT') == 1
    plot([went_bad_on_TGT went_bad_on_TGT], ylim,'--r','linewidth',2);
    legend('ship optical gauge','ET n(d)','ODM went bad May 15, 2018','location','best');
end
export_fig([plotdir cruise '_test_odm_R_et_' version_num],'-png','-nocrop','-r300');




% DSD.P_int_k = interp1(DSD.t, DSD.P, k.t,'nearest');

if strcmp(cruise,'PISTON') ~= 1
% find times that my time base (from raw files) overlaps with Kyla's
[ye, me, de, hhe, mme, sse] = datevec(DSD.t);
[yk, mk, dk, hhk, mmk, ssk] = datevec(t_k);
t_e_min = datenum(ye, me, de, hhe, mme, 0);
t_k_min = datenum(yk, mk, dk, hhk, mmk, 0);


time_et = round(t_e_min.*86400);
time_k = round(t_k_min.*86400);
noint_k = find(ismember(time_et, time_k) == 0);
int_k = find(ismember(time_et, time_k) == 1);
    

figure('position',[1,1,1200,400]);
subplot(1,3,1);
hold on;
plot(k.p, DSD.P(int_k),'o')
axis square
grid on;
plot([0 100],[0 100]);
xlabel('kyla');ylabel('ET'); 
xlim([0 100]); ylim(xlim);
title([cruise_title ' P [mm/h]']);

subplot(1,3,2);
hold on;
plot(k.w, DSD.W(int_k),'o')
axis square
grid on;
plot([0 6],[0 6]);
xlabel('kyla');ylabel('ET'); 
xlim([0 6]); ylim(xlim);
title([cruise_title ' W [g/m3]']);

subplot(1,3,3);
hold on;
plot(refl_k, DSD.Z(int_k),'o')
axis square
grid on;
plot([0 60],[0 60]);
xlabel('kyla');ylabel('ET'); 
xlim([0 60]); ylim(xlim);
title([cruise_title ' Z [dBZ]']);


export_fig([plotdir cruise '_test_odm_R_et_kyla_' version_num],'-png','-nocrop','-r300');

end

figure; 
hold on;
if strcmp(cruise,'SPURS2') == 1
    plot(k.p, DSD.P,'.')
elseif strcmp(cruise,'TGT') == 1
    plot(k.p, DSD.P(int_k),'.')
end
plot(DSD.praterain, DSD.P,'x')
axis square
plot([0 100],[0 100]);
xlabel('Klepp');ylabel('ET'); 
xlim([0 100]); ylim(xlim);
if strcmp(cruise,'PISTON') ~=1
    legend('Kyla Klepp v. ET n(d)','ET Klepp vs. ET n(d)', '1:1');
else
    legend('ET Klepp vs. ET n(d)', '1:1');
end  
title([cruise_title ' P [mm/h]']);
export_fig([plotdir cruise '_test_odm_R_et_' version_num],'-png','-nocrop','-r300');




if strcmp(cruise,'SPURS2') == 1
R_m_odm = interp1(t_m, R_m, DSD.t, 'linear');

for j = 1:3
    if j == 1
        wh = find(DSD.t < datenum(2017, 1, 1, 0, 0, 0));
        wh_title = 'SPURS-2 2016';
    elseif j == 2
        wh = find(DSD.t > datenum(2017, 1, 1, 0, 0, 0));
        wh_title = 'SPURS-2 2017';
    elseif j == 3
        wh = 1:length(DSD.t);
        wh_title = 'SPURS-2 2016 - 2017';
    end
    figure; 
    hold on;
    plot(k.p(wh), R_m_odm(wh),'.')
    plot(DSD.P(wh), R_m_odm(wh),'o');
    axis square
    plot([0 100],[0 100]);
    xlabel('ODM');ylabel('ship'); 
    xlim([0 100]); ylim(xlim);
    legend('Kyla','ET', '1:1');
    title([wh_title ' P [mm/h] ']);
    export_fig([plotdir cruise '_test_odm_R_ship_' sprintf('%i',j) '_' version_num],'-png','-nocrop','-r300');
end
end

if strcmp(cruise,'TGT') == 1
    figure; hold on;
    plot(DSD.t, DSD.praterain,'v', 'color',rgb('purple'));
    plot(t_k, R_k,'xb');
    plot(DSD.t, DSD.praterain,'oc', 'color',rgb('orange'));
    plot([went_bad_on_TGT went_bad_on_TGT], ylim,'--r','linewidth',2);
    legend('Klepp','Kyla','ET','OMD went bad May 15, 2018','location','best');
    xlim(et_xlim);
    grid on;
    title(cruise_title);
    ylabel('P, mm hr^{-1}');
    datetick('x','mmm-yyyy','keeplimits');
    export_fig([plotdir cruise '_test_odm_R_full_' version_num],'-png','-nocrop','-r300');
end

% % reflfact_k, reflfact
% figure; hold on;
% plot(t_k, reflfact_k, DSD.t, DSD.reflfact);
% legend('Kyla DSD','ET DSD','location','best');
% set(gca,'yscale','log');
% grid on;
% xlim(et_xlim);
% datetick('x','mm/dd','keeplimits');
% title(cruise_title);
% ylabel('z, mm^6 m^-3');
% export_fig([plotdir cruise '_test_odm_z_' version_num],'-png','-nocrop','-r300');

% reflfact_k, reflfact
figure; hold on;
plot(t_k, refl_k,'ob'); 
plot(DSD.t, DSD.Z,'xc');
% plot(DSD.t(DSD.nsumrain > 10), DSD.Z(DSD.nsumrain > 10),'oc');
legend('Kyla DSD','ET n(d)','location','best');
grid on;
xlim(et_xlim);
datetick('x','mm/dd','keeplimits');
title(cruise_title);
ylabel('Z, dBZ');
if strcmp(cruise,'TGT') == 1
    plot([went_bad_on_TGT went_bad_on_TGT], ylim,'--r','linewidth',2);
    legend('Kyla DSD','ET n(d)','ODM went bad May 15, 2018','location','best');
end
export_fig([plotdir cruise '_test_odm_Z_' version_num],'-png','-nocrop','-r300');


figure; 
plot(DSD.refl, DSD.Z,'x'); axis square; hold on; plot([-10 60],[-10 60]);
grid on;
xlabel('Z integral'); ylabel('Z sum');
xlim([-40 80]);
ylim(xlim);
title(cruise_title);
export_fig([plotdir cruise '_test_odm_Z_types_et_' version_num],'-png','-nocrop','-r300');


% lwc
figure; hold on;
plot(DSD.t, DSD.W,'o','color',rgb('green')); 
plot(DSD.t, DSD.lwc,'-c');
plot(DSD.t, 10*DSD.lwc,'p','color',rgb('orange'));
legend('ET mass','ET integral','ET integral*10','location','best');
grid on;
xlim(et_xlim);
datetick('x','mm/dd','keeplimits');
title(cruise_title);
ylabel('W, g m^{-3}');
if strcmp(cruise,'TGT') == 1
    plot(k.t, k.w, 'xb');
    plot([went_bad_on_TGT went_bad_on_TGT], ylim,'--r','linewidth',2);
    legend('ET mass','ET integral','ET integral*10','kyla','ODM went bad May 15, 2018','location','best');
end
export_fig([plotdir cruise '_test_odm_W_' version_num],'-png','-nocrop','-r300');


%% find bad data
if strcmp(cruise,'PISTON') ~= 1
    if strcmp(cruise,'SPURS2') == 1
        R_e = DSD.P;
    elseif strcmp(cruise,'TGT') == 1
        R_e = DSD.P(int_k);
    end
    %%% don't know what these are.... maybe thresholded on drops? n = 4, 48
    bad1 = find(R_e > 100 & R_k < 5);

    %%% don't know what's going on here either
    bad2 = find(R_e < (R_k + 0.01*R_k) & R_k > 12);


    bad = find(R_e./R_k > 2);
    good = find(R_e./R_k < 2);
    R_ratio = R_e./R_k;

    figure; hold on;
    scatter(t_k(good), k.p(good) - R_e(good), 15, R_e(good), 'filled');
    grid on;
    xlim(et_xlim);
    datetick('x','mm/dd','keeplimits');
    title(cruise_title);
    ylabel('P_{kyla} - P_{eliz}');
    c = colorbar;
    title(c,'P [mm/h]');
    if strcmp(cruise,'TGT') == 1
        plot([went_bad_on_TGT went_bad_on_TGT], ylim,'--r','linewidth',2);
        legend('data','ODM went bad May 15, 2018','location','best');
    end
    export_fig([plotdir cruise '_test_odm_R_kyla_et_diff_' version_num],'-png','-nocrop','-r300');

    figure;
    scatter(DSD.Dmax(R_ratio < 2), R_ratio(R_ratio < 2),15,DSD.nsumrain(R_ratio < 2))
    ylabel('P_{eliz} / P_{kyla}');
    xlabel('D_{MAX}, mm');
    grid on;
    c = colorbar;
    title(cruise_title);
    title(c,'# drops');
    export_fig([plotdir cruise '_test_odm_R_kyla_et_ratio_' version_num],'-png','-nocrop','-r300');
end


%% test fall speed parameterizations:
% fortran: Vfall(i)=9.65-10.3*EXP(-1.2*(dcm(i)*10.)/2.)
% matlab: vfall_mm =9.65-10.3*exp(-0.6*(D)); % D in mm

Dcm = DSD.D/10;
vfall_cm        = 9.65-10.3*exp(-1.2*(Dcm*10.)/2.);
vfall_cm2       = 9.65-10.3*exp(-0.6*(Dcm*10.));
vfall_cm3       = 9.65-10.3*exp(-6*Dcm);
vfallcm_cm       = 0.01*(965-1030*exp(-6*Dcm)); % converted to m/s from cm/s equation in Atlas et al. 1973 using data from Gunn and Kinzer 1949

figure; hold on;
plot(Dcm, DSD.vfall,'o','color',rgb('red'),'markersize',15);
plot(Dcm, vfall_cm,'o','color',rgb('orange'),'markersize',12);
plot(Dcm, vfall_cm2,'o','color',rgb('green'),'markersize',9);
plot(Dcm, vfall_cm3,'o','color',rgb('skyblue'),'markersize',6);
plot(Dcm, vfallcm_cm,'o','color',rgb('indigo'),'markersize',3);

ylabel('vfall m/s'); xlabel('D cm')
grid on;

legend('-0.6*(Dmm)','-1.2*(Dcm*10.)/2.','-0.6*(Dcm*10)','-6*Dcm','Atlas et al. 1973','location','best');
title('fall speed unit conversion');
export_fig([syndir 'fall_speeds'],'-png','-nocrop','-r300');