basin extension replaced

@miSim/lesserSinkNeighbor.m

@@ -1,11 +1,11 @@
 function obj = lesserSinkNeighbor(obj)
 % LESSERSINKNEIGHBOR  Capacity-Aware Lesser Sink Neighbor algorithm (CALSN).
 %
-% Implements the algorithm specified in new/new_routing_algorithm.pdf as the
+% Implements the algorithm given in the CALSN algorithm description document
-% selectable alternative to lesserNeighbor. UAVs 1 and N are base stations;
+% (ADD; kept outside the repository) as the selectable alternative to
-% every other UAV generates r = 1 unit of data per step, and all data must
+% lesserNeighbor. UAVs 1 and N are base stations; every other UAV generates
-% reach a base station. The algorithm outputs both a sparse maintained
+% r = 1 unit of data per step, and all data must reach a base station. The
-% topology and a routing table.
+% algorithm outputs both a sparse maintained topology and a routing table.
 %
 % 1. Base-station potentials phi via distributed Bellman-Ford over the
 %    current communication neighborhoods (obj.adjacency), with link cost
@@ -14,29 +14,31 @@ function obj = lesserSinkNeighbor(obj)
 % 2. Lesser sink neighbors L_i replace LNA's "j < i" rule:
 %        L_i = { j in N_i : phi_j < phi_i, or phi_j == phi_i and j < i }.
 %    The equal-potential lexicographic tie-break is a documented deviation
-%    from the spec's strict inequality: at an exact potential watershed
+%    from the ADD's strict inequality. It makes (phi, index) a total order,
-%    (e.g. a symmetric chain) the strict rule would split A_control into
+%    so every feasible link between two non-base UAVs has exactly one
-%    one tree per base station; the tie-break keeps the airborne graph
+%    deciding endpoint (the uphill one).
-%    connected, which the mission requires.
 % 3. LNA-style sparsity: one parent per connected component of the
 %    L_i-induced subgraph (only links with positive capacity in at least
 %    one direction count as edges), chosen by
 %        p = argmin_{j in C} ( phi_j + 1/(C_ij + eps) ).
 %    A_control(i,p) = A_control(p,i) = 1 (diagonal = I).
