%% 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 10.4: I2EM Backscattering from Periodic Sinusoidal Surface

%Description: Code computes sigma_0_vv, hh, and hv, for a sinusoidal
%periodic surface with a random smaller-scale component.

%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:
    %sigma_0_vv(theta_0, phi_0), sigma_0_hh, sigma_0_hv
    
    %To maintain reasonable execution time for this code, calculation of
    %sigma_0_hv is switched off here. However, it can be activated with
    %ease (check the function called Calc_scatt_coeff listed below).
    
%Book Reference: Section 10-4.2

%Matlab code:

function [sigma_0_vv sigma_0_hh sigma_0_vh] = 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

sig_ss_vv = zeros(ny,1); 
sig_ss_hh = sig_ss_vv; sig_ss_vh = sig_ss_vv;

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

sigma_0_vv = trapz(y,sig_ss_vv) / Gmm; % integrate
sigma_0_hh = trapz(y,sig_ss_hh) / Gmm;
sigma_0_vh = trapz(y,sig_ss_vh) / Gmm;

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

function [sig_ss_vv, sig_ss_hh, sig_ss_vh] = Calc_scatt_coeff(Zx,Zy,th,ph,eps, f, s, l, sp)

%-It calculates the scattering coefficients 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 backscattering coefficients due to small scale roughness

  % --- this part is not used anymore. It still does work ----%
% % % [sig_s_vv, sig_s_hh, sig_s_vvhh, sig_s_hv] = Calc_localScat_PO(thpr, conj(eps), f, s, l);
% % % 
% % % sig_ss_vv = dot(v, v_pr).^4 .* sig_s_vv + dot(v, h_pr).^4 .*sig_s_hh + ...
% % %     2 * dot(v, h_pr).^2 .* dot(v, v_pr).^2 .* sig_s_vvhh;
% % % 
% % % sig_ss_hh = dot(h, v_pr).^4 .* sig_s_vv + dot(h, h_pr).^4 .*sig_s_hh + ...
% % %     2 * dot(h, h_pr).^2 .* dot(h, v_pr).^2 .* sig_s_vvhh;
% % % 
% % % sig_ss_vh = dot(v, v_pr).^2 .* dot(v_pr, h).^2 .* sig_s_vv + dot(v, h_pr).^2 ...
% % %     .* dot(h_pr, h).^2 .*sig_s_hh + 2 * dot(v, v_pr) .* dot(h, h_pr) ... 
% % %     .* dot(h, v_pr) .* dot(h_pr, v) .* sig_s_vvhh + ( dot(v, v_pr).* ...
% % %     dot(h_pr, h) + dot(v, h_pr) .* dot(v_pr, h)).^2 .* sig_s_hv;
  %-----------------------------------------------------------%

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

%--calculate co-pol using I2EM_Bistat_model. Model is run in backscater
[sig_s_vv sig_s_hh] = I2EM_Bistat_model(f, s, l, thpr_deg, thpr_deg, 180, eps, sp, 1);



  %-- Calculation of x-pol can be done as shown below. However, it will
  %slow down the entire code considerably. It have been commented out here.
  %However, the user may choose to use it by removing the comments.
%--calculate x-pol using IEMX_model. 
%  [sig_s_hv] = IEMX_model(f, s, l, thpr_deg, eps, sp,1, 0);
sig_s_hv = 0;

 
 
sig_s_vv = 10.^(sig_s_vv/10); % transform to linear scale
sig_s_hh = 10.^(sig_s_hh/10);
sig_s_hv = 10.^(sig_s_hv/10);


%-- The following expressions keep only the strong terms.
sig_ss_vv = dot(v, v_pr).^4 .* sig_s_vv ;

sig_ss_hh = dot(h, h_pr).^4 .*sig_s_hh ;

sig_ss_vh = (dot(v, v_pr).* dot(h_pr, h) + ...
            dot(v, h_pr) .* dot(v_pr, h)).^2 .* sig_s_hv;

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

% This function is no longer in use now. 
function [sig_s_vv, sig_s_hh, sig_s_vvhh, sig_s_hv] = Calc_localScat_PO(th, eps, f, s, l)

