%% 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.4: Emissivity of a periodic Sinusoidal Surface

%Description: Code computes emissivity for a random surface superimposed on
%a periodic  sinusoidal surface.

%Input Variables:
    %eps: dielectric constant of the scattering medium (complex)
    %sp: correlation function of the random component:
            % 1- Exponential, 2- Gaussian
    %s: rms height of the random component (m)
    %l: correlation length of the random component (m)
    %Gmm: Spatial period of the periodic surface (m)
    %A: amplitude of periodic sinusoid (m)
    %f: frequency  in GHz
    %theta_0: Incidence angle (deg)
    %phi_0: Azimuth angle(Fig. 10-24) (deg)

%Output Products:
    %e_v(theta_0, phi_0): v-polarized emissivity
    %e_h(theta_0, phi_0): h-polarized emissivity    
   
%Book Reference: Section 12-4

%Matlab code:

function [ev eh] = Emissivity_I2EM_Periodic(theta_0,phi_0 ...
    , eps, f, s, l, sp, Gmm, A)

ny = 100; % number of segments or points used to discretize the surface
y = linspace(0, Gmm, ny);

Zx = 0; %partial derivative of periodic surface along x

ev_ss = zeros(ny,1); 
eh_ss = ev_ss;

for ii = 1: ny
%     Zy = 2*pi/Gmm *A * sin(2*pi/Gmm * y(ii));
    Zy = 2*pi/Gmm *A * cos(2*pi/Gmm * y(ii));
    alpha = atan(Zy);
    
    [ev_ss(ii), eh_ss(ii)] = Calc_scatt_coeff(Zx,Zy, ...
        theta_0, phi_0, eps, f, s, l, sp);
     
     ev_ss(ii) = ev_ss(ii) * sec(alpha); 
     eh_ss(ii) = eh_ss(ii) * sec(alpha); 
end

ev = trapz(y,ev_ss) / Gmm / cos(theta_0*pi/180); % integrate
eh = trapz(y,eh_ss) / Gmm / cos(theta_0*pi/180); % integrate


end
%-------------------------------------------------------------------------

function [ev_ss eh_ss] = Calc_scatt_coeff(Zx,Zy,th,ph,eps, f, s, l, sp)

%-It calculates the emissivity of the differential area at a
%-given point on the periodic surface due to rough surface scattering from
%-the small scale roughness. 


%-calculate the polarization vectors in both local and global coordinates
th_rad = th*pi/180; %present angle in radian
ph_rad = ph*pi/180;

[v,h,v_pr,h_pr,thpr] = vectors2(Zx,Zy,th_rad, ph_rad);

%-calculate local emissivity  due to small scale roughness

thpr_deg = thpr*180/pi; % need this angle in degrees

%--calculate local emissivity due to small-scale roughness using I2EM

% thpr_deg

if thpr_deg > 89.5 
    ev_s = 0; % shaded area most likely, so local emissivity calc is ignored
    eh_s = 0;
else   
[ev_s eh_s] = Calc_emissivity_I2EMmodel(f, s, l, thpr_deg, eps, sp);
end

ev_ss = dot(v, h_pr).^2 .* eh_s + dot(v, v_pr).^2 .* ev_s ;
eh_ss = dot(h, h_pr).^2 .* eh_s + dot(h, v_pr).^2 .* ev_s ;

ev_ss = ev_ss * cos(thpr);
eh_ss = eh_ss * cos(thpr);

% ev_ss, eh_ss

end

function [v,h,v_pr,h_pr,thpr] = vectors2(Zx,Zy,th, ph)
%
% Calculates the local coordinates and polarization vectors in global
% coordinate system. 

snt = sin(th);
cst = cos(th);
snp = sin(ph);
csp = cos(ph);

D0 = (1 + Zx.^2 + Zy.^2).^(-0.50);

%-- calculate polarization vectors, unit normal vector, and incident wave
%directions in global coordinates

n(1) = -Zx; % unit normal vector at point on periodic surface
n(2) = -Zy;
n(3) = 1;
n = n * D0;

n = n / norm(n);

h(1) = -snp;
h(2) = csp;
h(3) = 0;

v(1) = -cst.*csp ;
v(2) = -cst.*snp;
v(3) = -snt;

ni(1) = snt .* csp; %incident wave direction
ni(2) = snt .* snp;
ni(3) = -cst;
ni = ni / norm(ni);


