function [rain,D0,mum,lam,N0,Nt,lwc,sigm,Dm]= DSDretrieval(Z,Zdr,model)
%
% Do radar retrieval from Z(dBZ), Zdr(dB) measurements.
% Fix the maximum diameter Dmx=8 mm.
% model=1: exponential DSD
% model=2: C-G DSD 
% Retrieval based on S band radar (10.7 cm wavelength).

if nargin<3
   model=1;
end

%%%% initialize
data1=load('./scattftn80.dat');
d80=data1(:,1);
fh80=data1(:,2);
fv80=data1(:,3);
fd80=data1(:,4);
dd=d80(2)-d80(1);
Dmx=8;  %% fix the Dmx

cc1=10.19;
cc2=0.866;
c0=-0.1021;
c1=4.932;
c2=-0.9551;
c3=0.07934;
c4=-0.002362;
sigma=0.0;
phi=(sigma*pi/180.0);
phi2=phi^2;
epr=80.205+j*17.167;
K2=abs((epr-1)/(epr+2))^2;
wave=107.0;
v=c0+c1*d80+c2*d80.^2+c3*d80.^3+c4*d80.^4; %% terminal vel
%%%%

%%%% iteration initializing
dlm=0.1;
mt=30/dlm;
fp=-1.0;
lam1=1;  % initial value

zdr=10.^(Zdr/10);  %% convert to linear
refl=10.^(Z/10);    %% convert to linear
dbz=Z;
dbzdr=Zdr;
dmx=Dmx;

%%%%% solve DSD parameter lambda from measure Zdr
if ((refl>=1)&(dmx>=0.4)&(zdr>=1))
    if (dbzdr<0.1)  %% use empirical relation
        lwc=refl*10.^(polyval([-0.0411 0.373 -1.439 -2.9694],dbzdr));  % empirical relation of paper 2;
        D0=polyval([0.0226 -0.0785 0.8225 0.5811],dbzdr);  %-99;
        Dm=D0;
        mum=-99;  %-99;
        lam=-99;  %-99;
        N0=-99;  %-99;
        Nt=refl*10.*(polyval([-0.0558 0.5388 -1.8846 0.994],dbzdr));  %-99;
        rain=refl*10.^(polyval([-0.0307 0.2768 -1.1138 -1.9406],dbzdr));   %% @@@@
        sigm=-99;   %% @@@@
    else
        for n3=1:mt
            tf=0;
            lam=0+n3*dlm;
			switch (model)
			case 1
				mum=0; % exponential DSD
			case 2
				mum=-0.0201*lam^2+0.902*lam-1.7178; % C-G DSD 
			end

            zh=0.0;
            zv=0.0;
            nd=d80.^mum.*exp(-lam.*d80);
            nd(d80>Dmx-0.05)=0;  %% here is the same as fortran program   
            temp=((1-2*phi2).*fh80.^2+2*phi2*fh80.*fv80).*nd*dd;
            zh=sum(temp);
            temp=((1-2*phi2).*fv80.^2+2*phi2*fh80.*fv80).*nd*dd;
            zv=sum(temp);
            if (zv==0)
                continue;
            end
            rzdr=zh/zv;
            fc=zdr-rzdr;
            if (fc*fp<0)
                lam=lam1-abs(fp/(fc-fp))*dlm;
                tf=1;
                break;
            else
                fp=fc;
                lam1=lam;
            end
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% calculate rain microphysical characteristics
        if (tf==0)
            N0=nan;  %-99;
            mum=nan;  %-99;
            lam=nan;  %-99;
            D0=nan;  %-99;
            Dm=nan;  %-99;
            sigm=nan;  %-99;
            Nt=nan;  %-99;  %% 2.085*refl.*10.^(0.728*dbzdr.^2-2.066*dbzdr);
            lwc=5.589e-4*refl.*10.^(0.223*dbzdr.^2-1.124*dbzdr);    %% @@@@
            %lwc=2.432e-4*refl.*zdr.^(-4.26);
            rain=7.6e-3*refl.*10.^(0.164*dbzdr.^2-0.897*dbzdr);     %% @@@@
            %rkdp=5.304e-3*refl*zdr.^(-2.3);
            if(dbzdr<=0.35)
                D0=0.717+1.479*dbzdr-0.725*dbzdr.^2+0.171*dbzdr.^3;   %% @@@@
                %D0=0.679+1.606*dbzdr-0.8682*dbzdr.^2;
            end
        else
            
            switch mum
			case 1
				mum=0; % exponential DSD
			case 2
				mum=-0.0201*lam^2+0.902*lam-1.7178; % C-G DSD 
			end

            zh0=0.0;
            tkd=0.0;
            Dm0=0.0;
            D3=0.0;
            D4=0.0;
            rain0=0.0;

            nd=d80.^mum.*exp(-lam.*d80);
            nd(d80>Dmx-0.05)=0;
            zh0=sum(((1-2*phi2).*fh80.^2+2*phi2*fh80.*fv80).*nd*dd);

            D3=sum(nd.*d80.^3*dd);
            D4=sum(nd.*d80.^4*dd);
            rain0=sum(pi/6.0*3.6e-3.*v.*d80.^3.*nd*dd);
            Dm0=sum(nd*dd);

            %%%% estimate the spectrum width       %% @@@@
            Dm=D4/D3;
            sigm2=sum((d80-Dm).^2.*nd.*d80.^3*dd);
            sigm=(sigm2./D3).^(0.5);
            %%%%% 
            fhh2=refl*pi^4*K2/4.0/wave^4;
            N0=fhh2/zh0;
            Nt=Dm0*N0;
            rain=rain0*N0;
            lwc=1e-3*pi/6*N0*D3;        %% @@@@

            %%%% calculate the mean diameter D0
            nd=N0.*d80.^mum.*exp(-lam.*d80);
            nd(d80>(Dmx-0.05))=0;   %% here the bin of Dmx is mm+1, not mm
            mtp1=sum(nd.*d80.^3*dd);
            mtp2=0;
            m=0;
            mas=0;
            masp=0;
            if (mtp1>0)
                while (mas<0.5)
                    masp=mtp2/mtp1;
                    m=m+1;
                    mtp2=mtp2+nd(m).*d80(m).^3*dd;
                    mas=mtp2/mtp1;
                end
                D0=d80(m+1)-(mas-0.5)/(mas-masp)*dd;
            else
                D0=0;
            end

        end
    end
else
    rain=nan;  %-99;
    D0=nan;  %-99;
    Dm=nan;  %-99;
    mum=nan;  %-99;
    lam=nan;  %-99;
    N0=nan;  %-99;
    Nt=nan;  %-99;
    lwc=nan;   %% @@@@
    sigm=nan;   %% @@@@
end