%-calculates the noncoherent backscattering coefficients of the small scale roughness
%surface based on Physical Optics model

[rvv, rhh]=refl_coef(th, 1, eps); %calculates reflection coefficients at selected angle

wn = 2*pi * f / 0.3; % wavenumber

%-calculate coherent component at th = 0 deg

d0 = 0.01; % width of the illuminated area

th_err = 0.005 * pi/180;

if th <=th_err
    sig_s_vv_c = (4*wn^2*d0^2 * abs(rvv)^2 /pi) .*exp(-(2*wn*s*cos(th))^2);
    sig_s_hh_c = (4*wn^2*d0^2 * abs(rhh)^2 /pi) .*exp(-(2*wn*s*cos(th))^2);
else
    sig_s_vv_c = 0;
    sig_s_hh_c = 0;
end

%--calculate incoherent component
xx = (2 * wn * s* cos(th) )^2;
sum = 0;

for ii = 1:31
    sum = sum + (xx)^ii ./ factorial(ii);
end

cons = 2*wn * sin(th);
intg = 0;

% delx = 0.0001;
% for ii = 1: 100000
%     x = (ii-1) * delx;
%     rho = exp(-x^2/(l1^4 + x^2*l2^2)^0.5);
%     ss = rho*besselj(0,x*cons)*x*delx;
%     intg = intg + ss;
%  %   x, rho, ss
% end
% intg = intg * sum;



% load smallscale_corrfunc.mat;
% b = max(xmod); 
% a = 0;
% xi = linspace(a,b,10000);
% rho = interp1(xmod, tt_r, xi);

b = 1; 
a = 0;
xi = linspace(a,b,10000);
% rho = exp(-xi.^2./(l1.^4 + xi.^2 .*l2.^2).^0.5); %gaussian/exponential
% rho = exp(-xi.^2 / l1.^2); %gaussian
% rho = exp(-xi / l1); %Exponential

 rho = (1 + xi.^2 ./l.^2) .^(-1.5); % ulaby's


[JJ, ierr] = besselj(0, xi .*cons); 
yi = rho.* JJ .*xi;

if ierr ~= 0
    ierr
end

%I1 = quad(@(x) interp1(xi, yi, x, 'cubic',0), a, b, 1.0e-12);

I1 = trapz(xi, yi);

intg = I1 .*sum;


%--calculate the backscattering coefficients

sig_s_vv_n  = 2 *( wn*abs(rvv)*cos(th))^2 * exp(-(2*wn*s*cos(th))^2) * intg;
sig_s_hh_n  = 2 *( wn*abs(rhh)*cos(th))^2 * exp(-(2*wn*s*cos(th))^2) * intg;


sig_s_vv = sig_s_vv_c + sig_s_vv_n;
sig_s_hh = sig_s_hh_c + sig_s_hh_n;


sig_s_hv = 0;

sig_s_vvhh = sig_s_vv * abs(real(rvv * conj(rhh)))/abs(rvv)^2;

end
function [rho_v, rho_h]=refl_coef(the1, eps1, eps2)

% calculates the v and h-polarized reflection coefficients of a plane
% dielectric surface

n1 = sqrt(eps1);
n2 = sqrt(eps2);
costh2 = sqrt(1 - (n1.*sin(the1)./n2).^2);


rho_v = -(n2.*cos(the1) - n1.*costh2)./(n2.*cos(the1)+n1.*costh2);
rho_h = (n1.*cos(the1) - n2.*costh2)./(n1.*cos(the1)+n2.*costh2);
end


%Code 10.3: I2EM Bistatic Scattering from Single-Scale Random Surface

%Description: Code computes sigma_0_vv(thi, ths,phs) and
%sigma_0_hh(thi, ths,phs) for single-scale random surface

%Input Variables: 
    %er: complex dielectric constant of the scattering medium
    %thi: Incidence angle (deg)
    %ths: Scattering angle (deg)
    %phs: Relative azimuth angle (deg)
    
    %sp: type of correlation function: 1- exponential, 2- Gaussian
    %       3- x-power,  4- x-exponential
    %xx: coefficient (>1) needed for x-power and x-exponential correl. fnc.
    %sig: rms height (m)
    %L: correlation length (m)
    %fr: frequency (GHz)
    