costhpr = - dot(ni, n);
thpr = acos(costhpr);
sinthpr = sin(thpr);

%--calculate local coordinate axes
z_pr = n;

tmp = cross(n, ni);
y_pr = tmp /norm(tmp);

x_pr = cross(y_pr, z_pr);

h_pr = y_pr;

tmp = cross(h_pr, ni);
v_pr = tmp / norm(tmp);


end



%Code 12.2: I2EM Rough Surface Emissivity

%Description: Code computes the emissivity of a homogeneous medium with
%rough surface characterized by either the exponential correlation function
%or the Gaussian correlation function. 

%Input variables:
    %er: Dielectric constant of the surface
    %theta_d: incidence angle (degrees)
    %sig: rms height (m)
    %L: correlation length (m)
    %fr: frequency in GHz
    %sp: 1- exponential and 2- gaussian correlation function

%Output Products:
    %e_h and e_v: h- and v-polarized emissivities
    
% Book Reference: Section 12-3.1

%Matlab code:

function [ev eh] = Calc_emissivity_I2EMmodel(fr, sig, L, theta_d, er, sp)

%-- This code calculates the emission from rough surfaces using the I2EM
%model

sig = sig * 100; % transform to cm
L = L* 100; %transform to cm.


error = 1.0e8;

k = 2*pi *fr/30; % wavenumber in free space. Speed of light is in cm/sec
theta = theta_d .*pi/180; % transform to radian

ks = k * sig; % roughness parameter
kl = k * L;  

ks2 = ks .* ks; 
kl2 = kl.^2;

cs = cos(theta);
s = sin(theta);

s2 = s.^2;

%-- calculation of reflection coefficints
sq = sqrt(er - s2);

rv0 = (sqrt(er) - 1) ./(sqrt(er) + 1);
rh0 = -rv0;

rv = (er *cs - sq) ./(er*cs +sq);
rh = (cs - sq)./(cs + sq);



pol = 1;
refv = dblquad(@(ths,phs)emsv_integralfunc(ths, phs, pol, sp, k, L, ks, cs,s, kl, er, sq,rv,rh), 0.0, pi/2, 0, pi);

pol = 2;
refh = dblquad(@(ths,phs)emsv_integralfunc(ths, phs, pol, sp, k, L, ks, cs,s, kl, er, sq,rv,rh), 0.0, pi/2, 0, pi);

ev = 1- refv - exp(-ks2 .* cs.*cs) .* (abs(rv)).^2;
eh = 1- refh - exp(-ks2 .* cs.*cs) .* (abs(rh)).^2;

end

function ref = emsv_integralfunc(ths, phs, pol, sp, k,L, ks, cs,s, kl, er, sq, rv,rh)

error = 1.0e3;

cs2 = cs .^2;
ks2 = ks.^2;
s2 = s .*s;
kl2 = kl.* kl;

css = cos(ths);
css2 = css .* css;

ss = sin(ths);
ss2 = ss .* ss;

sf = sin(phs);
sf2 = sf .*sf;
csf= cos(phs);

sqs = sqrt(er - ss2);
rc = (rv - rh) ./2;
tv = 1 + rv;
th = 1 + rh; 
tv2 = tv .* tv;
th2 = th .* th;

%-- calc coefficients for surface correlation spectra
wvnb = k .* sqrt(s2 - 2 .* s .*ss .* csf + ss2); 

nr = length(ths);
thsmin = min(ths); % selected as the smallest angle in ths vector
costhsmin = cos(thsmin); 

%-- calculate number of spectral components needed 
n_spec = 1; 
while error > 1.0e-3,
  n_spec = n_spec + 1;
%   error = (ks2 .*(cs + css).^2 ).^n_spec ./ factorial(n_spec); 
%---- in this case we will use the smallest ths to determine the number of
%spectral components to use. It might be more than needed for other angles
%but this is fine. This option is used to simplify calculations.
  error = (ks2 .*(cs + costhsmin).^2 ).^n_spec ./ factorial(n_spec); 
end


%-- calculate expressions for the surface spectra
wn = spectrm(sp, kl2, L, wvnb, n_spec, nr);

%-- calculate fpq!

ff = 2 .*(s .* ss - (1 + cs .* css) .* csf)./(cs + css);

fvv = rv .* ff;
fhh = -rh .* ff;

