%% 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.10: Zenith Atmospheric Opacity - Cloudy Atmosphere

%Description: Using a Standard Atmosphere model, code computes
    %opacity of the entire atmosphere when containing a cloud layer
    % extending between heights z3 and z4, along the zenith
    %from sea-surface height z1=0 to z2=32km, at any
    %frequency 0<f<50 GHz.

%Input Variables
    %T0: sea-surface temperatre in 0C
    %P0: sea-surface total barometric pressure in mbar
    %rhov0: sea-surface water vapor density in g/m^3
    %f: frequency in GHz
    %z3: base height of cloud layer in km
    %z4: top height of cloud layer in km
    %cloud type: water or ice
    %mv: cloud water content in g/m^3
    
%Output Products:
    %Tau0: zenith opacity in Np
    % Tr: Transmissivity
    
%Book References: Sections 8-3 and 8-7

%MATLAB Code:


function [tau0 Tr]=ZenithAtmosOpacity_CloudyAtmos(T0,P0,rhov0,f,z3,z4,type,mv)
z1=0;
z2=32;
theta=0;

[tau0 Tr]= OpticalThickn_CloudyAtmos(T0, P0, rhov0, f, z1, z2, z3, z4, theta, type, mv);

end


%Code 8.8: Optical Thickness - Cloudy Atmosphere

%Description: Code computes the optical thickness under cloudy
    %conditions for an atmosphereic layer between heights z1 and z2,
    %containing a cloud between z3 and z4 for propagation along
    %direction theta relative to the zenith (0<theta<80) at any frquency
    %0<f<1000 GHz, using a Standard Atmospheric Model.  It is assumed that
    %the cloud layer is partially or totally within the atmospheric layer
    %of interest (z1 to z2).

%Input Variables
    %T0: Sea-surface temperatre in  degrees C
    %P0: sea-surface total barometric pressure in mbar
    %rhov0: sea-surface water vapor density in g/m^3
    %f: frequency in GHz
    %z1: base height of atmospheric layer in km
    %z2: top height of atmospheric layer in km
    %z3: base height of the cloud layer in km
    %z4: top height of the cloud layer in km
    %theta: zenith angle in degrees
    %cloud type: "water" or "ice"
    %mv: cloud waer content in g/m^3
    
%Output Products
    %tau: optical thickness in Np
    %Tr: Transmissivity (unitless)
%Book Reference: Sections 7-7

%MATLAB Code:

function [tau Tr]=OpticalThickn_CloudyAtmos(T0, P0, rhov0, f, z1, z2, z3, z4, theta, type, mv)  

%-- calculate the optical thickness of clear atmosphere 

  [tau_g Tr] = OpticalThickn_ClearAtmos(T0,P0,rhov0,f,z1,z2,theta); % 

%-- calculate the optical thickness of the cloud region
  
    int_tol = 1.0e-4; % set tolerance for the Quad integration function

    t = quad(@(z) calcCloud(z, T0,f, type, mv), z3, z4, int_tol);

    tau_c = t.*sec(deg2rad(theta));
    
    tau = tau_c + tau_g; % total optical thickness
    Tr = exp(-tau);
end

function [term]=calcCloud(z, T0, f,type,mv)

np = length(z);
term = zeros(np,1);
T0 = T0 + 273; % Transform to Kelvin

for ii=1:np
    %Calculate T
        if((z(ii)>=0)&&(z(ii)<11))
            T =T0-6.5*z(ii);
        elseif((z(ii)>=11)&&(z(ii)<20))
            T = T0-6.5*11;
        elseif((z(ii)>=20)&&(z(ii)<=32))
            T = (T0-6.5*11) +(z(ii)-20);
        end
        
        T = T - 273; %back to Celcius
        
    if (strcmp(type,'water'))    
        term(ii) = ExtincCoef_WaterCloud(T,mv,f); % in dB/km
    elseif (strcmp(type,'ice'))
        term(ii) = ExtincCoef_IceCloud(T,mv,f); % in dB/km
    else
        error('Invalid Cloud Type')
    end
    term(ii) = term(ii) / 4.34; % Transform to Np
    
end
end


%Code 8.7: Optical Thickness - Clear Atmosphere

%Description: Code computes the optical thickness under clear-sky
    %conditions for an atmosphereic layer between heights z1 and z2 for
    %propagation along direction theta relative to the zenith (0<theta<80)
    %at any ffrquency 0<f<1000 GHz, using a Standard Atmospheric Model.

%Input Variables
    %T0: Sea-surface temperature in degrees C
    %P0: sea-surface total barometric pressure in mbar
    %rhov0: sea-surface water vapor density in g/m^3
    %f: frequency in GHz
    %z1: base height of atmospheric layer in km
    %z2: top height of atmospheric layer in km
    %theta: zenith angle in degrees
    
%Output Products
    %tau: optical thickness in Np
    %Tr: Transmissivity (unitness)
%Book Reference: Sections 8-3

%MATLAB Code:

function [tau Tr]=OpticalThickn_ClearAtmos(T0, P0, rhov0, f, z1, z2, theta) 

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

    t = quad(@(z) calcTerm(z, T0,P0,rhov0,f), z1, z2, int_tol);
    tau = t.*sec(deg2rad(theta)); % optical thickness in Np
    Tr = exp(-tau ); 
end
%
function [term]=calcTerm(z, T0,P0,rhov0,f)

    T0 = T0+273; % transform the temperature to Kelvin
        
    np = length(z); % size of vector z
    term = zeros(np,1);
        
    for ii = 1:np

    %Calculate T
        if((z(ii)>=0)&&(z(ii)<11))
            T=T0-6.5*z(ii);
        elseif((z(ii)>=11)&&(z(ii)<20))
            T = T0-6.5*11;
        elseif((z(ii)>=20)&&(z(ii)<=32))
            T = (T0-6.5*11) +(z(ii)-20);
        end
    
    %Calculate P
        P = P0*exp(-z(ii)/7.7);
    
    %Calculate rhov
        rhov = rhov0*exp(-z(ii)/2.1);
    
        T = T - 273; % transform back to C
        term(ii) = AttenCoef_Gaseous(T,P,rhov,f);
        term(ii) = term(ii) /4.34; % transform the atten to Np
        
%        z(ii), term(ii)
    end
        
end


%Code 8.3: Gaseous Attenuation Coefficient

%Description: Code computes the total attenuation (absorption) coefficient
    %due to Oxygen and water vapor at any frequency 1<f<1000GHz,
    %temperature -100<t<50C, pressure 10^-5 mbar < P < 1013 mbar,
    %and water vapor density 0<rho_0<20g/m^3.

%Input Variables
    %T: temperature in in degree C
    %P: total barometric pressure in mbar
    %rho_0: water vapor density in g/m^3
    %f: frequency in GHz
    
%Output Products
    %kappag: gaseous attenuation coefficient in dB/km
    
%Book Reference: Section 8-2.3

%MATLAB Code:

function [kappag]=AttenCoef_Gaseous(T,P,rho_0,f)

%-- calculate atten. due to water vapor
   [kappaH2O]=AttenCoef_WaterVapor(T,P,rho_0,f);

% calculate atten. due to oxygen    
   [kappaO2]=AttenCoef_Oxygen(T,P,rho_0,f);
   
% total atten. in dB/km   
   kappag=kappaH2O+kappaO2;
end


%Code 8.1: Attenuation Coefficient of Water Vapor

%Description: Code computes the attenuation (absorption) coefficient of
%water vapor at any frequency 1<f<1000GHz, temperature -100<t<50C
%pressure 10e-5mbar<P<1013mbar, and water vapor density
%0<rho_0<20<g/m^3
%Input Variables:
    %T: temperature in Celsius
    %P: total barometric pressure in mbar
    %rho_0: water vapor density in g/m^3
    %f: frequency in GHz
%Output Products:
    %kappaH20: water vapor attenuation coefficient in dB/km

%Book Reference: Sections 8-2.3

%MATLAB Code:

function [kappaH2O]=AttenCoef_WaterVapor(T,P,rho_0,f)

    T = T + 273.15; % transform temperature to degrees Kelvin
    
    F=zeros(34,numel(f));
    term=zeros(34,numel(f));

    %Load Parameters from Parameter List
    load WATER.txt
    v=transpose(WATER(1:34,1));
    b([1:6],:)=transpose(WATER(1:34,[2:7]));

    %Calculate values of computed inputs
    theta=300/T;
    e=rho_0/(0.7223*theta);
    Pd= P - e;
    
    ks=0.444*theta^(7.5);
    kf=0.0145*theta^(4.5);
    Nc2 = e*(ks*e+kf*Pd)*10^-6.*f;

    delta = 0;
    
    gamma = b(3,:)*10^-3 .*(b(4,:).*e.*theta.^b(6,:) + Pd.*theta.^b(5,:));
    
    for ii=1:34
        F(ii,:)=f.*((1 -1i*delta)./(v(ii) -f -1i*gamma(ii))-(1+1i*delta)...
            ./(v(ii) + f +1i*gamma(ii)));
    end
    
    Sl = (b(1,:)./v).*e.*theta^(3.5).*exp(b(2,:)*(1-theta));

    for ii=1:34
        term(ii,:) = Sl(ii).*F(ii,:);
    end
    
    Nv2 = imag(sum(term))+Nc2;

    kappaH2O = 0.182*f.*Nv2;
end


%Code 8.2: Attenuation Coefficient of Oxygen

%Description: Code computes the attenuation (absorption) coefficient of
%Oxygen at any frequency 1<f<1000GHz, temperature -100<t<50C
%pressure 10e-5mbar<P<10e13mbar, and water vapor density
%0<rho_0<20<g/m^3

%Input Variables:
    %T: temperature in Celsius
    %P: total barometric pressure in mbar
    %rho_0: water vapor density in g/m^3
    %f: frequency in GHz
%Output Products:
    %kappa02: Oxygen attenuation coefficient in dB/km

%Book Reference: Sections 8-2.3