%Output Products:
    % sigma_0_vv, sigma_0_hh in dB

%Book Reference: Section 10-3.9

%Matlab Code: 

function [sigma_0_vv sigma_0_hh] = I2EM_Bistat_model(fr, sig, L, thi, ths, phs, er, sp, xx)


error = 1.0e8;
sig = sig * 100; % change from m to cm scale
L = L * 100;
mu_r = 1; % relative permeability

k = 2*pi *fr/30; % wavenumber in free space. Speed of light is in cm/sec
theta = thi .*pi/180; % transform to radian
phi = 0;
thetas = ths * pi/180;
phis = phs * pi/180;

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

ks2 = ks .* ks; 

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

sf = sin(phi);
cf = cos(phi);

ss = sin(thetas);
css = cos(thetas);

cfs = cos(phis);
sfs = sin(phis);

s2 = s * s;
% sq = sqrt(er - s2);

kx = k .* s .*cf;
ky = k .* s .*sf;
kz = k .* cs;

ksx = k .* ss .*cfs;
ksy = k .* ss .*sfs;
ksz = k .* css;

%-- reflection coefficients
rt = sqrt(er - s2);
Rvi = (er .*cs - rt) ./(er.*cs +rt);
Rhi = (cs - rt)./(cs + rt);


wvnb = k .* sqrt( (ss .*cfs - s .*cf).^2 + (ss .* sfs - s .* sf).^2 );

Ts = 1;

while error > 1.0e-8,
    Ts = Ts + 1;
    error = (ks2 .*(cs + css).^2 ).^Ts ./ factorial(Ts); 
end

%---------------- calculating roughness spectrum -----------

[wn rss] = roughness_spectrum(sp,xx, wvnb, sig, L, Ts);


%----------- compute R- transition ------------

Rv0 = (sqrt(er)-1) ./(sqrt(er)+1);
Rh0 = -Rv0;

Ft = 8 * Rv0.^2 * ss *(cs + sqrt(er - s2))./(cs .* sqrt(er - s2));
a1 = 0;
b1 = 0;
for n = 1:Ts
    a0 = (ks .*cs).^(2*n) ./factorial(n);
    a1 = a1 + a0 *wn(n);
    b1 = b1 + a0 * (abs(Ft./2 + 2.^(n+1) .*Rv0./cs .*exp(-(ks.*cs).^2))).^2 ...
        * wn(n);
end
St = 0.25 * (abs(Ft)).^2 * a1 ./ b1;

St0 = 1 ./ (abs(1 + 8 *Rv0./(cs .* Ft))).^2;

Tf = 1 - St ./St0;


%----------- compute average reflection coefficients ------------
%-- these coefficients account for slope effects, especially near the
%brewster angle. They are not important if the slope is small.

sigx = 1.1 .*sig/L;
sigy = sigx;
xxx = 3*sigx;

Rav = dblquad(@(Zx, Zy)Rav_integration(Zx, Zy, cs,s,er,s2,sigx, sigy),-xxx,xxx, -xxx, xxx );

Rah = dblquad(@(Zx, Zy)Rah_integration(Zx, Zy, cs,s,er,s2,sigx, sigy),-xxx,xxx, -xxx, xxx );

Rav = Rav ./(2*pi * sigx * sigy);
Rah = Rah ./(2*pi * sigx * sigy);



%-- select proper reflection coefficients

if thi == ths && phs==180, %i.e. operating in backscatter mode
    Rvt = Rvi + (Rv0 - Rvi) .*Tf;
    Rht = Rhi + (Rh0 - Rhi) .*Tf;

else        % in this case, it is the bistatic configuration and average R is used
%     Rvt = Rav + (Rv0 - Rav) .* Tf;
%     Rht = Rah + (Rh0 - Rah) .* Tf;
    Rvt = Rav;
    Rht = Rah;
end

fvv = 2 .* Rvt .*(s .* ss - (1 + cs .* css) .* cfs)./(cs + css);
fhh = -2 .* Rht .*(s .* ss - (1 + cs .* css) .* cfs)./(cs + css);

