
clear all;
close all;

load('./dsddata_20050513.mat')

%%%%
len=length(dsd_data(1,1,:));
dd=0.2;

Ct_all=squeeze(dsd_data(:,2,:));   %% Total # for every bin
D_all=squeeze(dsd_data(:,3,:));    %% Diameter: bin spaced by D 
D1_all=squeeze(dsd_data(:,4,:));   %% Arithmatical mean Diameter 
D2_all=squeeze(dsd_data(:,5,:));   %% Median Volume Diameter
Nd_all=squeeze(dsd_data(:,6,:));   %% # DSD

%%%% calculate Tct
Tct=sum(Ct_all,1)';
%%%% calculate Max diameter
Dx=max(D1_all,[],1)';
%%%%

for n=1:len
    %%%% time
    disd.day(n,1)=strread(files(n,1:6));  % year:month:day, hour:minute
    disd.hour(n,1)=strread(files(n,8:9));
    disd.minute(n,1)=strread(files(n,10:11));
    disd.time(n,1)=disd.hour(n,1)+disd.minute(n,1)/60;
end

disd.Tct=Tct;
disd.Dmax=Dx;

for n=1:len
    d=D_all(:,n);
    d1=D1_all(:,n);
    d(d1>0)=d1(d1>0);
    nd=Nd_all(:,n);
    [Z,Zdr,R,D0,Adp,Kdp,Rhohv,LWC,Re,Rc,Dm,Sigm,Ahh,Nt]= calc_variable(d,nd);
    disd.Z(n,1)=Z;
    disd.Zdr(n,1)=Zdr;
    disd.Kdp(n,1)=Kdp;
    disd.Rhohv(n,1)=Rhohv;
    disd.R(n,1)=R;
    disd.D0(n,1)=D0;
    disd.Ah(n,1)=Ahh;
    disd.Adp(n,1)=Adp;
    disd.W(n,1)=LWC;
    disd.Dm(n,1)=Dm;
    disd.Sigm(n,1)=Sigm;
    disd.Re(n,1)=Re;
    disd.Rc(n,1)=Rc;
    disd.Nt(n,1)=Nt;
end
disd.DSD=dsd_data;

%%%%%%%%%%%%%%%%%%%% figures
tim=disd.time;
R=disd.R;
figure(2)
set(gcf,'Position',[500 500 500 800])
set(gcf,'PaperPositionMode','auto')
subplot(3,1,1)
plot(tim,R,'+')
xlabel('Time, UTC')
ylabel('R, mm h^{-1}')
title('(a)')

subplot(3,1,2)
Nt=disd.Nt;
plot(tim,Nt,'+')
xlabel('Time, UTC')
ylabel('N_T, # m^{-3}')
title('(b)')

subplot(3,1,3)
Dm=disd.Dm;
plot(tim,Dm,'+')
ylim([0 6])
xlabel('Time, UTC')
ylabel('D_m, mm')
title('(c)')

%%%% radar variables
figure(3)
set(gcf,'Position',[500 500 500 800])
set(gcf,'PaperPositionMode','auto')
Z=disd.Z;
subplot(3,1,1)
plot(tim,Z,'+')
legend(['0',num2str(disd.day(1))])
xlabel('Time, UTC')
ylabel('Z_H, dBZ')

subplot(3,1,2)
Zdr=disd.Zdr;
plot(tim,Zdr,'+')
xlabel('Time, UTC')
ylabel('Z_{DR}, dB')

subplot(3,1,3)
RhoHV=disd.Rhohv;
plot(tim,RhoHV,'+')
ylim([0.98 1])
xlabel('Time, UTC')
ylabel('\rho_{hv}')


%%%% radar variables
figure(4)
set(gcf,'Position',[500 500 500 800])
set(gcf,'PaperPositionMode','auto')
subplot(3,1,1)
Kdp=disd.Kdp;
plot(tim,Kdp,'+')
ylim([0 2])
xlabel('Time, UTC')
ylabel('K_{DP}, ^\circkm^{-1}')

subplot(3,1,2)
Ah=disd.Ah;
plot(tim,Ah,'+')
xlabel('Time, UTC')
ylabel('A_{H}, dB km^{-1}')

subplot(3,1,3)
Adp=disd.Adp;
plot(tim,Adp,'+')
xlabel('Time, UTC')
ylabel('A_{DP}, dB km^{-1}')


