
clear all;
close all;

% load example radar data
load ('./KOUN_20050513-083020.mat');

Z=double(RawData.Zh);
Zdr=double(RawData.Zdr);
RhoHV=double(RawData.RhoHV);

na=length(Z(:,1));
ng=length(Z(1,:));

% radar raw data calibration based on disdrometer comparison
Z=Z+1.08;
Zdr=Zdr-0.36;

for model=1:2
    fprintf('\n... 000%%');
    for n1=1:na
        for n2=1:ng
            z=Z(n1,n2);
            zdr=Zdr(n1,n2);            
            [rain,D0,mu,lam,N0,Nt,lwc,sigm,Dm]= DSDretrieval(z,zdr,model);
            retr(model).rain(n1,n2)=rain;
            retr(model).D0(n1,n2)=D0;
            retr(model).mu(n1,n2)=mu;
            retr(model).lam(n1,n2)=lam;
            retr(model).N0(n1,n2)=N0;
            retr(model).Nt(n1,n2)=Nt;
            retr(model).lwc(n1,n2)=lwc;
            retr(model).sigm(n1,n2)=sigm;
            retr(model).Dm(n1,n2)=Dm;
        end
    fprintf('\b\b\b\b\b %3d%%',round(n1/na*100));
    end
end

stop

gate_start=1;
gate_end=length(RawData.Zh(1,:));
Az_deg=RawData.Az;
el_deg=0.5;
resolution=250;
r_min=0;
n1=gate_end;
Az_deg(n1+1)=Az_deg(1);

Az=Az_deg*pi/180;
el_rad=el_deg*pi/180;
range=(([gate_start:gate_end]-1).*resolution+r_min)./1000;

[r,az_rad] = meshgrid(range,Az);
x = r.*cos(el_rad).*sin(az_rad);  
y = r.*cos(el_rad).*cos(az_rad);
z = r.*sin(el_rad);

%%%% plot the figure, Rainrate
figure
set(gcf,'Position',[500 500 1000 450])
set(gcf,'PaperPositionMode','auto')
subplot(1,2,1)
zx=retr(1).rain;
zx(zx<=0)=nan;
zx=log(zx);
zx(n1+1,:)=zx(1,:);
pcolor(x,y,zx)
colormap('Jet')
caxis(log([0.1,100]));
colorbar('ytick',log([0.1,1,3,5,10,15,20,30,50,75,100]),'Yticklabel',{'0.1','1','3','5','10','15','20','30','50','75','100'})
shading flat
axis equal
xlabel('Zonal Distance (km)')
ylabel('Meridional Distance (km)')
axis([-150 150 -150 150])
title('(a) Exponential DSD: Rainrate, {mm} {h^-^1}')

subplot(1,2,2)
zx=retr(2).rain;
zx(zx<=0)=nan;
zx=log(zx);
zx(n1+1,:)=zx(1,:);
pcolor(x,y,zx)
colormap('Jet')
caxis(log([0.1,100]));
colorbar('ytick',log([0.1,1,3,5,10,15,20,30,50,75,100]),'Yticklabel',{'0.1','1','3','5','10','15','20','30','50','75','100'})
shading flat
axis equal
xlabel('Zonal Distance (km)')
ylabel('Meridional Distance (km)')
axis([-150 150 -150 150])
title('(b) C-G DSD: Rainrate, {mm} {h^-^1}')
%%%%

%%%% Dm
figure
set(gcf,'Position',[500 500 1000 450])
set(gcf,'PaperPositionMode','auto')
subplot(1,2,1)
zx=retr(1).Dm;
zx(n1+1,:)=zx(1,:);
pcolor(x,y,zx)
caxis([0 4])
colorbar
shading flat
axis equal
xlabel('Zonal Distance (km)')
ylabel('Meridional Distance (km)')
axis([-150 150 -150 150])
title('(a) Exponential DSD: D_m, {mm}')

subplot(1,2,2)
zx=retr(2).Dm;
zx(n1+1,:)=zx(1,:);
pcolor(x,y,zx)
caxis([0 4])
colorbar
shading flat
axis equal
xlabel('Zonal Distance (km)')
ylabel('Meridional Distance (km)')
axis([-150 150 -150 150])
title('(b) C-G DSD: Rainrate, D_m, {mm}')
%%%%

    
    

    