%------- Calculate the Fppup(dn) i(s) coefficients ----
[Fvvupi, Fhhupi] = Fppupdn_is_calculations(+1, 1, Rvi,Rhi,er,k,kz,ksz,s,cs,ss,css,cf,cfs,sfs);
[Fvvups, Fhhups] = Fppupdn_is_calculations(+1, 2, Rvi,Rhi,er,k,kz,ksz,s,cs,ss,css,cf,cfs,sfs);
[Fvvdni, Fhhdni] = Fppupdn_is_calculations(-1, 1, Rvi,Rhi,er,k,kz,ksz,s,cs,ss,css,cf,cfs,sfs);
[Fvvdns, Fhhdns] = Fppupdn_is_calculations(-1, 2, Rvi,Rhi,er,k,kz,ksz,s,cs,ss,css,cf,cfs,sfs);


qi = k .* cs;
qs = k .* css;

%----- calculating  Ivv and Ihh ----

Ivv = zeros(Ts, 1); Ihh = Ivv;
    
for n = 1:Ts
    Ivv(n) = (kz + ksz).^n .* fvv .* exp(-sig^2 .* kz .* ksz) + ...
        0.25*(Fvvupi .*(ksz-qi).^(n-1) .*exp(-sig^2 .*(qi.^2 - qi.*(ksz-kz)))+ ...
        Fvvdni .*(ksz+qi).^(n-1) .*exp(-sig^2 .*(qi.^2 + qi.*(ksz-kz)))+ ...
        Fvvups .*(kz+qs).^(n-1) .*exp(-sig^2 .*(qs.^2 - qs.*(ksz-kz)))+ ...
        Fvvdns .*(kz-qs).^(n-1) .*exp(-sig^2 .*(qs.^2 + qs.*(ksz-kz))));

    
    Ihh(n) = (kz + ksz).^n .* fhh .* exp(-sig^2 .* kz .* ksz) + ...
        0.25*(Fhhupi .*(ksz-qi).^(n-1) .*exp(-sig^2 .*(qi.^2 - qi.*(ksz-kz)))+ ...
        Fhhdni .*(ksz+qi).^(n-1) .*exp(-sig^2 .*(qi.^2 + qi.*(ksz-kz)))+ ...
        Fhhups .*(kz+qs).^(n-1) .*exp(-sig^2 .*(qs.^2 - qs.*(ksz-kz)))+ ...
        Fhhdns .*(kz-qs).^(n-1) .*exp(-sig^2 .*(qs.^2 + qs.*(ksz-kz))));
end


%-- Shadowing function calculations

if thi==ths && phs==180 %i.e. working in backscatter mode
    ct = cot(theta);
    cts = cot(thetas);
    rslp = rss;
    ctorslp = ct / sqrt(2) ./rslp;
    ctsorslp = cts / sqrt(2) ./rslp;
    shadf = 0.5 *(exp(-ctorslp.^2) ./ sqrt(pi)./ctorslp - erfc(ctorslp));
    shadfs = 0.5 *(exp(-ctsorslp.^2) ./ sqrt(pi)./ctsorslp - erfc(ctsorslp));
    ShdwS = 1./(1 + shadf + shadfs); 
else
 ShdwS = 1;
end

%------- calculate the values of sigma_note --------------

sigmavv = 0;sigmahh =0; 

for n = 1:Ts
    a0 = wn(n) ./factorial(n) .*sig.^(2*n);
    
    sigmavv = sigmavv+ abs(Ivv(n)).^2 .*a0;
    sigmahh = sigmahh+ abs(Ihh(n)).^2 .*a0;

end

sigmavv = sigmavv * ShdwS * k^2 ./2 * exp(-sig.^2 .*(kz.^2 +ksz.^2));  
sigmahh = sigmahh * ShdwS * k^2 ./2 * exp(-sig.^2 .*(kz.^2 +ksz.^2));  


ssv = 10 * log10(sigmavv); 
ssh = 10 * log10(sigmahh); 

sigma_0_vv = ssv;
sigma_0_hh = ssh;

end

function [vv, hh] = Fppupdn_is_calculations(ud, is, Rvi,Rhi,er,k,kz,ksz,s,cs,ss,css,cf,cfs,sfs)

