%% Source: Microwave Radar and Radiometric Remote Sensing, http://mrs.eecs.umich.edu
%% These MATLAB-based computer codes are made available to the remote
%% sensing community with no restrictions. Users may download them and
%% use them as they see fit. The codes are intended as educational tools
%% with limited ranges of applicability, so no guarantees are attached to
%% any of the codes. 
%%
%Code 8.6: Rain Attenuation Coefficient

%Description: Code computes the extinction coefficient (scattering plus
    %absorption) due to rain at any frequency 0<f<1000GHz,
    %rain rate 0<Rr<150mm/hr, and t>0C.

%Input Variables
    %Rr: rain rate in mm/hr
    %f: frequency in GHz
    %t: temperature in degree C
    
%Output Products
    %kapparain: rain attenuation coefficient in dB/km
    
%Book Reference: Sections 8-8.2

%MATLAB Code:

function [kapparain]=RainExtincCoef(Rr,f, t)

lmda = 0.3 ./ f; % wavelength in m

x_max = 2*pi ./ lmda *0.01; %maximum water droplet size 
x_min = 2*pi ./ lmda *1.0e-6;
[epsr epsi] = RelDielConst_PureWater(t,f); %calculate permittivity of water
eps = epsr -1i* epsi;

int_tol = 1.0e-2; % set tolerance for the Quad integration function

kapparain = quad(@(x) calcTermE(x, f, eps, Rr), x_min, x_max, int_tol);

kapparain = kapparain .* lmda.^3/8/ pi^2; % in Np/m

kapparain = kapparain * 4.3429*1000; % in dB/km

end

function [y] = calcTermE(x, f,eps, Rr)

np = length(x);

lmda = 0.3 ./f; %wavelength in meters
rad = x *lmda /(2*pi);

ext_mie= zeros(1, np); p_dist = ext_mie;

No = 8.0e6; 
b = 4100 * Rr^(-0.21);

for n = 1: np
   [Es Ea ext_mie(n) Eb eta_s_r eta_a_r eta_e_r eta_b_r] = ...
    Mie_Rayleigh_ScatteringOfSpheres(rad(n), f, eps, 1); 
    p_dist(n) = No .*exp(-b*2*rad(n)); %marchall palmer's function, normalized and scales!' 
end

y = x.^2 .*p_dist .* ext_mie;
    
end


%% Source: Microwave Radar and Radiometric Remote Sensing, http://mrs.eecs.umich.edu
%% These MATLAB-based computer codes are made available to the remote
%% sensing community with no restrictions. Users may download them and
%% use them as they see fit. The codes are intended as educational tools
%% with limited ranges of applicability, so no guarantees are attached to
%% any of the codes. 
%%
%Code 8.12: Mie and Rayleigh Scattering By Spherical Particle

%Description: Code computes the absorption, scattering, extinction, and
%backscattering efficiencies of a dielectric sphere according to bot the
%Mie solution and the Rayleigh approximation.

%Input Variables
    %r: radius of particle (meters)
    %f: frequency (GHz)
    %epsp: epsp1-j*epsp2 for particles
    %epsb = epsb1 -j epsb2 for background medium 
    
%Output Variables
    %Es: Mie Scattering Efficiency 
    %Ea: Mie Absorption Efficiency
    %Ee: Mie Extinction Efficiency
    %Eb: Mie Backscattering Efficiency
    
    %eta_s_r: Rayleigh Scattering Efficiency 
    %eta_a_r: Rayleigh Absorption Efficiency
    %eta_e_r: Rayleigh Extinction Efficiency
    %eta_b_r: Rayleigh Backscattering Efficiency

%Matlab Code

function [Es Ea Ee Eb eta_s_r eta_a_r eta_e_r eta_b_r] = ... 
    Mie_Rayleigh_ScatteringOfSpheres(r, f, epsp, epsb)

epsb1 = real(epsb);

np = sqrt(epsp); % index of refraction of spherical particle
nb = sqrt(epsb); % index of refraction of background medium

n = np ./nb; % relative index of refraction
%  n = np ./sqrt(epsb1);


chi = 20./3 *pi*r*f*sqrt(epsb1); % normalized circumference in background

%--- Calculate Rayleigh Approximation solution

BigK = (n.^2 - 1) ./(n.^2 + 2); 

eta_s_r = 8/3 .* chi.^4 .* abs(BigK).^2;

eta_a_r = 4 .* chi .* imag(-BigK); 

eta_e_r = eta_a_r + eta_s_r; % Extinction efficiency

eta_b_r = 4 * chi.^4 .* abs(BigK).^2;  % backscattering efficiency



%--- Calculate Mie Scattering solution

%Calculation of Es
l=1;
first = true;
runSum=linspace(0,0,numel(f));
oldSum=linspace(0,0,numel(f));

    %Values of W0 and W-1
    W_1 = sin(chi)+1i*cos(chi);
    W_2 = cos(chi)-1i*sin(chi);
    
    %Value of A0
    A_1 = cot(n*chi);
            
    %A_1=    (sin(real(n).*chi).*cos(real(n).*chi)+j*sinh(imag(n).*chi).*cosh(imag(n).*chi))...
    %        ./(sin(real(n).*chi).^2 + sinh(imag(n).*chi).^2);

while (first || ~endSum(oldSum, runSum, numel(chi)))
    W=(2*l-1)./chi .* W_1 - W_2;

    A = -l./(n.*chi) + (l./(n.*chi)-A_1).^-1;
    
    a = ((A/n + l./chi).*real(W)-real(W_1)) ./ ((A/n+l./chi).*W-W_1);
    b = ((n*A + l./chi).*real(W)-real(W_1)) ./ ((n*A+l./chi).*W-W_1);
    
    sumTerm = (2*l + 1).*(abs(a).^2+abs(b).^2);
    oldSum = runSum;
    runSum = runSum + sumTerm;
    
    %Increment Index varible
    l=l+1;
    
    %Increment W terms
    W_2 = W_1;
    W_1 = W;
    
    %Increment A Terms
    A_1=A;
    
    %Set first pass to false
    first = false;