-%    BASIN EXTENSION (documented deviation): with two potential minima
+% 4. Connectivity bridge (documented deviation; not in the ADD). Every UAV
-%    (the two bases), one-parent-per-component alone yields a spanning
+%    must stay connected to the whole fleet to exchange positions for
-%    forest with up to two disjoint base-rooted trees whenever a UAV's
+%    guidance, so A_control must be one connected graph. LNA's grouping
-%    lesser sink neighbors collapse into a single component (dense usable
+%    argument guarantees that whenever every feasible link has a deciding
-%    graphs). To restore LNA's global-connectivity property, each UAV
+%    endpoint. The one exception is the base-to-base link: both bases have
-%    selects one representative per component AND per base-station BASIN
+%    phi = 0 and select no parents, so nobody maintains it, yet it still
-%    (the base a neighbor's potential descent terminates at). Interior
+%    merges groups when other UAVs form their components. Steps 1-3
-%    UAVs are unaffected (their whole neighborhood shares one basin);
+%    therefore leave A_control either connected or split into exactly two
-%    watershed UAVs whose neighborhoods span both basins gain a second
+%    pieces, one containing each base. In the split case one bridge is
-%    parent, which is exactly the bridge between the two trees. A final
+%    added: of all feasible links between the two pieces (base-to-base
-%    stitch guard covers degenerate geometries: if A_control still has
+%    included), the one whose weaker direction has the highest SINR, i.e.
-%    multiple components, the lowest-cost feasible inter-component link
+%    the most margin under the comms barrier. A warning is issued if even
-%    is added (maintained for connectivity only; carries no flow).
+%    that link is inside the CBF maintenance band (sinrThreshold + 3 dB).
-% 4. Routing, children first (decreasing phi; equal phi processed in
+%    The bridge is maintained for connectivity only and carries no flow,
+%    so when present it is the one maintained link outside the ADD's
+%    invariant A_control(i,j) = 1, i ~= j, iff F(i,j) > 0 or F(j,i) > 0.
+% 5. Routing, children first (decreasing phi; equal phi processed in
 %    decreasing index order so tie-break children precede their parents):
 %    aggregate demand D_i = r + inflows, split capacity-proportionally
 %    over parents; saturated links are capped at capacity and the residual
@@ -72,13 +74,16 @@ function obj = lesserSinkNeighbor(obj)
         positions(kk, :) = obj.agents{kk}.pos;
     end
 
-    % Directed Shannon link capacities over the current communication
+    % Directed SINRs and Shannon link capacities over the current
-    % neighborhoods: C(i, j) = capacity of sending from i to j (receiver j)
+    % communication neighborhoods: S(i, j) and C(i, j) describe sending
+    % from i to j (receiver j)
+    S = zeros(nAgents);
     C = zeros(nAgents);
     for ii = 1:nAgents
         for jj = 1:nAgents
             if jj ~= ii && obj.adjacency(ii, jj)
-                C(ii, jj) = obj.bandwidth * log2(1 + obj.sinrLink(positions, jj, ii));
+                S(ii, jj) = obj.sinrLink(positions, jj, ii);
+                C(ii, jj) = obj.bandwidth * log2(1 + S(ii, jj));
             end
         end
     end
@@ -109,42 +114,7 @@ function obj = lesserSinkNeighbor(obj)
         end
     end
 
-    % ---- Base-station basins (for the basin extension) -------------------
+    % ---- Parent selection (one per lesser-sink component) ---------------
-    % Each finite-potential UAV descends its argmin neighbor chain; phi is
-    % strictly decreasing along it, so the chain terminates at a base. The
-    % basin label records which base that is.
-    descent = zeros(nAgents, 1);
-    for ii = 1:nAgents
-        if isBase(ii) || isinf(phi(ii))
-            continue;
-        end
-        bestMetric = inf;
-        for jj = 1:nAgents
-            if jj == ii || ~obj.adjacency(ii, jj)
-                continue;
-            end
-            metric = phi(jj) + 1 / (C(ii, jj) + EPS_CAP);
-            if metric < bestMetric
-                bestMetric = metric;
-                descent(ii) = jj;
-            end
-        end
-    end
-    basin = zeros(nAgents, 1);
-    basin(1) = 1;
-    basin(nAgents) = nAgents;
-    for ii = 1:nAgents
-        if isBase(ii) || isinf(phi(ii))
-            continue;
-        end
-        kk = ii;
-        while ~isBase(kk)
-            kk = descent(kk);
-        end
-        basin(ii) = kk;
-    end
-
-    % ---- Parent selection (one per lesser-sink component and basin) ------
     Acontrol = logical(eye(nAgents));
     isParent = false(nAgents); % isParent(i, p): p is a selected parent of i
     for ii = 1:nAgents
@@ -187,8 +157,7 @@ function obj = lesserSinkNeighbor(obj)
             end
         end
 
-        % One parent per connected component and per base-station basin
+        % One parent per connected component, selected by capacity and
-        % present in it (basin extension), selected by capacity and
         % potential: p = argmin( phi_j + 1/(C_ij + eps) )
         visited = false(1, lsnCount);
         for aa = 1:lsnCount
@@ -198,72 +167,58 @@ function obj = lesserSinkNeighbor(obj)
             component = bfs(subgraphAdjacency, aa);
             visited(component) = true;
 
-            for bb = [1, nAgents] % the two possible basin labels
+            bestMetric = inf;
-                bestMetric = inf;
+            bestParent = 0;
-                bestParent = 0;
+            for cc = 1:numel(component)
-                for cc = 1:numel(component)
+                jj = lsn(component(cc));
-                    jj = lsn(component(cc));
+                metric = phi(jj) + 1 / (C(ii, jj) + EPS_CAP);
-                    if basin(jj) ~= bb
+                if metric < bestMetric
-                        continue;
+                    bestMetric = metric;
-                    end
+                    bestParent = jj;
-                    metric = phi(jj) + 1 / (C(ii, jj) + EPS_CAP);
-                    if metric < bestMetric
-                        bestMetric = metric;
-                        bestParent = jj;
-                    end
-                end
-                if bestParent > 0
-                    isParent(ii, bestParent) = true;
-                    Acontrol(ii, bestParent) = true;
-                    Acontrol(bestParent, ii) = true;
                 end
             end
+            isParent(ii, bestParent) = true;
+            Acontrol(ii, bestParent) = true;
+            Acontrol(bestParent, ii) = true;
         end
     end
 
-    % ---- Connectivity stitch guard ---------------------------------------
+    % ---- Connectivity bridge ---------------------------------------------
-    % The basin extension bridges the two base-rooted trees through
+    % The unmaintained base-to-base link is the only way steps 1-3 can
-    % watershed UAVs whenever each base borders a routing basin. Degenerate
+    % leave A_control split, and then into exactly two pieces, one per
-    % geometries (e.g. a base whose only feasible neighbor is the other
+    % base. Bridge them with the feasible cross link whose weaker direction
-    % base) can still leave A_control split; add the lowest-cost feasible
+    % has the highest SINR. If the pieces are not joined by any feasible
-    % link between components until one remains. Stitched links are
+    % link, the feasibility graph itself is disconnected; validate()
-    % maintained for connectivity only and carry no routed flow. If the
-    % feasibility graph itself is disconnected no stitch exists; validate()
     % reports that case.
-    while true
+    pieceOfBase1 = false(1, nAgents);
-        comp = zeros(1, nAgents);
+    pieceOfBase1(bfs(Acontrol, 1)) = true;
-        nComp = 0;
+    if ~pieceOfBase1(nAgents)
-        for ii = 1:nAgents
+        bestSinr = -inf;
-            if comp(ii) == 0
-                nComp = nComp + 1;
-                comp(bfs(Acontrol, ii)) = nComp;
-            end
-        end
-        if nComp <= 1
-            break;
-        end
-        bestMetric = inf;
         bestU = 0;
         bestV = 0;
-        for uu = 1:nAgents
+        for uu = find(pieceOfBase1)
-            for vv = (uu + 1):nAgents
+            for vv = find(~pieceOfBase1)
-                if comp(uu) ~= comp(vv) && obj.adjacency(uu, vv)
+                if obj.adjacency(uu, vv)
-                    % Bidirectional link cost (finite even for infinite
+                    sinrBinding = min(S(uu, vv), S(vv, uu));
-                    % potentials, so isolated islands remain stitchable)
+                    if sinrBinding > bestSinr
-                    metric = 1 / (C(uu, vv) + EPS_CAP) + 1 / (C(vv, uu) + EPS_CAP);
+                        bestSinr = sinrBinding;
-                    if metric < bestMetric
-                        bestMetric = metric;
                         bestU = uu;
                         bestV = vv;
                     end
                 end
             end
         end
-        if bestU == 0
+        if bestU > 0
-            break; % feasibility graph disconnected: nothing to stitch with
+            Acontrol(bestU, bestV) = true;
+            Acontrol(bestV, bestU) = true;
+            % CBF maintenance band: constrainMotion holds maintained links
+            % at sinrThreshold + 3 dB (paddingDb; keep in sync)
+            if bestSinr < 10^((obj.sinrThreshold + 3.0) / 10)
+                warning("miSim:bridgeInCbfBand", ...
+                    "Connectivity bridge %d-%d has SINR %.1f dB, inside the CBF maintenance band (sinrThreshold %.1f dB + 3 dB)", ...
+                    bestU, bestV, 10 * log10(bestSinr), obj.sinrThreshold);
+            end
         end
-        Acontrol(bestU, bestV) = true;
-        Acontrol(bestV, bestU) = true;
     end
 
     % ---- Routing (children first: decreasing potential order) -----------

@miSim/miSim.m

@@ -32,10 +32,10 @@ classdef miSim
         bandwidth = 20e6;         % channel bandwidth (Hz) for thermal noise N = k_B*T*B
         % Topology-selection algorithm (requires useSinrComms; MATLAB sim only).
         % false -> lesser-neighbor; true -> Capacity-Aware Lesser Sink Neighbor
-        % (CALSN, see new/new_routing_algorithm.pdf): base-station potentials
+        % (CALSN, see the algorithm description document, ADD): base-station
-        % replace the index-based "lesser" rule, and the maintained topology
+        % potentials replace the index-based "lesser" rule, and the maintained
-        % comes with a routing table (flows F, fractions R) that sends every
+        % topology comes with a routing table (flows F, fractions R) that sends
-        % UAV's data toward base stations 1 and N.
+        % every UAV's data toward base stations 1 and N.
         useRoutingTopology = false;
         routingAdjacencyMatrix = false(0, 0); % directed flow links (i, j) = child i sends to parent j (F > 0)
         routingFlows = zeros(0, 0);      % F: flow on each directed link (units of r = 1 per UAV per step)

@miSim/sinrLink.m

@@ -60,7 +60,7 @@ function [sinr, grad] = sinrLink(obj, positions, rxIdx, txIdx)
     S      = P_t / (K * d_t^n);
 
     % Interference I = sum of received powers from every other agent at r
-    I = 0.0;
+    I = 0;
     for k = 1:nAgents
         if k == r || k == t
             continue;

@miSim/validate.m

@@ -20,7 +20,7 @@ function validate(obj)
     end
 
     % In routing (CALSN) mode the MAINTAINED topology itself must form a
-    % single connected component (basin extension + stitch guard invariant)
+    % single connected component (guaranteed by the connectivity bridge)
     if obj.useRoutingTopology && max(conncomp(graph(obj.constraintAdjacencyMatrix))) ~= 1
         error("Maintained routing topology (A_control) is not connected");
     end

@sigmoidSensor/plotParameters.m

@@ -46,7 +46,7 @@ function f = plotParameters(obj, f)
     nexttile(2, [1,  1]);
     grid("on");
     title("Tilt Membership Sigmoid");
-    xlabel("Tilt (deg)");
+    xlabel("Tilt ($^\circ$)");
     ylabel("Membership");
     lt = legend();
     hold("on");

test/test_miSim.m

@@ -1007,7 +1007,7 @@ classdef test_miSim < matlab.unittest.TestCase
             % corner, co-altitude at z = 2. Base 1 trails the fleet and
             % base 4 leads it, with the relays staggered between: relay 3
             % starts too far from base 1 for a direct link to be worth
-            % using, so its initial basin-1 route is MULTI-HOP via relay 2
+            % using, so its initial route to base 1 is MULTI-HOP via relay 2
             % and re-parents to the base directly as the fleet contracts.
             tc.agents = {agent; agent; agent; agent;};
             tc.collisionRanges = .25 * ones(size(tc.agents));
@@ -1128,7 +1128,7 @@ classdef test_miSim < matlab.unittest.TestCase
             % corner, co-altitude at z = 2. Base 1 trails the fleet and
             % base 4 leads it, with the relays staggered between: relay 3
             % starts too far from base 1 for a direct link to be worth
-            % using, so its initial basin-1 route is MULTI-HOP via relay 2
+            % using, so its initial route to base 1 is MULTI-HOP via relay 2
             % and re-parents to the base directly as the fleet contracts.
             tc.agents = {agent; agent; agent; agent;};
             tc.collisionRanges = .25 * ones(size(tc.agents));

test/test_sigmoidSensor.m

@@ -57,30 +57,6 @@ classdef test_sigmoidSensor < matlab.unittest.TestCase
             % Performance at great tilt should be 0
             tc.verifyEqual(tc.testClass.sensorPerformance([0, 0, h], [5, 5, 0]), 0, "AbsTol", 1e-9);
         end
-        function test_thesis_graphic(tc)
-            s = [sigmoidSensor, sigmoidSensor];
-
-            alphaDist = [2, 8];
-            betaDist = [25, 1];
-            alphaTilt = [50, 15]; % degrees
-            betaTilt = [0.2, 3];
-
-            s(1) = s(1).initialize(alphaDist(1), betaDist(1), alphaTilt(1), betaTilt(1));
-            s(2) = s(2).initialize(alphaDist(2), betaDist(2), alphaTilt(2), betaTilt(2));
-
-            set(groot, 'defaultTextInterpreter',          'latex');
-            set(groot, 'defaultAxesTickLabelInterpreter', 'latex');
-            set(groot, 'defaultLegendInterpreter',        'latex');
-            set(groot, 'defaultColorbarTickLabelInterpreter', 'latex');
-            set(groot, 'defaultAxesFontSize', 10);   % = \footnotesize in a 12pt doc
-            set(groot, 'defaultTextFontSize', 10);
-
-            % Plot
-            f = s(1).plotParameters;
-            f = s(2).plotParameters(f);
-
-            % f = exportThesisFig(f, "sigmoidmodel", 3/4);
-        end
     end
 
 end

test_thesis_plots.m

@@ -0,0 +1,230 @@
+classdef test_thesis_plots < matlab.unittest.TestCase
+    % test class that generates plots I use in my thesis
+    methods (Test)
+        function service_performance(tc)
+            % TX parameters
+            P_TX = 10e-3;
+            T_0 = 300;
+            BW = 20e6;
+            f_c = 2e9;
+            G_RX_dBi = 3;
+            beamwidthExponent = [16, 12, 32];
+            elevation = [30, 20, 45];
+            azimuth = [15, 325, 200];
+            lossExponent = 2;
+
+            % TX Positions
+            pos1 = [50,  50,  35];
+            pos2 = [70, 30,  30];
+            pos3 = [40, 75,  40];
+
+            % User Distribution Parameters
+            mu = [78, 66];
+            sig = reshape([2000, 1100; 1100, 2000], [1, 2, 2]);
+
+            sensor1 = rfSensor;
+            sensor1 = sensor1.initialize(P_TX, T_0, BW, f_c, G_RX_dBi, beamwidthExponent(1), elevation(1), azimuth(1), lossExponent);
+            sensor2 = rfSensor;
+            sensor2 = sensor2.initialize(P_TX, T_0, BW, f_c, G_RX_dBi, beamwidthExponent(2), elevation(2), azimuth(2), lossExponent);
+            sensor3 = rfSensor;
+            sensor3 = sensor3.initialize(P_TX, T_0, BW, f_c, G_RX_dBi, beamwidthExponent(3), elevation(3), azimuth(3), lossExponent);
+
+            % Compute sensor and system performance for a single user distribution
+            domain = rectangularPrism;
+            domain = domain.initialize([zeros(1, 3); 100 * ones(1, 3)], "DOMAIN", "Domain");
+
+            objectiveFunction = objectiveFunctionWrapper(mu, sig);
+            domain.objective = domain.objective.initialize(objectiveFunction, domain, 0.1, 0, 1e-6, mu, sig);
+
+            % Build target grid matching the domain objective grid and flatten user distribution
+            nX = size(domain.objective.values, 2);
+            nY = size(domain.objective.values, 1);
+            x_domain = linspace(0, 100, nX);
+            y_domain = linspace(0, 100, nY);
+            [Xd, Yd] = meshgrid(x_domain, y_domain);
+            targetPos_domain = [Xd(:), Yd(:), zeros(numel(Xd), 1)];
+            userDist = domain.objective.values(:) / sum(domain.objective.values, "all");
+
+            [SINR1, ~, sensor1, others] = sensor1.sensorPerformance(pos1, targetPos_domain, [pos2; pos3], {sensor2; sensor3});
+            sensor2 = others{1}; sensor3 = others{2};
+            [SINR2, ~, sensor2, others] = sensor2.sensorPerformance(pos2, targetPos_domain, [pos1; pos3], {sensor1; sensor3});
+            sensor1 = others{1}; sensor3 = others{2};
+            [SINR3, ~, sensor3, ~]      = sensor3.sensorPerformance(pos3, targetPos_domain, [pos1; pos2], {sensor1; sensor2});
+
+            maxSINR_linear = max([10.^(SINR1/10), 10.^(SINR2/10), 10.^(SINR3/10)], [], 2);
+
+            % Plot SINR from each UAV's perspective with a shared color scale,
+            % one top-down map per UAV, in the style of the sigmoid sensor figure
+            sharedClim = [-80, max([SINR1; SINR2; SINR3])];
+            SINR = {SINR1, SINR2, SINR3};
+            sensors = {sensor1, sensor2, sensor3};
+            positions = [pos1; pos2; pos3];
+            for ii = 1:3
+                f = plotSinrMap(x_domain, y_domain, reshape(SINR{ii}, size(Xd)), ...
+                    sensors, positions, ii, sharedClim);
+                if ii == 3 % only plot for #3 right now
+                    f = exportThesisFig(f, "sinrexample" + ii, 0.95, Colormap="parula", Recolor=false, Resolution=300);
+                end
+                close(f);
+            end
+        end
+        function sigmoid_sensor_performance(tc)
+            s = sigmoidSensor;
+            s = s.initialize(65, 1, 50, 0.2, 30, 120);
+            p = 50 * ones(1, 3);
+            mu = [50, 50];
+            sig = 25 * reshape(eye(2), [1, 2, 2]);
+            domain = rectangularPrism;
+            domain = domain.initialize([zeros(1, 3); 100 * ones(1, 3)], "DOMAIN", "Domain");
+            objectiveFunction = objectiveFunctionWrapper(mu, sig);
+            domain.objective = domain.objective.initialize(objectiveFunction, domain, 0.1, 0, 1e-6, mu, sig);
+
+            S = s.sensorPerformance(p, [domain.objective.X(:), domain.objective.Y(:), zeros(size(domain.objective.X(:)))]);
+            S = reshape(S, size(domain.objective.X));
+
+            % meshgrid puts y down the rows, so plot against the grid vectors
+            % and put y = 0 at the bottom (imagesc defaults to row 1 at the top)
+            x = domain.objective.X(1, :);
+            y = domain.objective.Y(:, 1);
+
+            f = figure;
+            imagesc(x, y, S);
+            set(gca, "YDir", "normal");
+            axis("image");
+            hold("on");
+            % Boresight direction: the arrow runs from the UAV to where the
+            % boresight meets the ground, so its length grows with tilt
+            % (altitude * tan(tilt)). Azimuth 0 = +Y, 90 = +X, as in sensorPerformance.
+            reach = p(3) * tand(s.tilt);
+            hit = p(1:2) + reach * [sind(s.azimuth), cosd(s.azimuth)];
+            quiver(p(1), p(2), hit(1) - p(1), hit(2) - p(2), 0, "k", "MaxHeadSize", 0.4);
+            hUav = plot(p(1), p(2), "k+");
+            hHit = plot(hit(1), hit(2), "kx");
+            hold("off");
+            legend([hUav, hHit], "UAV position", "Boresight-ground intersection", "Location", "northwest");
+            xlabel("$x$ (m)");
+            ylabel("$y$ (m)");
+            cb = colorbar;
+            cb.Label.String = "Sensor performance $s_n$";
+            clim([0, 1]);
+
+            % Full text width so it's included at 1:1 like the other figures;
+            % 0.8 fits the square axes plus colorbar without side margins
+            f = exportThesisFig(f, "sigmoidsensorexample", 0.8, Colormap="parula", Recolor=false, Resolution=300);
+            close(f);
+        end
+        function antenna_patterns(tc)
+            f = figure;
+            xi = -90:0.1:90; % degrees
+            n = [1, 4, 16, 64, 256];
+            for ii = 1:length(n)
+                gain = 10 .* log10(2 * (n(ii) + 1)) + 10 .* n(ii) .* log10(max(0, cosd(xi)));
+                plot(xi, gain, "LineWidth", 2);
+                if ii == 1
+                    hold("on");
+                end
+            end
+            hold("off");
+            grid("on");
+            ylim([-100, 35]);
+            xlim([-90, 90]);
+            xticks(-90:45:90);
+            % exportThesisFig switches all text to the LaTeX interpreter
+            title("Transmitting Antenna Gain Model");
+            ylabel("Transmitting Antenna Gain (dBi)");
+            xlabel("Elevation Angle from Antenna Boresight $\xi$ ($^\circ$)");
+            legend("$\delta = " + string(n) + "$");
+
+            exportThesisFig(f, "antennapatterns", 3/4);
+            close(f);
+        end
+        function sigmoid_sensor_membership_functions(tc)
+            s = [sigmoidSensor, sigmoidSensor];
+
+            alphaDist = [2, 8];
+            betaDist = [25, 1];
+            alphaTilt = [50, 15]; % degrees
+            betaTilt = [0.2, 3];
+
+            s(1) = s(1).initialize(alphaDist(1), betaDist(1), alphaTilt(1), betaTilt(1));
+            s(2) = s(2).initialize(alphaDist(2), betaDist(2), alphaTilt(2), betaTilt(2));
+
+            set(groot, 'defaultTextInterpreter',          'latex');
+            set(groot, 'defaultAxesTickLabelInterpreter', 'latex');
+            set(groot, 'defaultLegendInterpreter',        'latex');
+            set(groot, 'defaultColorbarTickLabelInterpreter', 'latex');
+            set(groot, 'defaultAxesFontSize', 10);   % = \footnotesize in a 12pt doc
+            set(groot, 'defaultTextFontSize', 10);
+
+            % Plot
+            f = s(1).plotParameters;
+            f = s(2).plotParameters(f);
+
+            f = exportThesisFig(f, "sigmoidmodel", 3/4);
+            close(f);
+        end
+    end
+
+end
+function f = plotSinrMap(x, y, SINR, sensors, positions, tx, cLimits)
+    % Top-down SINR map for transmitting UAV tx, matching the sigmoid sensor
+    % figure. Every UAV gets a marker, an arrow to the X where its boresight
+    % meets the ground, and its -3 dB (half-power) beam footprint; the
+    % transmitting UAV is drawn solid, the interferers dashed.
+    f = figure;
+    imagesc(x, y, SINR);
+    set(gca, "YDir", "normal");   % y = 0 at the bottom
+    axis("image");
+    hold("on");
+    for ii = 1:numel(sensors)
+        [hit, ring] = beamFootprint(sensors{ii}, positions(ii, :));
+        pos = positions(ii, :);
+        if ii == tx
+            style = "-";  marker = "k+";
+        else
+            style = "--"; marker = "ko";
+        end
+        hRing(ii) = plot(ring(:, 1), ring(:, 2), "k" + style); %#ok<AGROW>
+        quiver(pos(1), pos(2), hit(1) - pos(1), hit(2) - pos(2), 0, "k", ...
+            "LineStyle", style, "MaxHeadSize", 0.4);
+        hHit = plot(hit(1), hit(2), "kx");
+        hUav(ii) = plot(pos(1), pos(2), marker, "MarkerFaceColor", "w"); %#ok<AGROW>
+    end
+    hold("off");
+    rx = setdiff(1:numel(sensors), tx);
+    legend([hUav(tx), hUav(rx(1)), hHit, hRing(tx), hRing(rx(1))], ...
+        "Transmitting UAV", "Interfering UAVs", "Boresight-ground intersection", ...
+        "Transmitter $-3$ dB contour", "Interferer $-3$ dB contour", ...
+        "Location", "southoutside", "NumColumns", 2);
+    xlim([min(x), max(x)]);   % footprints can extend past the area
+    ylim([min(y), max(y)]);
+    xlabel("$x$ (m)");
+    ylabel("$y$ (m)");
+    cb = colorbar;
+    cb.Label.String = "SINR (dB)";
+    clim(cLimits);
+end
+
+function [hit, ring] = beamFootprint(sensor, pos)
+    % Where the boresight meets the ground (z = 0), and the ground curve where
+    % the gain is 3 dB below boresight: the intersection of the half-power cone
+    % with the ground, as in rfSensor/plot. Azimuth 0 = +Y, 90 = +X.
+    tlt = sensor.tilt;
+    az = sensor.azimuth;
+    hit = pos(1:2) + pos(3) * tand(tlt) * [sind(az), cosd(az)];
+
+    % Rotate nadir to the boresight
+    Ry = [cosd(tlt), 0, -sind(tlt); 0, 1, 0; sind(tlt), 0, cosd(tlt)];
+    Rz = [sind(az), -cosd(az), 0; cosd(az), sind(az), 0; 0, 0, 1];
+    R = Rz * Ry;
+    b = R * [0; 0; -1];
+    u = R * [1; 0; 0];
+    v = R * [0; 1; 0];
+
+    ha = sensor.halfAngle();
+    phi = linspace(0, 2*pi, 720)';
+    dirs = cosd(ha) .* b' + sind(ha) .* (cos(phi) .* u' + sin(phi) .* v');
+    t = -pos(3) ./ dirs(:, 3);
+    t(t <= 0) = NaN;   % cone generators that never reach the ground
+    ring = pos(1:2) + t .* dirs(:, 1:2);
+end

util/exportThesisFig.m

@@ -5,8 +5,12 @@ function f = exportThesisFig(f, name, aspect, opts)
     %   exportThesisFig(f, "sigmoidmodel")                 % full width, 4:3
     %   exportThesisFig(f, "traj", 1, Width=230)           % half-width subfigure
     %   exportThesisFig(f, "traj", 3/4, LineStyles=true)   % distinguishable in grayscale
+    %   exportThesisFig(f, "perf", 0.8, Width=330)         % imagesc plot -> name.png
     %
-    % The figure is modified in place.
+    % Figures containing images (imagesc, image) are exported as a lossless
+    % 400 dpi PNG instead: MATLAB's vector PDF export draws them with a
+    % visible diagonal seam. The PNG records its dpi, so LaTeX still places it
+    % at Width points. The figure is modified in place.
     arguments (Input)
         f (1, 1) {mustBeA(f, "matlab.ui.Figure")};
         name (1, 1) string;
@@ -15,11 +19,17 @@ function f = exportThesisFig(f, name, aspect, opts)
         opts.FontSize (1, 1) double {mustBePositive} = 10;    % \footnotesize in a 12pt doc
         opts.LineWidth (1, 1) double {mustBePositive} = 2; % data lines, points
         opts.AxesLineWidth (1, 1) double {mustBePositive} = 0.75;
-        opts.MarkerSize (1, 1) double {mustBePositive} = 4;
+        opts.MarkerSize (1, 1) double {mustBePositive} = 4;        % markers on lines
+        opts.PointMarkerSize (1, 1) double {mustBePositive} = 16;  % marker-only series, e.g. plot(x, y, "k+")
+        opts.PointMarkerLineWidth (1, 1) double {mustBePositive} = 2.5;
         opts.Palette (:, 3) double = okabeIto();              % colorblind- and print-safe
         opts.Recolor (1, 1) logical = true;                   % also override explicitly set colors
         opts.LineStyles (1, 1) logical = false;               % cycle line styles per series
-        opts.Raster (1, 1) logical = false;                   % for very dense plots
+        opts.Raster (1, 1) logical = ~isempty(findall(f, 'Type', 'image')); % PNG; also for very dense plots
+        % dpi, for Raster. MATLAB silently writes an all-black PNG when the image is
+        % too large (3909 px wide failed, 2750 px worked), so full width stays at 400
+        opts.Resolution (1, 1) double {mustBePositive} = 400;
+        opts.Colormap = [];                                   % for images/surfaces, e.g. "parula"; [] keeps the figure's
         opts.OutDir (1, 1) string = matlab.project.rootProject().RootFolder;
     end
     arguments (Output)
@@ -41,6 +51,17 @@ function f = exportThesisFig(f, name, aspect, opts)
         if opts.LineStyles
             ax.LineStyleOrder = styles;
         end
+        if ~isempty(opts.Colormap)
+            colormap(ax, opts.Colormap);
+        end
+    end
+
+    % Colorbars (image and surface plots) get the same frame as the axes. Their
+    % Label isn't reached by findall above, so set it directly.
+    for cb = findall(f, 'Type', 'colorbar')'
+        cb.LineWidth = opts.AxesLineWidth;
+        cb.Label.Interpreter = 'latex';
+        cb.Label.FontSize = opts.FontSize;
     end
 
     % Data series. SeriesIndex is the series' position in its axes' color
@@ -57,11 +78,21 @@ function f = exportThesisFig(f, name, aspect, opts)
         isLine = isprop(s, 'Color');
         isScatter = isa(s, 'matlab.graphics.chart.primitive.Scatter');
 
-        if isLine || isScatter
+        % Marker-only series (no connecting line) mark individual points; a
-            s.LineWidth = opts.LineWidth;
+        % thick data-line width would fill in markers like + and x
-        end
+        isPoint = isLine && isprop(s, 'Marker') && strcmp(s.LineStyle, 'none') ...
-        if isprop(s, 'MarkerSize')
+            && ~strcmp(s.Marker, 'none');
-            s.MarkerSize = opts.MarkerSize;
+
+        if isPoint
+            s.LineWidth = opts.PointMarkerLineWidth;
+            s.MarkerSize = opts.PointMarkerSize;
+        else
+            if isLine || isScatter
+                s.LineWidth = opts.LineWidth;
+            end
+            if isprop(s, 'MarkerSize')
+                s.MarkerSize = opts.MarkerSize;
+            end
         end
         if opts.LineStyles && isLine && ~strcmp(s.LineStyle, 'none')
             s.LineStyle = styles(mod(k - 1, numel(styles)) + 1);
@@ -89,9 +120,9 @@ function f = exportThesisFig(f, name, aspect, opts)
     f.PaperSize = [w, h];
 
     if opts.Raster
-        renderer = {'-image', '-r300'};
+        % PNG rather than '-dpdf -image', which embeds a lossy JPEG
+        print(f, fullfile(opts.OutDir, name + ".png"), '-dpng', "-r" + opts.Resolution);
     else
-        renderer = {'-vector'};
+        print(f, fullfile(opts.OutDir, name + ".pdf"), '-dpdf', '-vector');
     end
-    print(f, fullfile(opts.OutDir, name + ".pdf"), '-dpdf', renderer{:});
 end