fq = 28e9; % 28 GHz

tx = txsite("Name","South Uncanoonuc (BS)", ...
    "Latitude",42.983723, ...
    "Longitude",-71.587173, ...
    "TransmitterPower",1, ...
    "TransmitterFrequency",fq);
show(tx)

rxBedford = rxsite("Name","Bedford Town Center", ...
    "Latitude",42.946193, ...
    "Longitude",-71.516234);

rxStA = rxsite("Name","St. Anselm College", ...
    "Latitude",42.987386, ...
    "Longitude",-71.507475);

rxGPD = rxsite("Name","Goffstown Police Dept", ...
    "Latitude",43.009335, ...
    "Longitude",-71.539083);

rxs = [rxBedford, rxStA, rxGPD];
show(rxs)

los(tx,rxs)

% Place antennas on structures at receiver sites. Assume 6 m utility poles for Bedford
% and St. Anselm sites, and 15 m antenna pole at Goffstown Police Department.
rxBedford.AntennaHeight = 6;
rxStA.AntennaHeight = 6;
rxGPD.AntennaHeight = 15;

% Increase height of antenna at base station until line-of-sight is achieved with all receiver sites
tx.AntennaHeight = 10;
while ~all(los(tx,rxs))
    tx.AntennaHeight = tx.AntennaHeight + 5;
end

% Display line-of-sight
los(tx,rxs)
disp("Antenna height required for line-of-sight: " + tx.AntennaHeight + " m")

% Design reflector-backed crossed dipole antenna
txElement = reflectorCrossedDipoleElement(fq);

% Define array size
ntxrow = 8;
ntxcol = 12;

% Define element spacing
lambda = physconst("lightspeed")/fq;
drow = lambda/2;
dcol = lambda/2;

% Create 8-by-12 antenna array
tx.Antenna = phased.URA("Size",[ntxrow ntxcol], ...
    "Element",txElement, ...
    "ElementSpacing",[drow dcol]);

% Plot pattern on the map
pattern(tx)

rxElement = reflectorDipoleElement(fq);

% Define array size
nrxrow = 3;
nrxcol = 3;
    
% Define element spacing
lambda = physconst("lightspeed")/fq;
drow = lambda/2;
dcol = lambda/2;

% Create antenna array
rxarray = phased.URA("Size",[nrxrow nrxcol], ...
    "Element",rxElement, ...
    "ElementSpacing",[drow dcol]);

% Assign array to each receiver site and point toward base station
for rx = rxs
    rx.Antenna = rxarray;
    rx.AntennaAngle = angle(rx, tx);
    pattern(rx,fq)
end

steeringVector = phased.SteeringVector("SensorArray",tx.Antenna);
for rx = rxs
    % Compute steering vector for receiver site
    [az,el] = angle(tx,rx);
    sv = steeringVector(fq,[az;el]);
    
    % Update base station radiation pattern
    tx.Antenna.Taper = conj(sv);
    pattern(tx)
    
    % Compute signal strength (dBm)
    ss = sigstrength(rx,tx,"freespace");
    disp("Signal strength at " + rx.Name + ":")
    disp(ss + " dBm")
end

steeringVector = phased.SteeringVector("SensorArray",tx.Antenna);

