%% 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 12.7: Coherent Emissivity of a two-layer structure

%Description: Code computes emissivity of a perfectly homogeneous layer
%(such as ice or oil) above another homogeneous medium, bounded by
%perfectly smooth planar surfaces.

%Input Variables: 
    %eps2: dielectric constant of middle layer
    %eps3: dielectric constant of bottom layer
    %theta_i: incidence angle in air (deg)
    %d: layer thickness (m)
    %f: frequency (GHz)
    
%Output Products:
    %e_v_coh(theta_i): v-polarized emissivity
    %e_h_coh(theta_i): h-polarized emissivity
    
%Book Reference: Section 12-12.1

%Matlab code

function [e_v_coh e_h_coh] = CohEmissivity_TwoLayer(eps2, eps3, theta_i, d, f) 

eps1 = 1;
theta1 = theta_i*pi/180; %transform to radians
n2 = sqrt(eps2); %calc index of refraction
n1 = 1.0;

theta2 = asin(abs(n1/n2) * sin(theta1)); % incidence angle in medium 2 (refracted)

k = 2*pi *f/ 0.3; %wavenumber 
%-- calculate reflectivies at the two interfaces

%- 12 interface
[rhoh12 rhov12 gammah gammav tauh tauv Th Tv] = ReflTransm_PlanarBoundary(eps1, eps2, theta_i);

%- 23 interface
[rhoh23 rhov23 gammah gammav tauh tauv Th Tv] = ReflTransm_PlanarBoundary(eps2, eps3, theta2*180/pi);


alpha = k* imag(sqrt(eps2));
beta = k * real(sqrt(eps2));

gamma = 1i * alpha +  beta;

rho_v = (rhov12 + rhov23 * exp(-2*1i*gamma*d*cos(theta2))) ./ ...
    (1 + rhov12 * rhov23 * exp(-2*1i*gamma*d*cos(theta2))); 

rho_h = (rhoh12 + rhoh23 * exp(-2*1i*gamma*d*cos(theta2))) ./ ...
    (1 + rhoh12 * rhoh23 * exp(-2*1i*gamma*d*cos(theta2))); 

e_v_coh = 1- (abs(rho_v)^2);
e_h_coh = 1- (abs(rho_h)^2);

end


%Code 2.3: Oblique Reflection and Transmission @ Planar Boundry
%Description: Code computes the reflection coefficients, transmission
%coefficients, reflectivities and transmissivities for incidence in
%medium (medium 1) upon the planar boundary of a lossless or lossy
%medium (medium 2) at any incidence angle, for both h and v polarizations

%Input Variables:
    %eps1: eps1r -j*eps1i: relative dielectric constant of medium 1
    %eps2 = eps2r-j*eps2i: relative dielectric constant of medium 2
    %theta1d: incidence angle in medium 1 in degrees

%Output Products:
    %rhoh: reflection coefficient for h pol
    %rhov: reflection coefficient for v pol
    %gammah:reflectivity for h pol
    %gammav: reflectivity for v pol
    %tauh: transmission coefficient for h pol
    %tauv: transmission coefficient for v pol
    %Th: transmissivity for h pol
    %Tv:transmissivity for v pol
%Book Reference: Sections 2-7 & 2-8

%Example call: [rhoh rhov gammah gammav tauh tauv Th Tv] = ReflTransm_PlanarBoundary(eps1, eps2, theta1d)
%Computes 

%MATLAB Code

function [rhoh rhov gammah gammav tauh tauv Th Tv] = ReflTransm_PlanarBoundary(eps1, eps2, theta1d)
    
    theta1 = degtorad(theta1d);
    
    sin_theta2 = sqrt(eps1)/sqrt(eps2).*sin(theta1);
    cos_theta2 = sqrt(1 - sin_theta2.^2);
    
    rhoh = (sqrt(eps1).*cos(theta1)-sqrt(eps2).*cos_theta2) ./ (sqrt(eps1).*cos(theta1) + sqrt(eps2).*cos_theta2);
    rhov = (sqrt(eps1).*cos_theta2-sqrt(eps2).*cos(theta1)) ./ (sqrt(eps1).*cos_theta2 + sqrt(eps2).*cos(theta1));

    tauh = 1 + rhoh;
    tauv = (1 + rhov).*(cos(theta1)./cos_theta2);
       

    gammah = abs(rhoh).^2;
    gammav = abs(rhov).^2;
        
    Th = 1-gammah;
    Tv = 1-gammav;
    
  end