if is==1
    Gqi = ud .* kz;
    Gqti = ud .*k .*sqrt(er-s.^2);
    qi = ud .* kz;
    
    c11 = k .* cfs .*(ksz - qi);
    c21 = cs .*(cfs .*(k^2 .*s.*cf.*(ss .*cfs - s .* cf) + Gqi.*(k .* css - qi)) ...
        + k^2 .*cf .* s .*ss .*sfs^2);
    c31 = k.*s.*(s.*cf.*cfs.*(k.*css-qi) - Gqi.*(cfs.*(ss.*cfs -s.*cf)+ ss .*sfs^2));
    c41 = k .*cs.*(cfs.*css.*(k.*css - qi) + k .*ss.*(ss.*cfs-s.*cf));
    c51 = Gqi.*(cfs .*css.*(qi-k.*css) - k .*ss.*(ss.*cfs-s.*cf));
    
    c12 = k .* cfs .*(ksz - qi);
    c22 = cs .*(cfs .*(k^2 .*s.*cf.*(ss .*cfs - s .* cf) + Gqti.*(k .* css - qi)) ...
        + k^2 .*cf .* s .*ss .*sfs^2);
    c32 = k.*s.*(s.*cf.*cfs.*(k.*css-qi) - Gqti.*(cfs.*(ss.*cfs -s.*cf)- ss .*sfs^2));
    c42 = k .*cs.*(cfs.*css.*(k.*css - qi) + k .*ss.*(ss.*cfs-s.*cf));
    c52 = Gqti.*(cfs .*css.*(qi-k.*css) - k .*ss.*(ss.*cfs-s.*cf));    
end

if is==2
    Gqs = ud .* ksz;
    Gqts = ud .*k .*sqrt(er-ss.^2);
    qs = ud .* ksz;
    
    c11 = k .* cfs .*(kz + qs);
    c21 = Gqs .*(cfs.*(cs.*(k.*cs+qs)-k.*s.*(ss .*cfs-s.*cf))-k.*s.*ss.*sfs.^2);
    c31 = k .*ss.*(k.*cs.*(ss.*cfs - s.*cf)+ s.*(kz+qs));
    c41 = k.*css.*(cfs.*(cs.*(kz+qs)-k.*s.*(ss.*cfs-s.*cf))-k.*s.*ss.*sfs.^2);
    c51 = -css .*(k.^2 .*ss .*(ss.*cfs -s.*cf)+ Gqs.*cfs.*(kz+qs));
         
    c12 = k .* cfs .*(kz + qs);
    c22 = Gqts .*(cfs.*(cs.*(kz+qs)-k.*s.*(ss .*cfs-s.*cf))-k.*s.*ss.*sfs.^2);
    c32 = k .*ss.*(k.*cs.*(ss.*cfs - s.*cf)+ s.*(kz+qs));
    c42 = k.*css.*(cfs.*(cs.*(kz+qs)-k.*s.*(ss.*cfs-s.*cf))-k.*s.*ss.*sfs.^2);
    c52 = -css .*(k.^2 .*ss .*(ss.*cfs -s.*cf)+ Gqts.*cfs.*(kz+qs));
end

q=kz;
qt = k .*sqrt(er - s.^2);

vv = (1+Rvi) .*( -(1-Rvi) .*c11 ./q + (1+Rvi) .*c12 ./ qt) + ...
    (1 - Rvi) .*( (1-Rvi) .*c21 ./q - (1+Rvi) .*c22 ./ qt) + ...
    (1+Rvi) .*( (1-Rvi) .*c31 ./q - (1+Rvi) .*c32 ./er ./qt) + ...
    (1 - Rvi) .*( (1+Rvi) .*c41 ./q - er.*(1 - Rvi) .*c42 ./ qt) + ...
    (1+Rvi) .*( (1+Rvi) .*c51 ./q - (1-Rvi) .*c52 ./ qt);

