

function [C, T] = RadarEchoClass_3input (Z,Zdr,RhoHV,Hbt,resolution,el_deg)

% This applies only 3 inputs of Park et al. 2009
% Hbt=[Hb,Ht];  % bottom/top of melting layer, % km
% Rbt=[Rbb,Rb,Rt,Rtt];  % km
% resolution: in meter
% Input Z,Zdr,RhoHV are vectors or 2-d matrice with the same order

if (nargin<3)
    error('Input parameters are less than 3')
end
if (nargin<4)
    Hbt=[-99 -99];
    %Rbt=[-99 -99 -99 -99];
end
if (nargin<5)
    resolution=250;
    el_deg=0.5;
end

Rbt = range_ML(Hbt,el_deg);

[na,nb]=size(Z);
class=zeros(na,nb);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% initialization
f11=-0.5+2.5e-3*5+7.5e-4*5^2;
f12=-0.5+2.5e-3*40+7.5e-4*40^2;
f13=-0.5+2.5e-3*10+7.5e-4*10^2;
f14=-0.5+2.5e-3*45+7.5e-4*45^2;
f15=-0.5+2.5e-3*50+7.5e-4*50^2;
f16=-0.5+2.5e-3*75+7.5e-4*75^2;
f17=-0.5+2.5e-3*55+7.5e-4*55^2;
f18=-0.5+2.5e-3*80+7.5e-4*80^2;

f21=0.68-4.81e-2*20+2.92e-3*20^2;
f22=0.68-4.81e-2*25+2.92e-3*25^2;
f23=0.68-4.81e-2*45+2.92e-3*45^2;
f24=0.68-4.81e-2*55+2.92e-3*55^2;
f25=0.68-4.81e-2*50+2.92e-3*50^2;
f26=0.68-4.81e-2*60+2.92e-3*60^2;

f31=1.42+6.67e-2*45+4.85e-4*45^2;
f32=1.42+6.67e-2*50+4.85e-4*50^2;

% (1)GC/AP (2)BS (3)DS (4)WS (5)CR (6)GR (7)BD (8)RA (9)HR (10)RH
A.XZ=[15  5  5 25  0 25 20  5 40 45;
    20 10 10 30  5 35 25 10 45 50;
    70 20 35 40 20 50 45 45 55 75;
    80 30 40 50 25 55 50 50 60 80];
A.XZdr=[-4.0 0.0 -0.3 0.5 0.1 -0.3 f21-0.3 f11-0.3 f12-0.3 -0.3;
      -2.0 2.0  0.0 1.0 0.4  0.0  f22     f13     f14      0.0;
       1.0  10  0.3 2.0 3.0  f15  f31     f23     f24      f16;
       2.0  12  0.6 3.0 3.3  f17+0.3  f32+1   f25+0.5 f26+0.5  f18+0.5];
A.Xrho=[0.5  0.3 0.95  0.88 0.95  0.9 0.92 0.95 0.92 0.85;
      0.6  0.5 0.98  0.92 0.98 0.97 0.95 0.97 0.95  0.9;
      0.9  0.8  1.0 0.95  1.0  1.0  1.0  1.0  1.0  1.0;
     0.95 0.83 1.01 0.985 1.01 1.01 1.01 1.01 1.01 1.01];

A.WZ =  [0.2 0.4 1.0 0.6 1.0 0.8 0.8 1.0 1.0 1.0]';
A.WZdr=[0.4 0.6 0.8 0.8 0.6 1.0 1.0 0.8 0.8 0.8]';
A.Wrho=[1.0 1.0 0.6 1.0 0.4 0.4 0.6 0.6 0.6 0.6]';

