
clear all;
close all;

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

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

% Use 3 inputs to do the classification 
% the melting layer is located at 2.6-3.1 km
[class, class_t] = RadarEchoClass_3input(Z,Zdr,RhoHV,[2.6,3.1]); 

%%%% plot HydroClass
figure
set(gcf,'Position',[500 500 600 450])
set(gcf,'PaperPositionMode','auto')
Zx=class;

rmax=150;
resolution=250;
Az_deg=RawData.Az;
el_deg=0.5;
gate_start=1;
gate_end=length(Zx(1,:));
r_min=0;

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);

pcolor(x,y,Zx)
colormap ([1    1    1; ...
           0.20 1.00 1.00; ...
           0.00 0.00 0.80; ...
           0.00 0.80 0.00; ...
           0.00 0.60 0.00; ...
           1.00 1.00 0.20; ...
           1.00 0.00 0.00; ...
           0.80 0.00 0.00; ...
           1.00 0.00 0.80; ...
           0.60 0.20 0.80; ...
           0.50 0.00 0.50]);
caxis([0 10]);
colorbar('ytick',[0.5:9/10:9.5],'Yticklabel',{'NE','GC/AP','BS','DS','WS','CR','GR','BD','RA','HR','HA'})
shading flat
axis([-rmax rmax -rmax rmax])
xlabel('Zonal Distance (km)')
ylabel('Meridional Distance (km)')
title(' Radar Echo Class')