hh = (1 + Rhi) .*( (1-Rhi) .*c11 ./q - er.*(1+Rhi) .*c12 ./ qt) - ...
    (1 - Rhi) .*( (1-Rhi) .*c21 ./q - (1+Rhi) .*c22 ./ qt) - ...
    (1 + Rhi) .*( (1-Rhi) .*c31 ./q - (1+Rhi) .*c32 ./qt) - ...
    (1 - Rhi) .*( (1+Rhi) .*c41 ./q - (1 - Rhi) .*c42 ./ qt) - ...
    (1 + Rhi) .*( (1+Rhi) .*c51 ./q - (1-Rhi) .*c52 ./ qt);

end


function Rav = Rav_integration(Zx, Zy, cs,s,er,s2,sigx, sigy)

A = cs + Zx .* s;
B = er .* (1 + Zx.^2 + Zy.^2);
CC = s2 - 2.*Zx .*s .*cs + Zx.^2 .* cs^2 + Zy.^2;

Rv = (er.*A - sqrt(B-CC))./(er.*A + sqrt(B-CC)); 

pd = exp(-Zx.^2 ./(2*sigx.^2) -Zy.^2 ./(2*sigy.^2));
Rav = Rv .* pd;

end

function Rah = Rah_integration(Zx, Zy, cs,s,er,s2,sigx, sigy)

A = cs + Zx * s;
B = er .* (1 + Zx.^2 + Zy.^2);
CC = s2 - 2.*Zx .*s .*cs + Zx.^2 .* cs^2 + Zy.^2;

Rh = (A - sqrt(B-CC))./(A + sqrt(B-CC)); 

pd = exp(-Zx.^2./(2*sigx.^2) -Zy.^2./(2*sigy.^2));
Rah = Rh .* pd;

end

function [wn rss] = roughness_spectrum(sp,xx, wvnb, sig, L, Ts)

wn = zeros(Ts,1);
%-- exponential correl func
if sp ==1,
    for n = 1: Ts
        wn(n) = L.^2 / n.^2 .* (1+(wvnb.*L/n).^2).^(-1.5);
    end
    rss= sig ./L;
end

%-- gaussian correl func
if sp ==2,
    for n = 1: Ts
        wn(n) =  L.^2/(2*n) .* exp(-(wvnb.*L).^2 ./(4*n));
    end
    rss = sqrt(2) * sig./L;
end

%-- x-power correl func
if sp == 3,
   
    for n =1: Ts
      if wvnb==0, 
          wn(n) = L.^2 ./(3 * n - 2);
      else
          wn(n) = L.^2 * (wvnb.*L).^(-1 + xx*n) .* besselk(1 -xx*n, wvnb*L) ...
              ./ (2.^(xx*n -1) .* gamma(xx .*n)); 
      end
    end
%    if xx == 1.5
        rss  = sqrt(xx *2) .*sig ./L;
%    else
%        rss = 0;
%    end
end

%-- x- exponential correl func

if sp == 4,
    for n =1: Ts
        tmp = quad(@(z) x_exponential_spectrum(z,wvnb,L,n,xx), 0,9);
        wn(n) = L.^2 ./ n.^(2/xx) .* tmp;
    end
    rss  =  sig ./ L;
end

end

function [tmp] = x_exponential_spectrum(z,wvnb,L,n,xx)

tmp = exp(-abs(z).^xx) .* besselj(0, z*wvnb*L./(n.^(1/xx))).*z;
end


function sigvh = IEMX_model(fr, sig, L, theta_d, er, sp,xx, auto)

sig = sig * 100 ; % change to cm scale
L = L * 100; % change to cm scale;


%- fr: frequency in GHz
%- sig: rms height of surface in cm
%- L: correlation length of surface in cm
%- theta_d: incidence angle in degrees
%- er: relative permittivity
%- sp: type of surface correlation function

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+ 0.001);

s2 = s.^2;

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

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

rvh = (rv - rh) ./2;

%-- rms slope values
sig_l = sig/L;
if sp==1                %-- exponential correl func
    rss = sig_l;
end
if sp==2                %-- Gaussian correl func
    rss = sig_l * sqrt(2);
end

if sp==3                 %-- 1.5-power spectra correl func
    rss = sig_l * sqrt(2*xx);
end

%--- Selecting number of spectral components of the surface roughness
if auto == 0 
n_spec = 15; % number of terms to include in the surface roughness spectra
end