fvh = -2 * rc .* sf;
fhv = 2 .* rc .* sf;

%-- calculate Fpq and Fpqs -----
cm1 = s .*(ss - s .*csf) ./(cs2 .* css);
T = (sq .*(cs + sq) + cs.*(er.*cs + sq)) ./ (er .*cs .*(cs+sq) + sq .*(er.*cs +sq));
cm2 = css .* sq ./cs ./sqs - 1;
ex = exp(-ks2 .*cs .*css);
de = 0.5 .*exp(-ks2 .*(cs2 + css2)); 

if pol == 1
  Fvv = (er - 1) .*s2 .*tv2 .*cm1 ./ er.^2;
  Fhv = (T .*s .*s -1.+ cs./css + (er .*T .*cs .*css .*(er .*T - s .*s) - ...
     sq .* sq) ./ (T .*er .*sq .* css)) .*(1 - rc .*rc) .*sf;
 
  Fvvs = -cm2 .*sq .*tv2 .*(csf - s .*ss) ./cs2 ./er - cm2 .* sqs .*tv2 .*csf ./er ...
    -(css .*sq ./cs ./sqs ./er -1) .* ss .* tv2 .* (s - ss .*csf) ./cs;
  Fhvs = -(ss.*ss ./T -1 + css ./cs + (cs .*css .*(1- ss.*ss.*T) - ...
    T .*T.*sqs.*sqs) ./(T .*sqs.*cs)).*(1-rc .*rc).*sf;

   %-- calculate the bistatic field coefficients ---

  svv = zeros(n_spec, nr); 
  for n = 1:n_spec
    Ivv = fvv .*ex .*(ks .*(cs + css)).^n + (Fvv .*(ks .*css).^n + Fvvs .*(ks.*cs).^n) ./2;
    Ihv = fhv .*ex .*(ks .*(cs + css)).^n + (Fhv .*(ks .*css).^n + Fhvs .*(ks.*cs).^n) ./2;

    wnn = wn(n,:)./factorial(n);
    vv = wnn .* (abs(Ivv)).^2;
    hv = wnn .* (abs(Ihv)).^2;
    svv(n,:) = de .*(vv + hv) .* ss ./4/pi ./cs;
  end    

ref = sum(svv,1); % adding all n terms stores in different rows
end

if pol == 2
  Fhh = -(er -1) .* th2 .* cm1;
  Fvh = (s .*s ./T -1.+ cs./css + (cs .*css .*(1- s .*s .*T) - ...
    T .*T .*sq .* sq) ./ (T .*sq .* css)) .*(1 - rc .*rc) .*sf;

  Fhhs = cm2 .* sq .*th2 .*(csf - s .*ss) ./cs2 + cm2 .* sqs .* th2 .* csf + ...
     cm2 .*ss .*th2 .*(s - ss.*csf)./cs;
  Fvhs = -(T .* ss .* ss - 1 + css./cs + (er .*T .*cs .*css .*(er.*T - ss .*ss) ...
    - sqs .* sqs) ./(T .* er .*sqs .*cs)) .*(1- rc .*rc).*sf;


  for n = 1:n_spec
   Ihh = fhh .*ex .*(ks .*(cs + css)).^n + (Fhh .*(ks .*css).^n + Fhhs .*(ks.*cs).^n) ./2;
   Ivh = fvh .*ex .*(ks .*(cs + css)).^n + (Fvh .*(ks .*css).^n + Fvhs .*(ks.*cs).^n) ./2;
  
    wnn = wn(n,:)./factorial(n);
    hh = wnn .* (abs(Ihh)).^2;
    vh = wnn .* (abs(Ivh)).^2;
    shh(n,:) = de .*(hh + vh) .* ss ./4/pi ./cs;
  end    
  
ref = sum(shh, 1); % adding all n terms stores in different rows
end

end

%------------------------------------------------------------------
%------------------------------------------------------------------
function wn = spectrm(sp, kl2, L, wvnb, np, nr)

wn = zeros(np,nr);

if sp == 1  % exponential
    for n = 1: np       
        wn(n, :) = n* kl2 ./(n.^2 + (wvnb.*L).^2).^1.5;
    end
end

if sp == 2  %  gaussian
    for n = 1: np
        wn(n,:) = 0.5 * kl2 ./n .* exp(-(wvnb .*L).^2/(4*n)) ;
    end
end

end
