function [pl,output] = tirempl(ds, Z, fqs, varargin)
%TIREMPL   Path loss using Terrain Integrated Rough Earth Model (TIREM)
%   L = TIREMPL(R,Z,F) returns the path loss L in dB for a signal with
%   frequency F in Hz when it is propagated over terrain. The terrain is
%   specified by numeric vectors R and Z, which define distances and
%   corresponding elevation values (both in meters) along the great circle
%   path between the transmitter and receiver. The first and last values
%   correspond to the transmitter and receiver locations, respectively.
%   Terrain elevation is measured as height above mean sea level.
%
%   L = tirempl(R,Z,F,Name,Value) returns the path loss with additional
%   options specified by one or more Name-Value pairs.
%
%   [L,OUTPUT] = tirempl(___) returns the path loss L and an output
%   structure containing information about the TIREM analysis.
%
%   This function requires access to the external TIREM library. Use 
%   <a href = "matlab:help('tiremSetup')">tiremSetup</a> to set up access to TIREM.
% 
%   Path loss is calculated from free-space loss, terrain diffraction,
%   ground reflection, refraction through the atmosphere, tropospheric
%   scatter, and atmospheric absorption. The TIREM model combines physics
%   with empirical data to provide path loss estimates.
%
%   The TIREM model is valid from 1 MHz to 1000 GHz.
%
%   tirempl Name-Value pairs:
%
%   TransmitterAntennaHeight - Height of the transmitter antenna in meters,
%      specified as a numeric scalar in the range 0 to 30000. The height is
%      measured from ground elevation to the center of the antenna. The
%      default value is 10.
%
%   ReceiverAntennaHeight - Height of the receiver antenna in meters, 
%      specified as a numeric scalar in the range 0 to 30000. The height is
%      measured from ground elevation to the center of the antenna. The
%      default value is 1.
%
%   AntennaPolarization - Polarization of the transmitter and receiver 
%      antennas, specified as 'horizontal' or 'vertical'. Both antennas are
%      assumed to have the same polarization. This value is used to
%      calculate path loss due to ground reflection. The default value is
%      'horizontal'.
%
%   GroundConductivity - Conductivity of the ground in Siemens per meter 
%      (S/m), specified as a numeric scalar in the range 0.00005 to 100.
%      This value is used to calculate path loss due to ground reflection.
%      The default value is 0.005, which corresponds to average ground.
%
%   GroundPermittivity - Relative permittivity of the ground, specified as 
%      a numeric scalar in the range 1 to 100. Relative permittivity is
%      expressed as a ratio of absolute material permittivity to the
%      permittivity of vacuum. This value is used to calculate path loss
%      due to ground reflection. The default value is 15, which corresponds
%      to average ground.
%
%   AtmosphericRefractivity - Atmospheric refractivity near ground in 
%      N-units, specified as a numeric scalar in the range 250 to 400. This
%      value is used to calculate path loss due to refraction through the
%      atmosphere and tropospheric scatter. The default value is 301, which
%      corresponds to average atmospheric conditions.
%
%   Humidity - Absolute air humidity near ground in g/m^3, specified as a 
%      numeric scalar in the range 0 to 110. This value is used to
%      calculate path loss due to atmospheric absorption. The default value
%      is 9, which corresponds to the absolute humidity of air at 15
%      degrees Celsius and 70% relative humidity.
%
%   % Example: Calculate path loss with terrain obstruction
%
%   % Calculate path loss over flat terrain. Define the terrain profile for
%   % distances up to 10 km with step size of 100 m.
%   freq = 28e9;
%   R = 0:100:10000;
%   Z = zeros(1,numel(R));
%   Lterrain1 = tirempl(R,Z,freq,...
%      'TransmitterAntennaHeight',5, ...
%      'ReceiverAntennaHeight',5)
%
%   % Calculate path loss assuming free-space propagation. The result is
%   % similar to the path loss using TIREM for line-of-sight conditions.
%   c = physconst('lightspeed');
%   Lfs = fspl(R(end),c./freq)
%
%   % Introduce a terrain obstruction and recompute path loss
%   Z(40:60) = [0:1:10 9:-1:0];
%   Lterrain2 = tirempl(R,Z,freq,...
%      'TransmitterAntennaHeight',5, ...
%      'ReceiverAntennaHeight',5)
%
% See also: tiremSetup, fspl

%   Copyright 2018 The MathWorks, Inc.

% Validate required inputs
[ds,Z] = tirem.internal.Validators.validateTerrainProfile(ds,Z);
fqs = tirem.internal.Validators.validateFrequency(fqs);

% Process optional Name-Value pairs
p = tirem.internal.Validators.InputParser;
p.parse(varargin{:});
inputs = p.Results;

% Validate and get Name-Value pairs
txHeight = tirem.internal.Validators.validateTransmitterAntennaHeight(inputs.TransmitterAntennaHeight);
rxHeight = tirem.internal.Validators.validateReceiverAntennaHeight(inputs.ReceiverAntennaHeight);
refractivity = tirem.internal.Validators.validateAtmosphericRefractivity(inputs.AtmosphericRefractivity);
pol = tirem.internal.Validators.validateAntennaPolarization(inputs.AntennaPolarization);
conductivity = tirem.internal.Validators.validateGroundConductivity(inputs.GroundConductivity);
permittivity = tirem.internal.Validators.validateGroundPermittivity(inputs.GroundPermittivity);
humidity = tirem.internal.Validators.validateHumidity(inputs.Humidity);

% Perform TIREM analysis for each frequency
analyzer = tirem.internal.PropagationAnalyzer;
numFqs = numel(fqs);
pl = zeros(1,numFqs);
returnOutputs = (nargout > 1);
for fqInd = 1:numFqs
    analyzer.analyze(txHeight, rxHeight, fqs(fqInd), Z, ds, false, ...
        refractivity, conductivity, permittivity, humidity, pol, 0, 0, 0);
    pl(fqInd) = analyzer.TotalLoss;
    
    % Package TIREM analysis outputs if requested
    if returnOutputs
        output(fqInd) = analyzer.Outputs; %#ok<AGROW>
    end
end