if auto == 1
   n_spec = 1; 
   while error > 1.0e-8,
    n_spec = n_spec + 1;
    error = (ks2 .*( 2*cs).^2 ).^n_spec ./ factorial(n_spec); 
   end

end

%-- calculating shadow consideration in single scat (Smith, 1967)

ct = cot(theta+ 0.001);
farg = ct /sqrt(2) ./rss;
gamma = 0.5 *(exp(-farg.^2) / 1.772 / farg - erfc(farg));
Shdw = 1 ./ (1 + gamma);

%-- calculating multiple scattering contribution
%------ a double integration function

svh = dblquad(@(r,phi)xpol_integralfunc(r, phi, sp,xx, ks2, cs,s, kl2, L, er, rss, rvh, n_spec), 0.1, 1, 0, pi);

sigvh = 10*log10(svh .* Shdw);

end

function y = xpol_integralfunc(r, phi, sp, xx, ks2, cs,s, kl2, L, er, rss, rvh, n_spec)

cs2 = cs .^2;

r2 = r.^2;
nr = length(r);

sf = sin(phi);
csf = cos(phi);
rx = r .* csf;
ry = r .* sf;

%-- calculation of the field coefficients
rp = 1 + rvh;
rm = 1 - rvh; 

q = sqrt(1.0001 - r2);
qt = sqrt(er - r2);

a = rp ./q;
b = rm ./q;
c = rp ./qt;
d = rm ./qt;

%--calculate cross-pol coefficient
B3 = rx .* ry ./cs;
fvh1 = (b-c).*(1- 3*rvh) - (b - c./er) .* rp; 
fvh2 = (a-d).*(1+ 3*rvh) - (a - d.*er) .* rm;
Fvh = ( abs( (fvh1 + fvh2) .*B3)).^2;


%-- calculate shadowing func for multiple scattering 
au = q ./r ./1.414 ./rss;
fsh = (0.2821./au) .*exp(-au.^2) -0.5 .*(1- erf(au));
sha = 1./(1 + fsh); 

%-- calculate expressions for the surface spectra
wn = spectrm1(sp, xx, kl2, L, rx, ry, s, n_spec, nr);
wm = spectrm2(sp, xx, kl2, L, rx, ry, s, n_spec, nr);


%--compute VH scattering coefficient
acc = exp(-2* ks2 .*cs2) ./(16 .* pi);
vhmnsum = zeros(1,nr);
for n = 1: n_spec
    for m = 1: n_spec
        vhmnsum = vhmnsum + wn(n,:).*wm(m,:) .*(ks2*cs2).^(n+m) ...
            ./factorial(n)./factorial(m); 
    end
end

VH = 4 * acc .* Fvh .* vhmnsum .*r;
y = VH .* sha;

end

function wn = spectrm1(sp, xx, kl2, L, rx, ry, s, np, nr)

wn = zeros(np,nr);

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

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

if sp== 3 % x-power
    for n = 1: np
        wn(n,:) = kl2 ./(2.^(xx*n-1) .*gamma(xx*n)).* ( ( (rx-s).^2 ...
          + ry.^2).*L).^(xx*n-1) .* besselk(-xx*n+1, L*((rx-s).^2 + ry.^2));
    end
end

end

function wm = spectrm2(sp, xx, kl2, L, rx, ry, s, np, nr)

wm = zeros(np,nr);

if sp == 1  % exponential
    for n = 1: np
        wm(n,:) = n* kl2 ./(n.^2 + kl2 *((rx+s).^2+ry.^2)).^1.5;
    end
end

if sp == 2  %  gaussian
    for n = 1: np
        wm(n,:) = 0.5 * kl2 ./n .* exp(-kl2*((rx+s).^2 + ry.^2)/(4*n)) ;
    end
end

if sp== 3 % x-power
    for n = 1: np
        wm(n,:) = kl2 ./(2.^(xx*n-1) .*gamma(xx*n)).* ( ( (rx+s).^2 ...
          + ry.^2).*L).^(xx*n-1) .* besselk(-xx*n+1, L*((rx+s).^2 + ry.^2));
    end
end
end
