function [total, surf_specular, sub_specular, volume] = simple_two_layer(x, eps1, eps2, s1, s2, A)
% Models the backscatter of surface, volume, and a subsurface with these properties:
% x: inc angle (rad)
% eps1: effective dialectric constant of the surface material
% eps2: effective dialectric constant of the subsurface material
% s1: rms roughness of surface (rad)
% s2: rms roughness of subsurface (rad)
% A: albedo

%Set surface and subsurface correlation functions.
%1 = exponential; 2 = gaussian; 3 = hagfors
surf_corr = 1;
sub_corr = 1;

surf_specular = specular(x,eps1,s1, surf_corr); %Surface backscatter -- no different from a 1 layer model.

n1 = sqrt(eps1); %refractive index of surf
n2 = sqrt(eps2); %refractive index of sub

x2 = asin(sin(x)/n1); %inc angle of the beam in the medium. adjusted by Snell's law
 
[~, ~, ~, T2] = transmission_function(eps1, eps2, 1, 1, x); %Using T2 because it depends only on surface dialectric


sub_specular = T2 .* specular(x2, eps2, s2, sub_corr);

%Next, calculate the volume scattering
%T = choose VV polarization

volume = volume_scatter(A, x);

total = surf_specular + sub_specular + volume;

end


function [specular] = specular(x,diel,s, corr)
%Calculates specular reflection off of a surface
%Input parameters:
%   diel = dielectric constant
%   rms = rms height (meters)
%   L = correlation length (meters)
%   theta = incidence angle (radians)
%   corr: 1 = exponential; 2 = gaussian; 3 = hagfors

if x >= pi/2
    specular = NaN;
    return
end

ref = ((1 - sqrt(diel))/(1 + sqrt(diel)))^2; %Fresnel reflection coefficient

switch corr
    case 1 %Exponential
        specular = ( (3.*(ref)) ./ ((s.^2).* (cos(x).^4) )) .* exp( (-(sqrt(6) .* tan(x)) ) ./ ((s) ));
    case 2 %Gaussian
        specular = ((ref) ./ ((s.^2).* (cos(x).^4) )) .* exp( (-(tan(x).^2)) ./ ((s.^2) ));
    case 3 %Hagfors
        C = 1./(s.^2); %width parameter
        specular = ((ref.*C) ./ 2) .* ((((cos(x).^4) + C.*(sin(x).^2))).^(-3/2)); 
end
end


