Iridium Satellite Spot Beam Coverage on the US

% Iridium satellite frequency
freqLow = 1616e6;                                       % Hz
freqHigh = 1626.5e6;                                    % Hz
freq = (freqLow + freqHigh)/2;                          % Hz

% Determine the wavelength
lambda = physconst('LightSpeed')/freq;

% Define array size in row and column
nrow = 12;
ncol = 10;

% Define row and column element spacing
drow = lambda/2;
dcol = lambda/2;

% Configure Iridium 120-element MMA
iridiumArray = phased.URA(...
    'Size',[nrow ncol],...
    'ElementSpacing',[drow dcol],...
    'Lattice','Triangular');
iridiumArray.Element.BackBaffled = true;

% View Iridium MMA
figure
viewArray(iridiumArray) 

% Plot 3D directivity pattern of Iridium MMA
figure
pattern(iridiumArray,freq)

% Define start and end time
startTime = datetime(2021,3,17,22,52,0);               % 17 March 2021 10:52 PM UTC
stopTime = startTime + minutes(1);                     % 17 March 2021 10:53 PM UTC
sampleTime = 60;                                       % s

% Configure satellite scenario
sc = satelliteScenario(startTime,stopTime,sampleTime);

% Define the semimajor axis of the orbit of the satellite
semiMajorAxis = (6371 + 780)*1e3;                       % m

% Define the shape of the orbit of the satellite
eccentricity = 0;

% Define the inclination angle
inclination = 86.4;                                     % degrees

% Define the argument of periapsis 
argOfPeriapsis = 0;                                     % degrees

% Define the right ascension of ascending node
raan = 55;                                              % degrees

% Define the true anomaly
trueAnomaly = 38;                                       % degrees

% Configure Iridium satellite
iridiumSat = satellite(sc, semiMajorAxis, ...
    eccentricity, ...
    inclination, ...
    raan, ...
    argOfPeriapsis, ...
    trueAnomaly, ...
    Name = "Iridium Satellite", ...
    OrbitPropagator = "two-body-keplerian");

% Satellite scenario
v = satelliteScenarioViewer(sc,ShowDetails=false); 

% Number of spot beams
numSpotBeam = 48;

% Define 48 spot beam indices in 4 tiers
tier1 = [16,32,48];                                     % inner-most beams
tier2 = [13:15,29:31,45:47];                            % inner-middle beams
tier3 = [8:12,24:28,40:44];                             % outer-middle beams
tier4 = [1:7,17:23,33:39];                              % outer-most beams

% Define transmission parameters
maxRange = 5e6;                                         % m
transmitterPower = 11;                                  % watts
thetaoffset = [11,19,33,45];                            % degrees

% Define signal strength range
sigstrengthlow = -105;                                  % dBm
sigstrengthhigh = -95;                                  % dBm
sigstrengths = sigstrengthlow:sigstrengthhigh; 
sigcolors = [sigstrengthlow,sigstrengthhigh];

% Configure steering vector for spot beam beamforming
steeringvector = phased.SteeringVector("SensorArray",iridiumArray);

% Obtain Iridium satellite position 
iridiumpos = states(iridiumSat,startTime,"CoordinateFrame","geographic");

% Configure site viewer to view spot beams
viewer = siteviewer('Terrain','none');

% Loop beam index to plot 48 spot beams
for beamIdx = 1:numSpotBeam
    % Configure phi/theta angles for each spot beam
    if any(ismember(beamIdx,tier1))
        k = find(beamIdx == tier1);
        phi = 360/numel(tier1)*(k-1);                   % degrees
        theta = thetaoffset(1);                         % degrees
    elseif any(ismember(beamIdx,tier2))
        k = find(beamIdx == tier2);
        phi = 360/numel(tier2)*(k-1);                   % degrees
        theta = thetaoffset(2);                         % degrees
    elseif any(ismember(beamIdx,tier3))
        k = find(beamIdx == tier3);
        phi = 360/numel(tier3)*(k-1);                   % degrees
        theta = thetaoffset(3);                         % degrees
    else
        k = find(beamIdx == tier4);
        phi = 360/numel(tier4)*(k-1);                   % degrees
        theta = thetaoffset(4);                         % degrees
    end

    % Convert steering angle from phi/theta form to azimuth/elevation form
    azel = phitheta2azel([phi;theta]);

    % Set transmit beamforming weights for a spot beam
    sv = steeringvector(freq,azel);
    iridiumArray.Taper = sv';

    % Configure the transmit site for a spot beam
    spotbeam = txsite("Name","Iridium-Next Satellite Tx",...
        "Latitude",iridiumpos(1), ...
        "Longitude",iridiumpos(2), ...
        "AntennaHeight",iridiumpos(3),...
        "Antenna",iridiumArray,...
        "AntennaAngle",[0;-90],...
        "TransmitterFrequency",freq,...
        "TransmitterPower",transmitterPower);

    % Plot 3D beam pattern in animation
    pattern(spotbeam)

    % Display coverage map
    coverage(spotbeam,"freespace",...
        'MaxRange',maxRange,...
        'SignalStrengths',sigstrengths,...
        'ColorLimits',sigcolors)
end