function path_loss_db = two_ray_loss(d, fc, ht, hr, epsilon, polarization)
% two_ray_loss - path loss for a two ray channel
%
% Usage:
%   path_loss = two_ray_loss(d, fc, hb, hm, epsilon, polarization)
%
% Inputs:
%   d - distance (meters) between (base of) transmiter and receiver
%   fc - carrier frequency
%   ht - transmitter antenna height (meters)
%   hr - receiver antenna height (meters)
%   epsilon - relative dielectric constant of grond (default: 15)
%   polarization - either 'vertical' or 'horizontal' (default: 'horizontal')
%
% Output:
%   path_loss in dB

%% set default parameters
if nargin < 5 || isempty(epsilon)
    epsilon = 15;
end

if nargin < 6 || isempty(polarization)
    polarization = 'horizontal';
end

if ~strcmp(polarization, 'horizontal') && ~strcmp(polarization, 'vertical')
    error('Polarization must be either ''horizontal'' or ''vertical'', not %s', polarization);
end

%% constants and such
c = 3e8;       % speed of light
lambda = c/fc; % wavelength

%% set up geometry
tx_loc = [0 ht];        % tx antenna
rx_loc = [d hr];        % rx antenna
ref_rx_loc = [d -hr];   % reflected rx antenna

%% compute path lengths
l_los = norm(tx_loc - rx_loc);
l_ref = norm(tx_loc - ref_rx_loc);

%% angle of reflection
% the slope of the line from tx to reflected rx equal the tangent of the
% angle between that ray and teh ground plane
slope = (ht + hr)/d;
ref_angle = atan(slope);

%% reflection cefficient
if strcmp(polarization, 'horizontal')
    Z = sqrt(epsilon - cos(ref_angle)^2);
else
    Z = 1/epsilon*sqrt(epsilon - cos(ref_angle)^2);
end

R = (sin(ref_angle) - Z)/(sin(ref_angle) + Z);

%% path loss
% phase difference due to delay difference
d_phi = 2*pi/lambda*(l_ref - l_los);

% inverse of path loss on linear scale 
path_gain_lin = (lambda/(4*pi))^2 * abs(1/l_los + R*exp(-1j*d_phi)/l_ref)^2;

% turn into path loss in db
path_loss_db = -10*log10(path_gain_lin);

