

function C = RadarEchoClass (Z,Zdr,RhoHV);

% 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

[na,nb]=size(Z);

%%%%
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*75+7.5e-4*75^2;
f16=-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*60+2.92e-3*60^2;
f25=0.68-4.81e-2*50+2.92e-3*50^2;
f26=0.68-4.81e-2*65+2.92e-3*65^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;

XZ=[15  5  5 25  0 20 20  5 40 45;
    20 10 10 30  5 25 25 10 45 50;
    70 20 35 40 20 45 45 45 60 75;
    80 30 40 50 25 55 50 50 65 80];
XZdr=[-4.0 0.0 -0.3 0.5 0.0 -0.3 f21-0.3 f11-0.3 f12-0.3 -0.3;
      -2.0 2.0  0.0 1.0 0.3  0.0 f22     f13     f14      0.0;
       1.0  10  0.3 2.0 3.0  0.3 f31     f23     f24      f15;
       2.0  12  0.6 3.0 3.3  0.6 f32+1   f25+0.5 f26+0.5  f16+0.5];
Xrho=[0.5  0.3 0.95  0.8 0.95  0.9 0.92 0.95 0.92 0.85;
      0.6  0.5 0.98  0.9 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.97 1.01 1.01 1.01 1.01 1.01 1.01];

for n1=1:na
     for n2=1:nb
         for n3=1:10
            tz=Z(n1,n2);
            if (tz<XZ(1,n3))
                PZ(n1,n2,n3)=0;
            elseif (tz<XZ(2,n3))
                PZ(n1,n2,n3)=0.5;
            elseif (tz<XZ(3,n3))
                PZ(n1,n2,n3)=1;
            elseif (tz<XZ(4,n3))
                PZ(n1,n2,n3)=0.5;
            else
                PZ(n1,n2,n3)=0;
            end
            tzdr=Zdr(n1,n2);            
            if (tzdr<XZdr(1,n3))
                PZdr(n1,n2,n3)=0;
            elseif (tzdr<XZdr(2,n3))
                PZdt(n1,n2,n3)=0.5;
            elseif (tzdr<XZdr(3,n3))
                PZdr(n1,n2,n3)=1;
            elseif (tzdr<XZdr(4,n3))
                PZdr(n1,n2,n3)=0.5;
            else
                PZdr(n1,n2,n3)=0;
            end
            trho=RhoHV(n1,n2);
            if (trho<Xrho(1,n3))
                Prho(n1,n2,n3)=0;
            elseif (trho<Xrho(2,n3))
                Prho(n1,n2,n3)=0.5;
            elseif (trho<Xrho(3,n3))
                Prho(n1,n2,n3)=1;
            elseif (trho<Xrho(4,n3))
                Prho(n1,n2,n3)=0.5;
            else
                Prho(n1,n2,n3)=0;
            end
        end
    end
end
Ai=PZ+PZdr+Prho;
[temp,C]=max(Ai,[],3);
C(temp==0)=0;
class=C;

%%%%% correct the classification based on the melting layer
% below the melting layer
Ai1=Ai;
Ai1(:,:,[3,4,5,6])=0;  %only class 1, 2, 7, 8, 9, 10;
% within the melting layer
Ai2=Ai;
Ai2(:,:,[1,5,8,9])=0;  %only class 2, 3, 4, 6, 7, 10;
% above the melting layer
Ai3=Ai;
Ai3(:,:,[2,3,4,7,8,9])=0;  %only class 1, 5, 6, 10;

Ai0=Ai1;  % initiation
for n1=1:na
    x=class(n1,:);

    icetype=find((x==4));  % melting
    icetype=icetype(icetype>300);  
    if length(icetype)>=1
        imax=max(icetype);
        imin=min(icetype);
    else
        imax=9999;
        imin=9999;
    end
    ice=[imin, imax]
    for n2=1:nb
        if n2>imin
            Ai0(n1,n2,:)=Ai2(n1,n2,:);
            if n2>imax
                Ai0(n1,n2,:)=Ai3(n1,n2,:);
            end
        end
    end
end
[temp,C]=max(Ai0,[],3);
C(temp==0)=0;