% Compute steering vector for receiver site
[az,el] = angle(tx,rxs);
sv = steeringVector(fq,[az el]');

% Update base station radiation pattern
tx.Antenna.Taper = conj(sum(sv,2));
pattern(tx)

% Compute signal strength (dBm)
for rx = rxs
    ss = sigstrength(rx,tx,"freespace");
    disp("Signal strength at " + rx.Name + ":")
    disp(ss + " dBm")
end

% Assume that propagation path travels through 25 m of foliage
foliageDepth = 25;
L = 1.33*((fq/1e9)^0.284)*foliageDepth^0.588; % Weissberger model for d > 14
disp("Path loss due to foliage: " + L + " dB")

% Assign foliage loss as static SystemLoss on each receiver site
for rx = rxs
    rx.SystemLoss = L;
end

% Compute signal strength with foliage loss
for rx = rxs
    rx.SystemLoss = L;
    ss = sigstrength(rx,tx,"freespace");
    disp("Signal strength at " + rx.Name + ":")
    disp(ss + " dBm")
end

% Compute signal strength including propagation through gas and rain. Use
% the "+" operator to add the propagation models to create a composite
% model including both atmospheric effects.
weatherpm = propagationModel("gas") + propagationModel("rain");
for rx = rxs
    ss = sigstrength(rx,tx,weatherpm);
    disp("Signal strength at " + rx.Name + ":")
    disp(ss + " dBm")
end

fq = 3.5e9; % 3.5 GHz

% Create antenna array for base station
lambda = physconst("lightspeed")/fq;
drow = lambda/2;
dcol = lambda/2;
tx.TransmitterFrequency = fq;
tx.Antenna = phased.URA("Size",[ntxrow ntxcol], ...
    "Element",reflectorCrossedDipoleElement(fq), ...
    "ElementSpacing",[drow dcol]);

% Create antenna array for receiver sites
lambda = physconst("lightspeed")/fq;
drow = lambda/2;
dcol = lambda/2;
rxarray = phased.URA("Size",[nrxrow nrxcol], ...
    "Element",reflectorDipoleElement(fq), ...
    "ElementSpacing",[drow dcol], ...
    "ArrayNormal","x");
for rx = rxs
    rx.Antenna = rxarray;
end

% Compute steering vector for receiver site
steeringVector = phased.SteeringVector("SensorArray",tx.Antenna);
[az,el] = angle(tx,rxs);
sv = steeringVector(fq,[az el]');

% Update base station radiation pattern
tx.Antenna.Taper = conj(sum(sv,2));
pattern(tx,'Size',4000)
    
% Recompute loss due to foliage
L = 1.33*((fq/1e9)^0.284)*foliageDepth^0.588; % Weissberger model for d > 14

% Assign foliage loss as static SystemLoss on each receiver site
for rx = rxs
    rx.SystemLoss = L;
end
disp("Path loss due to foliage: " + L + " dB")

% Add weather-based path loss to the Longley-Rice propagation model
pm = propagationModel('longley-rice') + weatherpm;

% Compute receiver gain from peak antenna gain and system loss
G = pattern(rxarray, fq);
rxGain = max(G(:)) - L;

coverage(tx, ...
    'PropagationModel',pm, ...
    'ReceiverGain',rxGain, ...
    'ReceiverAntennaHeight',6, ...
    'SignalStrengths',-84:-50)

% Compute signal strength with foliage loss and weather
for rx = rxs
    ss = sigstrength(rx,tx,pm);
    disp("Signal strength at " + rx.Name + ":")
    disp(ss + " dBm")
end

function element = reflectorCrossedDipoleElement(fq, showAntenna)
%reflectorCrossedDipoleElement   Design reflector-backed crossed dipole antenna element

if nargin < 2
    showAntenna = false;
end

lambda = physconst("lightspeed")/fq;
offset = lambda/50;
gndspacing = lambda/4;
gndLength = lambda;
gndWidth = lambda;

% Design crossed dipole elements
d1 = design(dipole,fq);
d1.Tilt = [90,-45];
d1.TiltAxis = ["y","z"];
d2 = copy(d1);
d2.Tilt = 45;
d2.TiltAxis = "x";

% Design reflector
r = design(reflector,fq);
r.Exciter = d1;
r.GroundPlaneLength = gndLength;
r.GroundPlaneWidth = gndWidth;
r.Spacing = gndspacing;
r.Tilt = 90;
r.TiltAxis = "y";
if showAntenna
    show(r)
end

% Form the crossed dipole backed by reflector
refarray = conformalArray;
refarray.ElementPosition(1,:) = [gndspacing 0 0];
refarray.ElementPosition(2,:) = [gndspacing+offset 0 0];
refarray.Element = {r, d2};
refarray.Reference = "feed";
refarray.PhaseShift = [0 90];
if showAntenna
    show(refarray);
    view(65,20)
end

% Create custom antenna element from pattern
[g,az,el] = pattern(refarray,fq);
element = phased.CustomAntennaElement;
element.AzimuthAngles = az;
element.ElevationAngles = el;
element.MagnitudePattern = g;
element.PhasePattern = zeros(size(g));
end

function element = reflectorDipoleElement(fq)
%reflectorDipoleElement   Design reflector-backed dipole antenna element

% Design reflector and exciter, which is vertical dipole by default
element = design(reflector,fq);
element.Exciter = design(element.Exciter,fq);

% Tilt antenna element to radiate in xy-plane, with boresight along x-axis
element.Tilt = 90;
element.TiltAxis = "y";
element.Exciter.Tilt = 90;
element.Exciter.TiltAxis = "y";
end