A.Rbb=Rbt(1);
A.Rb=Rbt(2);
A.Rt=Rbt(3);
A.Rtt=Rbt(4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

TP=['NE';'GC'; 'BS'; 'DS'; 'WS'; 'CR'; 'GR'; 'BD'; 'RA'; 'HR'; 'RH'];
for n1=1:na
    for n2=1:nb
        z=Z(n1,n2);
        zdr=Zdr(n1,n2);
        rhohv=RhoHV(n1,n2);
        A.R=n2*resolution/1000;   % km
        tmp = classification_pixel (z,zdr,rhohv,A);
        class(n1,n2)=tmp;
        T{n1,n2}=TP(tmp+1,:);
    end
end
C=class;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% subroutine %%%%%%%%%%%%%%%%%%%%%%% 
function class = classification_pixel (z,zdr,rhohv,A)

XZ=A.XZ;
XZdr=A.XZdr;
Xrho=A.Xrho;
WZ=A.WZ;
WZdr=A.WZdr;
Wrho=A.Wrho;
R=A.R;
Rbb=A.Rbb;
Rb=A.Rb;
Rt=A.Rt;
Rtt=A.Rtt;


PZ=zeros(10,1);
PZdr=zeros(10,1);
Prho=zeros(10,1);
for n3=1:10
    tz=z;
    if (tz<XZ(1,n3))
        PZ(n3,1)=0;
    elseif (tz<XZ(2,n3))
        PZ(n3,1)=(tz-XZ(1,n3))./(XZ(2,n3)-XZ(1,n3));
    elseif (tz<XZ(3,n3))
        PZ(n3,1)=1;
    elseif (tz<XZ(4,n3))
        PZ(n3,1)=(XZ(4,n3)-tz)./(XZ(4,n3)-XZ(3,n3));
    else
        PZ(n3,1)=0;
    end

    tz=zdr;
    if (tz<XZdr(1,n3))
        PZdr(n3,1)=0;
    elseif (tz<XZdr(2,n3))
        PZdr(n3,1)=(tz-XZdr(1,n3))./(XZdr(2,n3)-XZdr(1,n3));
    elseif (tz<XZdr(3,n3))
        PZdr(n3,1)=1;
    elseif (tz<XZdr(4,n3))
        PZdr(n3,1)=(XZdr(4,n3)-tz)./(XZdr(4,n3)-XZdr(3,n3));
    else
        PZdr(n3,1)=0;
    end
    
    tz=rhohv;
    if (tz<Xrho(1,n3))
        Prho(n3,1)=0;
    elseif (tz<Xrho(2,n3))
        Prho(n3,1)=(tz-Xrho(1,n3))./(Xrho(2,n3)-Xrho(1,n3));
    elseif (tz<Xrho(3,n3))
        Prho(n3,1)=1;
    elseif (tz<Xrho(4,n3))
        Prho(n3,1)=(Xrho(4,n3)-tz)./(Xrho(4,n3)-Xrho(3,n3));
    else
        Prho(n3,1)=0;
    end
end

Ai=(PZ.*WZ+PZdr.*WZdr+Prho.*Wrho)./sum(WZ+WZdr+Wrho);

if (Rbb<0)
    % no melting layer information
else
    % correct the classification based on the melting layer
    % (1)GC/AP (2)BS (3)DS (4)WS (5)CR (6)GR (7)BD (8)RA (9)HR (10)RH
    
    %0<R<Rbb ?GC/AP, BS, BD, RA, HR, RH;
    if (R<Rbb)
       Ai([3,4,5,6])=0;  %only class 1, 2, 7, 8, 9, 10;
    end

    %Rbb<R<Rb ?GC/AP, BS, WS, GR, BD, RA, HR, RH
    if (Rbb<R)&&(R<Rb)
       Ai([3,5])=0;  
    end 
    
    %Rb<R<Rt ?GC/AP, BS, DS, WS, GR, BD, RH;
    if (Rb<R)&&(R<Rt)
       Ai([5,8,9])=0; 
    end 
    
    %Rt<R<Rtt ?GC/AP, BS, DS, WS, CR, GR, BD, RH; 
    if (Rt<R)&&(R<Rtt)
       Ai([8,9])=0; 
    end 
    
    %R>Rtt ?DS, CR, GR, RH.
    if (R>Rtt)
       Ai([1,2,4,7,8,9])=0; 
    end 

end

[temp,C]=max(Ai);
if (temp<=0) 
    C=0; 
end
class=C;
    


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rbt = range_ML(Hbt,el_deg)

% Hbt=[Hb,Ht];  % bottom/top of melting layer, % km
% Rbt=[Rbb,Rb,Rt,Rtt];  % km

% a=6374 km, earth radius
% use ae=4/3*a rule 
% h=ae*(cos(el)/cos(el+s/ae)-1)

if Hbt(1)>0
    
    a=6374; % km, earth radius
    ae=4/3*a; % 4/3 rule 

    el1=(el_deg+0.5)/180*pi;
    el2=el_deg/180*pi;
    el3=(el_deg-0.5)/180*pi;

    Hb=Hbt(1);
    Ht=Hbt(2);

    Rbb=(acos(cos(el1)/(Hb/ae+1))-el1)*ae;
    Rb=(acos(cos(el2)/(Hb/ae+1))-el2)*ae;
    Rt=(acos(cos(el2)/(Ht/ae+1))-el2)*ae;
    Rtt=(acos(cos(el3)/(Ht/ae+1))-el3)*ae;

    Rbt=[Rbb,Rb,Rt,Rtt];

else
    
    Rbt=[-99 -99 -99 -99];

end