end
Es = 2./(chi).^2 .* runSum;

%Calculation of Ee
l=1;
first = true;
runSum=linspace(0,0,numel(f));
oldSum=linspace(0,0,numel(f));

    %Values of W0 and W-1
    W_1 = sin(chi)+1i*cos(chi);
    W_2 = cos(chi)-1i*sin(chi);
    
    %Value of A0
    A_1 = cot(n*chi);
            
    %A_1 =   (sin(real(n).*chi).*cos(real(n).*chi)+j*sinh(imag(n).*chi).*cosh(imag(n).*chi))...
    %       ./(sin(real(n).*chi).^2 + sinh(imag(n).*chi).^2);

while (first || ~endSum(oldSum, runSum, numel(chi)))
    W=(2*l-1)./chi .* W_1 - W_2;

    A = -l./(n.*chi) + (l./(n.*chi)-A_1).^-1;
    
    a = ((A/n + l./chi).*real(W)-real(W_1)) ./ ((A/n+l./chi).*W-W_1);
    b = ((n*A + l./chi).*real(W)-real(W_1)) ./ ((n*A+l./chi).*W-W_1);
    
    sumTerm = (2*l + 1).*real(a+b);
    oldSum = runSum;
    runSum = runSum + sumTerm;
    
    %Increment Index varible
    l=l+1;
    
    %Increment W terms
    W_2 = W_1;
    W_1 = W;
    
    %Increment A Terms
    A_1=A;
    
    %Set first pass to false
    first = false;
end
Ee = 2./(chi).^2 .* runSum;

%Calculation of Eb
l=1;
first = true;
runSum=linspace(0,0,numel(f));
oldSum=linspace(0,0,numel(f));

    %Values of W0 and W-1
    W_1 = sin(chi)+1i*cos(chi);
    W_2 = cos(chi)-1i*sin(chi);
    
    %Value of A0
    A_1 = cot(n*chi);
            

while (first || ~endSum(oldSum, runSum, numel(chi)))
    W=(2*l-1)./chi .* W_1 - W_2;

    A = -l./(n.*chi) + (l./(n.*chi)-A_1).^-1;
    
    a = ((A/n + l./chi).*real(W)-real(W_1)) ./ ((A/n+l./chi).*W-W_1);
    b = ((n*A + l./chi).*real(W)-real(W_1)) ./ ((n*A+l./chi).*W-W_1);
    
    sumTerm = (-1)^l .* (2*l + 1) .*(a-b);
    oldSum = runSum;
    runSum = runSum + sumTerm;
    
    %Increment Index varible
    l=l+1;
    
    %Increment W terms
    W_2 = W_1;
    W_1 = W;
    
    %Increment A Terms
    A_1=A;
    
    %Set first pass to false
    first = false;
end

Eb = 1./(chi).^2 .* abs(runSum).^2;


Ea=Ee-Es;

end



function [stop]=endSum(A0, A1, num)
    stop=true;
    pDiff = abs((A1-A0)./A0) .* 100;
    
    for t=1:num,
        if((pDiff(t)>=0.001) || (A0(t) ==0))
            stop=false;
        end
    end

end


%% Source: Microwave Radar and Radiometric Remote Sensing, http://mrs.eecs.umich.edu
%% These MATLAB-based computer codes are made available to the remote
%% sensing community with no restrictions. Users may download them and
%% use them as they see fit. The codes are intended as educational tools
%% with limited ranges of applicability, so no guarantees are attached to
%% any of the codes. 
%%
%Code 4.1: Relative Dielectric Constant of Pure Water
%Description: Code computes the real and imaginary parts of the relative
    %dielectric constant of water at any temperature 0<t<30 and frequency
    %0<f<1000 GHz. Uses the double-Debye model.
%Input Variables:
    %t: temperature in C
    %f: frequency in GHz
%Output Products:
    %epsr: real part of relative dielectric constant
    %epsi: imaginary part of relative dielectric constant
%Book Reference: Section 4.1
%MATLAB Code:  RelDielConst_PureWater.m

%Example call: [A B] = RelDielConst_PureWater(t,f)
%Computes the real and imaginary components of the permitivity of pure
    %water based on the temperature value (t) in degrees C and frequency
    %vector (f) and assigns them to vectors A and B respectively

%MATLAB Code

function [epsr epsi] = RelDielConst_PureWater(t,f)
   
    %Model Paramaters
        a=[0.63000075e1 0.26242021e-2 0.17667420e-3 0.58366888e3 0.12634992e3 ...
            0.69227972e-4 0.30742330e3 0.12634992e3 0.37245044e1 0.92609781e-2]; %Eqn parameters

        epsS = 87.85306*exp(-0.00456992*t);
        epsOne = a(1)*exp(-a(2)*t);
        tau1 = a(3)*exp(a(4)/(t+a(5)));
        tau2 = a(6)*exp(a(7)/(t+a(8)));
        epsInf = a(9) + a(10)*t;

        %Permitivity Calculation
        eps = ((epsS-epsOne)./(1-j*2*pi.*f.*tau1)) + ((epsOne-epsInf)./(1-j*2*pi.*f.*tau2)) + (epsInf);

        %Extracts the real and imaginary components of relative permitivity
        epsr = real(eps);
        epsi = imag(eps);
  end