%Description: Attenuation Coefficient of Oxygen
function [kappaO2]= AttenCoef_Oxygen(T,P,rho_0,f)

    T = T + 273.15; % transform temperature to degrees Kelvin
    
    F=zeros(44,numel(f));
    term=zeros(44,numel(f));    

    %Load Parameters from Parameter List
    load OXYGEN.txt
    v=transpose(OXYGEN(1:44,1));
    a(1:6,:)=transpose(OXYGEN(1:44,[2:7]));
    a(1,:)=a(1,:)*10^-6;

    %Calculate values of computed inputs
    theta=300/T;
    e=rho_0/(0.7223*theta);
    Pd=P-e;
    
    gamma0 = 0.56e-3*P*theta^0.8;
    S0 = 6.14e-5*Pd*theta^2;
    F0 = -f./(f+ 1i*gamma0);
    
    Sn = 1.4e-12*Pd^2*theta^3.5;
    Fn2 = f./(1 + 1.9e-5*f.^1.5);
    
    Nn = S0.*F0 + 1i*Sn*Fn2;

    delta = 0;
%     delta = (a(5,:)+ a(6,:) .* theta) .*P .* theta.^(0.8);
    
    gamma = a(3,:)*10^-3 .*(Pd.*theta.^a(4,:) + 1.10*e*theta);
    
    for ii=1:44
        F(ii,:)=f.*((1-1i*delta)./(v(ii)-f-1i*gamma(ii))-(1+1i*delta)...
            ./(v(ii)+f+ 1i*gamma(ii)));
    end
    
    Sk = (a(1,:)./v).*Pd.*theta^(3).*exp(a(2,:)*(1-theta));

    for ii=1:44
        term(ii,:) = Sk(ii).*F(ii,:);
    end
    Nd2 = imag(sum(term))+imag(Nn);

    kappaO2 = 0.182*f.*Nd2;
end



%Code 8.4: Water Cloud Attenuation Coefficient

%Description: Code computes the extinction coefficient (scattering plus
    %absorption) for a water cloud at any frequency 0 < f < 50 GHz,
    %temperature 0 < t < 50 C, and cloud water content 0 < mv <5 g/m^3.

%Input Variables
    %T: temperature in degree C
    %mv: cloud water content in g/m^3
    %f: frequency in GHz
    
%Output Products
    %kappawc: cloud extinction coefficient in dB/km
    
%Book Reference: Section 8-7.3

%MATLAB Code:

function [kappawc]=ExtincCoef_WaterCloud(T,mv,f)

    lmda_0 = 30 ./f; % wavelength in cm
    
    %Dielecric Constant of pure water
    
    [epsr epsi] = RelDielConst_PureWater(T,f);
    
    epsw = epsr- 1i*epsi;
    
    k= imag(-(epsw-1)./(epsw+2));

    kappawc = 0.434 * 6*pi/lmda_0 * mv * k;
end


%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


%Code 8.5: Ice Cloud Attenuation Coefficient

%Description: Code computes the extinction coefficient (scattering plus
    %absorption) for an ice cloud at any frequency 0 < f < 70 GHz,
    %temperature -50 < t < 0 C, and coud water content 0 < mv < 5g/m^3.

%Input Variables
    %T: temperature in degree C
    %mv: cloud water content in g/m^3
    %f: frequency in GHz
    
%Output Products
    %kappaic: cloud attenuation coefficient in dB/km
    
%Book Reference: Section 8-7.3

%MATLAB Code:

function [kappaic]= ExtincCoef_IceCloud(T,mv,f)

    lmda_0 = 30 ./f; % wavelength in cm

    %Dielecric Constant of pure ice
    [epsr epsi]=RelDielConst_PureIce(T,f);
    
    epsic = epsr - 1i*epsi;
    Nic2 = imag( -(epsic-1)./(epsic+2));

    kappaic=   0.434 * 6*pi/lmda_0 * mv * Nic2;
end


%Code 4.3: Relative Dielectric Constant of Pure Ice
%Description: Code computes the real and imaginary parts of the relative
%dielectric constant of pure ice at any temperature -40C<T<0C and 
%frequency 1<f<1000GHz
%Input Variables:
    %T: Temperature in degree 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-3
%MATLAB Code: RelDielConst_PureIce.m

%Example call: [A B] = RelDielConst_PureIce(T,f)
%Computes the real and imaginary components of the permitivity of pure ice
    %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_PureIce(T,f)

T = T + 273; % represent temperature in Kelvin

theta = (300/T) - 1;

B1 = 0.0207;
B2 = 1.16e-11;
b = 335;

alpha = (0.00504 + 0.0062*theta)*exp(-22.1*theta);

betaM = (B1/T) * exp(b/T) / (exp(b/T)-1)^2  + B2.*f.^2;
delBeta = exp(-9.963 + 0.0372*(T-273.16));

beta = betaM + delBeta;

% epsr = 3.1884 + 9.1e-4*(T-273) .*f./f;

epsr = 3.1884 + 9.1e-4 *(T-273);

epsi = alpha./f + beta.*f;
end
