function [Z,Zdr,R,D0,Adp,Kdp,Rhohv,LWC,Re,Rc,Dm,Sigm,Ahh,Nt]= calc_variable(D,Nd)

% Given DSD, calculate the C-band parameters (R,D0,Z,Zdr,Adp,Kdp,Rhohv,...)
% this is based on S band (10 cm)
% This code is for 41 bins [0.1,0.3,0.5,...8.1], D should be at least 0.1 mm spaced
% Unit: Z(dBZ), Zdr(dB), R(mm/h), D0(mm), Adp(dB/km), Kdp(deg/km)...
% Rhohv, LWC(g/m^3), Re(g/m^3/s), Rc(g/m^3/s), Dm(mm), Sigm(mm), Ahh(dB/km)

% The scattering amplitude is for S-band 10 cm
fort18=load('./Sband_Scatter.dat');

% fort18 is the scattering amplitude data for D = [0.1:0.1:10]; 
% here D should include bins which space 0.2 mm. i.e. [0.1,0.3,0.5,...];
% so the interpolation is performed
D=D(:);
Nd=Nd(:);
D=[-1;D];
Nd=[0;Nd];  % enlarge the range to do the interpolation

d=fort18(:,1);
nd=interp1(D,Nd,d);  % interpolate 
d=d(:);
nd=nd(:);
len=length(nd);
dd=d(2)-d(1);

% backward scattering
bar=fort18(:,2);   % *1000 is to convert to mm
bai=fort18(:,3);
bbr=fort18(:,4);
bbi=fort18(:,5);
% forward scattering
far=fort18(:,6);
fai=fort18(:,7);
fbr=fort18(:,8);
fbi=fort18(:,9);
%%%%
bam=abs(bar+i.*bai);
bbm=abs(bbr+i.*bbi);
fam=abs(far+i.*fai);
fbm=abs(fbr+i.*fbi);
faa=angle(fbr+i.*fai).*180./pi;
fba=angle(fbr+i.*fbi).*180./pi;
baa=angle(bar+i.*bai).*180./pi;
bba=angle(-bbr-i.*bbi).*180./pi;

c0=-0.1021; c1=4.932; c2=-0.9551; c3=0.07934; c4=-0.002362; 
v=c0+c1.*d+c2.*d.^2+c3.*d.^3+c4.*d.^4;    %% terminal velocity of Ed. 
fh=bam;
fv=bbm;
fhc=bar+i.*bai;
fvc=bbr+i.*bbi;
fd=far-fbr;
ah=fai;
av=fbi;
%%%%
wave=100;  %% mm
epr=79.669+j*18.2257;    %% for S-band 10 cm
phi=0*pi/180;
K2=abs((epr-1)/(epr+2))^2;
fhh2=pi^4*K2/4/wave^4;
kx=2*pi/wave;
%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% calculate rain microphysical characteristics
Fh2=sum(fh.^2.*nd).*dd;
Fv2=sum(fv.^2.*nd).*dd;
Fhv=sum(fhc.*fvc.*nd).*dd;
Fd=sum(fd.*nd).*dd;
Ah=sum(ah.*nd).*dd;
Av=sum(av.*nd).*dd;
ath=(1-2*phi^2)*Ah+2*phi^2*sqrt(Ah*Av);
atv=(1-2*phi^2)*Av+2*phi^2*sqrt(Ah*Av);

temp1=((1-2*phi^2)*fh.^2+2*phi^2.*fhc.*fvc).*nd*dd; 
temp2=((1-2*phi^2)*fv.^2+2*phi^2.*fhc.*fvc).*nd*dd;
temp3=wave*(1-2*phi^2)*fd.*nd*dd*180/pi/1e+3;
zh=abs(sum(temp1,1)')./fhh2; 
zv=abs(sum(temp2,1)')./fhh2;  

Z=10*log10(zh);
Zdr=10*log10(zh./zv);
Kdp=sum(temp3,1)';
Rhohv=abs(Fhv)./sqrt(Fh2.*Fv2);  
Adp=4.343e-3*2*wave*(ath-atv);    %% unit (dB/km)
Ahh=4.343e-3*2*wave*ath;    %% unit (dB/km)

%%%%
D3=sum(nd.*d.^3*dd);   % moment 3
D4=sum(nd.*d.^4*dd);
R=sum(pi/6.0*3.6e-3.*v.*d.^3.*nd*dd);
LWC=1e-3*pi/6.*D3;
%%%% estimate the spectrum width     
Nt=sum(nd.*dd);   % moment 0
Dm=D4/D3;
sigm2=sum((d-Dm).^2.*nd.*d.^3*dd);
Sigm=(sigm2./D3).^(0.5);
%%%% calculate the mean diameter D0
D0=D0find(d,nd);
%%%%
%%%% calculate the Re, Rc, Vtm
temp=3.55e-7.*nd.*d.^(8/5).*dd;
Re=sum(temp,1)';

temp=pi/4*(1e-6).*nd.*(d.^2).*v.*dd;  %% use fitted Vt
Rc=sum(temp,1)';    %% based on fitted Vt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



function D0 = D0find (D,Nd)

% D0 = D0find (D,Nd); calculate the mean volume diameter D0
% assume dd=0.2;    
    
dd=0.2;
mtemp=0;
m3=Nd.*D.^3.*dd;
mtemp1=sum(m3);
m=0;
mas=0;
if (mtemp1>0)
    while (mas<0.5)
        masp=mtemp/mtemp1;
        m=m+1;
        mtemp=mtemp+m3(m);
        mas=mtemp/mtemp1;
    end
    D0=D(m)-(mas-0.5)/(mas-masp)*dd;
else
    D0=0;
end

