From ded6f2c95c4996209f13b201d6160c93b947aa66 Mon Sep 17 00:00:00 2001 From: rohithreddyvennam Date: Wed, 27 May 2026 12:59:23 -0700 Subject: [PATCH 01/10] Add simulator scripts and beampattern figures (INFOCOM 2026) --- .github/workflows/ci.yml | 28 +++ .gitignore | 27 +++ Readme.md | 121 +++++++++- arraylink/__init__.py | 33 +++ arraylink/array_geometry.py | 186 +++++++++++++++ arraylink/beamforming.py | 307 ++++++++++++++++++++++++ arraylink/channel.py | 210 +++++++++++++++++ arraylink/utils.py | 77 ++++++ configs/arraylink_1km.yaml | 56 +++++ configs/upa_baseline.yaml | 48 ++++ data/README.md | 45 ++++ environment.yml | 19 ++ scripts/fig04_gain_vs_arrays.py | 104 ++++++++ scripts/fig06_mimo_boundaries.py | 138 +++++++++++ scripts/fig09_beampattern_sim.py | 339 +++++++++++++++++++++++++++ scripts/fig10_hardware_validation.py | 182 ++++++++++++++ scripts/fig11_2d_beampattern.py | 295 +++++++++++++++++++++++ scripts/fig12_mimo_dof.py | 158 +++++++++++++ scripts/fig13_throughput.py | 131 +++++++++++ tests/__init__.py | 0 tests/test_array_geometry.py | 122 ++++++++++ tests/test_beamforming.py | 117 +++++++++ tests/test_channel.py | 128 ++++++++++ tests/test_smoke.py | 53 +++++ 24 files changed, 2922 insertions(+), 2 deletions(-) create mode 100644 .github/workflows/ci.yml create mode 100644 .gitignore create mode 100644 arraylink/__init__.py create mode 100644 arraylink/array_geometry.py create mode 100644 arraylink/beamforming.py create mode 100644 arraylink/channel.py create mode 100644 arraylink/utils.py create mode 100644 configs/arraylink_1km.yaml create mode 100644 configs/upa_baseline.yaml create mode 100644 data/README.md create mode 100644 environment.yml create mode 100644 scripts/fig04_gain_vs_arrays.py create mode 100644 scripts/fig06_mimo_boundaries.py create mode 100644 scripts/fig09_beampattern_sim.py create mode 100644 scripts/fig10_hardware_validation.py create mode 100644 scripts/fig11_2d_beampattern.py create mode 100644 scripts/fig12_mimo_dof.py create mode 100644 scripts/fig13_throughput.py create mode 100644 tests/__init__.py create mode 100644 tests/test_array_geometry.py create mode 100644 tests/test_beamforming.py create mode 100644 tests/test_channel.py create mode 100644 tests/test_smoke.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..f3d668d --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,28 @@ +name: CI + +on: + push: + branches: ["**"] + pull_request: + branches: [main] + +jobs: + test: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v4 + + - name: Set up micromamba + uses: mamba-org/setup-micromamba@v2 + with: + environment-file: environment.yml + cache-environment: true + + - name: Run unit tests + run: pytest tests/test_channel.py tests/test_array_geometry.py tests/test_beamforming.py -v + shell: bash -el {0} + + - name: Run smoke tests + run: pytest tests/test_smoke.py -v --timeout=120 + shell: bash -el {0} diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..7c6c8d9 --- /dev/null +++ b/.gitignore @@ -0,0 +1,27 @@ +__pycache__/ +*.py[cod] +*.egg-info/ +dist/ +build/ +.eggs/ + +# Simulation outputs +paper_figures/ +results/ +*.h5 +*.dat + +# Temporary files +*.pkl +/tmp/ + +# IDE +.vscode/ +.idea/ + +# OS +.DS_Store + +# Pytest +.pytest_cache/ +.coverage diff --git a/Readme.md b/Readme.md index c942247..71b5a35 100644 --- a/Readme.md +++ b/Readme.md @@ -1,3 +1,120 @@ -# ArrayLink Repository +# ArrayLink Simulator -Code will be updated soon. +Near-field LoS MIMO simulator for the ArrayLink distributed ground station (INFOCOM 2026). + +ArrayLink uses 16 phased-array panels spread across a km-scale aperture to form a +distributed ground station that exploits spatial multiplexing in the radiative near-field +of LEO/MEO satellites. + +## Installation + +Requires [mamba](https://mamba.readthedocs.io/) or [micromamba](https://mamba.readthedocs.io/en/latest/installation/micromamba-installation.html). + +```bash +git clone https://github.com/ucsdwcsng/ArrayLink.git +cd ArrayLink +mamba env create -f environment.yml +mamba activate arraylink +``` + +For the interactive 2D beam pattern (Fig. 11 `--interactive` flag): +```bash +mamba install -c conda-forge plotly +``` + +## Quick Start + +```python +from arraylink.channel import compute_channel_matrix, mimo_region_bounds +from arraylink.array_geometry import build_ground_station +from arraylink.beamforming import total_array_gain_dbi, parabolic_gain_dbi + +# MIMO feasibility range for the ArrayLink aperture at 28 GHz +lam = 3e8 / 28e9 +r_min, r_max = mimo_region_bounds(d_tx=2000, d_rx=1.0, wavelength=lam) +print(f"MIMO region: {r_min/1e3:.0f} km – {r_max/1e3:.0f} km") + +# Build the ArrayLink ground station (16 panels, center-dense layout) +gnd = build_ground_station( + mode="arraylink", + subarray_shape=(32, 32), + element_spacing=lam / 2, + Nx=4, Ny=4, + Lx=1414.2, Ly=1000.0, +) +print(f"Total antenna elements: {len(gnd)}") +``` + +## Reproducing Paper Figures + +Run any figure script from the repo root: + +| Script | Figure | Description | +|--------|--------|-------------| +| `scripts/fig04_gain_vs_arrays.py` | Fig. 4 | Array gain vs number of panels | +| `scripts/fig06_mimo_boundaries.py` | Fig. 6 | Theoretical MIMO region boundaries | +| `scripts/fig09_beampattern_sim.py` | Fig. 9 | Simulation setup: UPA/ArrayLink positions and gain vs θ / distance | +| `scripts/fig10_hardware_validation.py` | Fig. 10 | Hardware experiment validation | +| `scripts/fig11_2d_beampattern.py` | Fig. 11 | 2D beam pattern heatmap | +| `scripts/fig12_mimo_dof.py` | Fig. 12 | MIMO degrees of freedom vs distance | +| `scripts/fig13_throughput.py` | Fig. 13 | Throughput comparison | + +```bash +# Generate all figures (outputs to paper_figures/) +for s in scripts/fig*.py; do python $s; done + +# Interactive 2D beam pattern +python scripts/fig11_2d_beampattern.py --interactive +``` + +### Hardware data (Fig. 10) + +Fig. 10 shows theory and simulation curves even without hardware data. To add the +measured curves, place `hardware_metrics.pkl` at `data/hardware_metrics.pkl`. See +[data/README.md](data/README.md) for the expected format. + +## Running Tests + +```bash +pytest tests/ -v +``` + +## Package Structure + +``` +arraylink/ + channel.py # Channel matrix, singular values, MIMO bounds, spectral efficiency + array_geometry.py # UPA, center-dense, and uniform array placement + beamforming.py # DC weights, beam pattern, array/dish gain formulas + utils.py # Coordinate transforms and helpers +scripts/ # One script per paper figure +configs/ # YAML parameter files for ArrayLink and UPA baseline +tests/ # Unit tests (42) + smoke tests (8) +data/ # Hardware experiment data (not included; see data/README.md) +environment.yml # Mamba/conda environment +``` + +## Key Assumptions + +| Parameter | Value | Notes | +|-----------|-------|-------| +| Carrier frequency | 28 GHz | Ka-band (λ ≈ 10.7 mm) | +| Hardware experiment frequency | 27 GHz | Matches lab setup | +| Panel element gain | 6 dBi | Isotropic radiator baseline | +| Dish efficiency | 0.6 (60%) | Standard aperture efficiency | +| MIMO threshold τ | 0.1 | σ₂/σ₁ > τ for spatial multiplexing | +| Panels | 16 (4×4 layout) | Center-dense over √2 km × 1 km | +| Elements per panel | 32×32 = 1024 | Half-wavelength spacing | +| Center-dense exponent γ | 3.0 | Power-law placement transform | + +## Citation + +If you use this simulator, please cite: + +```bibtex +@inproceedings{arraylink2026, + title = {ArrayLink: Distributed Ground Station for Near-Field LoS MIMO Satellite Communications}, + booktitle = {IEEE INFOCOM}, + year = {2026}, +} +``` diff --git a/arraylink/__init__.py b/arraylink/__init__.py new file mode 100644 index 0000000..38b6285 --- /dev/null +++ b/arraylink/__init__.py @@ -0,0 +1,33 @@ +"""ArrayLink near-field LoS MIMO simulator.""" +from .channel import ( + compute_channel_matrix, + compute_singular_values, + singular_value_ratio, + degrees_of_freedom, + spectral_efficiency, + mimo_region_bounds, + theoretical_singular_value_ratio, +) +from .array_geometry import ( + upa_positions, + center_dense_positions, + uniform_grid_positions, + place_subarrays, + build_ground_station, +) +from .beamforming import ( + parabolic_gain_dbi, + total_array_gain_dbi, + marginal_gain_dbi, + dc_weights, + compute_distances_batched, + compute_beam_pattern, + compute_beam_pattern_numpy, +) +from .utils import ( + spherical2cartesian, + cartesian2spherical, + linear2db, + db2linear, + generate_grid_points, +) diff --git a/arraylink/array_geometry.py b/arraylink/array_geometry.py new file mode 100644 index 0000000..50d91b5 --- /dev/null +++ b/arraylink/array_geometry.py @@ -0,0 +1,186 @@ +""" +Antenna array placement functions. + +Supports: + - Uniform Planar Arrays (UPA / monolithic) + - ArrayLink center-dense distributed placement + - Uniform grid distributed placement + - Subarray expansion (place an M×N panel at each base location) + +All positions returned as (N, 3) NumPy arrays in km (or the unit passed in). +""" +import numpy as np + + +# --------------------------------------------------------------------------- +# Monolithic UPA +# --------------------------------------------------------------------------- + +def upa_positions(M, N, element_spacing, center=None): + """ + Generate a M×N Uniform Planar Array centred at the origin (z=0 plane). + + Parameters + ---------- + M, N : element counts along x and y + element_spacing : element pitch (same unit as returned positions) + center : optional (3,) offset; defaults to [0,0,0] + + Returns + ------- + coords : (M*N, 3) antenna positions + """ + xs = (np.arange(M) - (M - 1) / 2.0) * element_spacing + ys = (np.arange(N) - (N - 1) / 2.0) * element_spacing + X, Y = np.meshgrid(xs, ys, indexing='ij') + coords = np.stack([X.ravel(), Y.ravel(), np.zeros(M * N)], axis=1) + if center is not None: + coords = coords + np.asarray(center) + return coords + + +# --------------------------------------------------------------------------- +# ArrayLink: center-dense panel base placement +# --------------------------------------------------------------------------- + +def center_dense_positions(Mx, My, Gx, Gy, gamma=3.0, seed=None, jitter_frac=0.0): + """ + Generate 2-D center-dense panel base positions using a power-law transform. + + Starting from a uniform Mx×My grid on [-Gx/2, Gx/2] × [-Gy/2, Gy/2], + apply s -> sign(s)*|s|^gamma to concentrate panels toward the centre, + then optionally add random jitter. + + Paper reference: Sec. III-D, Fig. 7b. + + Parameters + ---------- + Mx, My : number of panels along x and y + Gx, Gy : aperture dimensions (km) + gamma : power-law exponent; higher = more centre-concentrated (default 3.0) + seed : random seed for jitter reproducibility + jitter_frac : jitter as fraction of minimum inter-panel spacing (default 0 = no jitter) + + Returns + ------- + bases : (Mx*My, 3) panel base positions (z=0) + """ + tx = np.linspace(0, 1, Mx) + ty = np.linspace(0, 1, My) + sx = 2 * tx - 1 # -> [-1, 1] + sy = 2 * ty - 1 + sx = np.sign(sx) * np.abs(sx) ** gamma # power-law transform + sy = np.sign(sy) * np.abs(sy) ** gamma + x_pos = (Gx / 2.0) * sx + y_pos = (Gy / 2.0) * sy + + X, Y = np.meshgrid(x_pos, y_pos, indexing='ij') + coords = np.stack([X.ravel(), Y.ravel(), np.zeros(Mx * My)], axis=1) + + if jitter_frac > 0.0 and (Mx > 1 or My > 1): + rng = np.random.default_rng(seed) + min_dx = np.min(np.diff(x_pos)) if Mx > 1 else Gx + min_dy = np.min(np.diff(y_pos)) if My > 1 else Gy + jitter_scale = jitter_frac * min(min_dx, min_dy) + coords[:, :2] += rng.uniform(-jitter_scale, jitter_scale, (Mx * My, 2)) + + return coords + + +def uniform_grid_positions(Nx, Ny, Lx, Ly, center=None): + """ + Generate a uniform Nx×Ny grid spanning Lx×Ly, centred at origin. + + Parameters + ---------- + Nx, Ny : number of panels along x and y + Lx, Ly : physical extent (km) + center : optional (3,) offset + + Returns + ------- + bases : (Nx*Ny, 3) positions (z=0) + """ + xs = np.linspace(-Lx / 2.0, Lx / 2.0, Nx) + ys = np.linspace(-Ly / 2.0, Ly / 2.0, Ny) + X, Y = np.meshgrid(xs, ys, indexing='ij') + coords = np.stack([X.ravel(), Y.ravel(), np.zeros(Nx * Ny)], axis=1) + if center is not None: + coords = coords + np.asarray(center) + return coords + + +# --------------------------------------------------------------------------- +# Subarray placement +# --------------------------------------------------------------------------- + +def place_subarrays(base_positions, subarray_shape, element_spacing): + """ + Expand each base position into a subarray of shape (M, N). + + Parameters + ---------- + base_positions : (K, 3) panel base positions + subarray_shape : (M, N) elements per panel + element_spacing : element pitch within a panel + + Returns + ------- + all_positions : (K*M*N, 3) full antenna positions + """ + M, N = subarray_shape + offsets = [] + for i in range(M): + for j in range(N): + dx = (i - (M - 1) / 2.0) * element_spacing + dy = (j - (N - 1) / 2.0) * element_spacing + offsets.append([dx, dy, 0.0]) + offsets = np.array(offsets) # (M*N, 3) + + all_pos = [] + for base in base_positions: + all_pos.append(base + offsets) + return np.vstack(all_pos) + + +# --------------------------------------------------------------------------- +# Full ground-station factory +# --------------------------------------------------------------------------- + +def build_ground_station(mode, subarray_shape, element_spacing, + Nx, Ny, Lx, Ly, + gamma=3.0, seed=None, jitter_frac=0.0, + center=None): + """ + Build ground-station antenna positions. + + Parameters + ---------- + mode : 'arraylink' (center-dense) or 'uniform' + subarray_shape : (M, N) elements per panel + element_spacing: element pitch within a panel (km) + Nx, Ny : number of panels along x and y + Lx, Ly : aperture dimensions (km) + gamma : power-law exponent for 'arraylink' mode + seed : RNG seed for jitter + jitter_frac : jitter as fraction of min spacing + center : optional (3,) ground-station centre offset (km) + + Returns + ------- + ant_positions : (Nx*Ny*M*N, 3) all antenna positions + base_positions : (Nx*Ny, 3) panel base positions + """ + if mode == 'arraylink': + bases = center_dense_positions(Nx, Ny, Lx, Ly, gamma=gamma, + seed=seed, jitter_frac=jitter_frac) + elif mode == 'uniform': + bases = uniform_grid_positions(Nx, Ny, Lx, Ly) + else: + raise ValueError(f"Unknown mode '{mode}'. Choose 'arraylink' or 'uniform'.") + + if center is not None: + bases = bases + np.asarray(center) + + ants = place_subarrays(bases, subarray_shape, element_spacing) + return ants, bases diff --git a/arraylink/beamforming.py b/arraylink/beamforming.py new file mode 100644 index 0000000..525c60f --- /dev/null +++ b/arraylink/beamforming.py @@ -0,0 +1,307 @@ +""" +Beamforming weight generation and beam pattern computation. + +Heavy computation (satellite-scale) uses: + - numexpr for fast exp(-j*psi*d) + - numpy.memmap for out-of-core distance matrices + - h5py for compressed beam-pattern storage + +For small/quick runs (unit tests, smoke tests), pure-numpy fallbacks are used +when numexpr is not available. +""" +import os +import gc +import tempfile +import numpy as np + +try: + import numexpr as ne + _HAS_NUMEXPR = True +except ImportError: + _HAS_NUMEXPR = False + +try: + import h5py + _HAS_H5PY = True +except ImportError: + _HAS_H5PY = False + + +# --------------------------------------------------------------------------- +# Parabolic dish reference +# --------------------------------------------------------------------------- + +def parabolic_gain_dbi(D, frequency, efficiency=0.6): + """ + Directional gain of a parabolic dish antenna (dBi). + + G = eta * (pi*D/lambda)^2 + + Paper reference: Sec. II-A, Eq. before (2). + + Parameters + ---------- + D : dish diameter (metres) + frequency : carrier frequency (Hz) + efficiency: aperture efficiency eta (default 0.6, typical 0.5-0.7) + + Returns + ------- + gain_dbi : float + """ + c = 3e8 + lam = c / frequency + gain_linear = efficiency * (np.pi * D / lam) ** 2 + return 10.0 * np.log10(gain_linear) + + +# --------------------------------------------------------------------------- +# Gain aggregation model +# --------------------------------------------------------------------------- + +def total_array_gain_dbi(N, G_pa_dbi): + """ + Total gain of N coherently combined phased-array panels. + + G_total = 10*log10(N) + G_PA (paper Eq. (2), ignoring sync loss delta) + + Parameters + ---------- + N : number of panels (int or array) + G_pa_dbi: gain of a single panel (dBi) + + Returns + ------- + G_total_dbi : float or array + """ + return 10.0 * np.log10(np.asarray(N, dtype=float)) + G_pa_dbi + + +def marginal_gain_dbi(N_array, G_pa_dbi): + """ + Marginal gain from adding one more panel: G(N+1) - G(N). + + Parameters + ---------- + N_array : 1-D array of panel counts + G_pa_dbi : gain of a single panel (dBi) + + Returns + ------- + marginal : array of same length as N_array + """ + N = np.asarray(N_array, dtype=float) + G = total_array_gain_dbi(N + 1, G_pa_dbi) - total_array_gain_dbi(N, G_pa_dbi) + return G + + +# --------------------------------------------------------------------------- +# Delay-compensation (DC) weights +# --------------------------------------------------------------------------- + +def dc_weights(satellite_loc, rx_coords, wavelength): + """ + Compute delay-compensation beamforming weights. + + w_k = exp(-j * 2*pi/lambda * d_k) / ||w|| + + where d_k is the distance from the satellite to the k-th receive antenna. + + Convention: beam_pattern_numpy computes af @ conj(w) where + af_k = exp(-j*psi*d_k_eval). With w_k = exp(-j*psi*d_k_focal), + conj(w_k) = exp(+j*psi*d_k_focal), so the product is + exp(j*psi*(d_k_focal - d_k_eval)) which sums coherently (= N) at + boresight. Using +j here (phase conjugate of h) doubles the phase + instead of cancelling it and gives completely wrong results for + distributed arrays with aperture >> lambda. + + Parameters + ---------- + satellite_loc : (3,) or (1,3) satellite position (km) + rx_coords : (M, 3) receive antenna positions (km) + wavelength : carrier wavelength (km) + + Returns + ------- + w : (M, 1) complex normalised weight vector + """ + sat = np.asarray(satellite_loc, dtype=float).reshape(1, 3) + rx = np.asarray(rx_coords, dtype=float) + d = np.linalg.norm(rx - sat, axis=1) # (M,) + w = np.exp(-1j * 2 * np.pi / wavelength * d) + w = w / np.linalg.norm(w) + return w.reshape(-1, 1) + + +# --------------------------------------------------------------------------- +# Distance computation (batched, mmap-backed) +# --------------------------------------------------------------------------- + +def compute_distances_batched(tx_points, rx_antennas, batch_size=1000, + tmp_file=None): + """ + Compute Euclidean distances from every tx point to every rx antenna. + + Uses a memory-mapped file to avoid loading the full (Ntot×M) matrix + into RAM at once — essential for satellite-scale grids. + + Parameters + ---------- + tx_points : (Ntot, 3) transmit/satellite grid points + rx_antennas : (M, 3) receive antenna positions + batch_size : number of tx points processed per batch + tmp_file : optional path for the mmap file; uses tempfile if None + + Returns + ------- + distances : np.memmap of shape (Ntot, M), float64 + Caller is responsible for deleting tmp_file when done. + tmp_file : str, path to the mmap file + """ + tx = np.asarray(tx_points, dtype=np.float64) + rx = np.asarray(rx_antennas, dtype=np.float64) + Ntot, M = tx.shape[0], rx.shape[0] + + if tmp_file is None: + fd, tmp_file = tempfile.mkstemp(suffix='.dat', prefix='arraylink_dist_') + os.close(fd) + + distances = np.memmap(tmp_file, dtype=np.float64, mode='w+', shape=(Ntot, M)) + + for start in range(0, Ntot, batch_size): + end = min(start + batch_size, Ntot) + dx = tx[start:end, 0:1] - rx[:, 0] # (batch, M) + dy = tx[start:end, 1:2] - rx[:, 1] + dz = tx[start:end, 2:3] - rx[:, 2] + if _HAS_NUMEXPR: + distances[start:end] = ne.evaluate("sqrt(dx**2 + dy**2 + dz**2)") + else: + distances[start:end] = np.sqrt(dx**2 + dy**2 + dz**2) + + distances.flush() + return distances, tmp_file + + +# --------------------------------------------------------------------------- +# Beam pattern (h5py-backed, batched over the first axis of the grid) +# --------------------------------------------------------------------------- + +def compute_beam_pattern( + grid_shape, satellite_points, rx_coords, + weights, wavelength, + gain_variations=True, + outer_batch_size=10, + min_limit_db=-50, + output_file=None, + dist_batch_size=1000, +): + """ + Compute the beam pattern over a 3-D satellite grid. + + For each point (r, theta, phi) or (x, y, z) on the grid, evaluates + + P(point) = |A(point) * w*|^2 in dB + + where A is the array factor vector exp(-j*2*pi/lambda * d_k). + + When gain_variations=True the amplitude is scaled by the free-space + path loss (1/d), matching the near-field model. + + Results are stored in an HDF5 file to keep memory bounded. + + Parameters + ---------- + grid_shape : (D, T, P) tuple matching satellite_points first 3 dims + satellite_points : (D, T, P, 3) Cartesian satellite grid (km) + rx_coords : (M, 3) receive antenna positions (km) + weights : (M, 1) complex beamforming weights + wavelength : carrier wavelength (km) + gain_variations : include 1/d amplitude weighting (default True) + outer_batch_size : rows of dim-0 processed per iteration + min_limit_db : floor clamp for the dB output + output_file : path for HDF5 output; uses tempfile if None + dist_batch_size : inner batch size for distance computation + + Returns + ------- + output_file : str, path to the HDF5 file containing dataset 'beam_pattern_dB' + """ + if not _HAS_H5PY: + raise ImportError("h5py is required for compute_beam_pattern. " + "Install with: pip install h5py") + + D, T, P = grid_shape + psi = 2 * np.pi / wavelength + + if output_file is None: + fd, output_file = tempfile.mkstemp(suffix='.h5', prefix='arraylink_bp_') + os.close(fd) + + with h5py.File(output_file, 'w') as f: + dset = f.create_dataset( + 'beam_pattern_dB', shape=(D, T, P), + dtype=np.float64, chunks=True, compression='gzip' + ) + + for start in range(0, D, outer_batch_size): + end = min(start + outer_batch_size, D) + batch_pts = satellite_points[start:end] # (batch, T, P, 3) + flat_pts = batch_pts.reshape(-1, 3) # (batch*T*P, 3) + + dist_flat, tmp_dist = compute_distances_batched( + flat_pts, rx_coords, batch_size=dist_batch_size + ) + + # Array factor: exp(-j*psi*d) shape (batch*T*P, M) + if _HAS_NUMEXPR: + af = ne.evaluate("exp(-1j * psi * dist_flat)") + else: + af = np.exp(-1j * psi * dist_flat) + + if gain_variations: + # Scale by absolute satellite distance to account for path loss + abs_dist = np.linalg.norm(flat_pts, axis=1, keepdims=True) # (N,1) + abs_dist = np.where(abs_dist == 0, np.finfo(float).eps, abs_dist) + af = af / (dist_flat / abs_dist) # relative amplitude weighting + + # Beam pattern: (batch*T*P, 1) -> dB + bp_linear = af @ np.conj(weights) # (N, 1) + bp_db = 20.0 * np.log10(np.abs(bp_linear) + np.finfo(float).eps) + bp_db = np.clip(bp_db, min_limit_db, None) + dset[start:end] = bp_db.reshape(end - start, T, P) + + del dist_flat, af, bp_linear, bp_db + os.remove(tmp_dist) + gc.collect() + + return output_file + + +# --------------------------------------------------------------------------- +# Quick (pure-numpy, no h5py) beam pattern for small grids / tests +# --------------------------------------------------------------------------- + +def compute_beam_pattern_numpy(satellite_points, rx_coords, weights, wavelength, + min_limit_db=-50): + """ + Pure-numpy beam pattern computation for small grids (tests / quick mode). + + Parameters + ---------- + satellite_points : (Ntot, 3) Cartesian satellite grid (km) + rx_coords : (M, 3) receive antenna positions (km) + weights : (M, 1) complex weights + wavelength : carrier wavelength (km) + + Returns + ------- + bp_db : (Ntot,) beam pattern in dB + """ + d = np.linalg.norm( + satellite_points[:, None, :] - rx_coords[None, :, :], axis=-1 + ) # (Ntot, M) + psi = 2 * np.pi / wavelength + af = np.exp(-1j * psi * d) # (Ntot, M) + bp = af @ np.conj(weights) # (Ntot, 1) + bp_db = 20.0 * np.log10(np.abs(bp).ravel() + np.finfo(float).eps) + return np.clip(bp_db, min_limit_db, None) diff --git a/arraylink/channel.py b/arraylink/channel.py new file mode 100644 index 0000000..6dfc98b --- /dev/null +++ b/arraylink/channel.py @@ -0,0 +1,210 @@ +""" +Free-space LoS channel matrix and near-field MIMO analysis. + +All positions in the same unit (km recommended for satellite scale; +metre scale is fine for hardware experiments — just keep lambda consistent). + +Mathematical reference: ArrayLink paper (INFOCOM 2026), Sec. III-B/C and Appendix. +""" +import numpy as np + + +# --------------------------------------------------------------------------- +# Channel matrix +# --------------------------------------------------------------------------- + +def compute_channel_matrix(tx_pos, rx_pos, wavelength): + """ + Compute the Nr x Nt free-space spherical-wave channel matrix. + + h_ij = (lambda / (4*pi*d_ij)) * exp(-j*2*pi*d_ij / lambda) + + Assumption: amplitude ~1/d (near-field model); for far-field validation + only phase terms matter and amplitude differences are negligible. + + Parameters + ---------- + tx_pos : (Nt, 3) transmit antenna positions + rx_pos : (Nr, 3) receive antenna positions + wavelength : carrier wavelength (same units as positions) + + Returns + ------- + H : (Nr, Nt) complex channel matrix + """ + tx = np.asarray(tx_pos, dtype=float) # (Nt, 3) + rx = np.asarray(rx_pos, dtype=float) # (Nr, 3) + + # pairwise distances (Nr, Nt) + diff = rx[:, None, :] - tx[None, :, :] # (Nr, Nt, 3) + d = np.linalg.norm(diff, axis=-1) # (Nr, Nt) + + H = (wavelength / (4 * np.pi * d)) * np.exp(-1j * 2 * np.pi * d / wavelength) + return H + + +def compute_singular_values(H, normalise=True): + """ + SVD of H; optionally normalise so ||s||_2 = 1. + + Returns singular values in descending order. + """ + s = np.linalg.svd(H, compute_uv=False) # descending + if normalise and s[0] > 0: + s = s / np.linalg.norm(s) + return s + + +def singular_value_ratio(H): + """ + Return sigma_min / sigma_max from the normalised SVD of H. + For a 2×2 channel this equals sigma_2 / sigma_1. + Useful for checking MIMO feasibility (threshold tau ~ 0.1). + """ + s = compute_singular_values(H, normalise=True) + if s[0] == 0: + return 0.0 + return float(s[-1] / s[0]) + + +def degrees_of_freedom(H, threshold=0.1): + """ + Number of spatial streams whose normalised singular value exceeds `threshold`. + + threshold=0.1 matches the feasibility criterion in Eq. (7) of the paper. + """ + s = compute_singular_values(H, normalise=True) + return int(np.sum(s >= threshold)) + + +def spectral_efficiency(H, snr_linear): + """ + MIMO Shannon capacity under uniform power allocation (bits/s/Hz). + + C = sum_k log2(1 + SNR/Nt * sigma_k^2) + + Parameters + ---------- + H : (Nr, Nt) channel matrix + snr_linear : total receive SNR (linear) + """ + Nt = H.shape[1] + s = np.linalg.svd(H, compute_uv=False) + return float(np.sum(np.log2(1.0 + (snr_linear / Nt) * s**2))) + + +# --------------------------------------------------------------------------- +# MIMO region boundaries (2×2 closed-form, paper Sec. III-C) +# --------------------------------------------------------------------------- + +def mimo_region_bounds(d_tx, d_rx, wavelength, tau=0.1): + """ + Analytical MIMO feasibility boundaries for a 2×2 LoS system. + + Equations (9) and (11) from the paper. + + Parameters + ---------- + d_tx, d_rx : transmit and receive antenna spacings (any consistent unit) + wavelength : carrier wavelength (same unit) + tau : singular-value ratio threshold (default 0.1) + + Returns + ------- + r_min, r_max : distance boundaries (same unit as d_tx/d_rx/wavelength) + + Notes + ----- + At r < r_min the channel is unstable (oscillating singular-value ratio). + At r_min <= r <= r_max MIMO is well-conditioned. + At r > r_max the channel degrades toward rank-1. + """ + dtdrl = d_tx * d_rx / wavelength + r_min = (np.pi / (2 * np.arctan(1.0 / tau))) * dtdrl + r_max = (np.pi / (2 * np.arctan(tau))) * dtdrl + return float(r_min), float(r_max) + + +def theoretical_singular_value_ratio(r, d_tx, d_rx, wavelength): + """ + Closed-form sigma_2/sigma_1 vs distance for a 2×2 perpendicular LoS channel. + + From Eq. (4)-(5): Delta = 2*pi*d_tx*d_rx / (lambda * r) + sigma_1,2 = sqrt(2 +/- 2*|cos(Delta/2)|) + ratio = sigma_2 / sigma_1 + + Parameters + ---------- + r : distance (scalar or array, same unit as d_tx/d_rx/wavelength) + d_tx, d_rx : antenna spacings + wavelength : carrier wavelength + + Returns + ------- + ratio : sigma_2 / sigma_1 (in [0, 1]) + """ + r = np.asarray(r, dtype=float) + delta = 2 * np.pi * d_tx * d_rx / (wavelength * r) + s1 = np.sqrt(2 + 2 * np.abs(np.cos(delta / 2))) + s2 = np.sqrt(2 - 2 * np.abs(np.cos(delta / 2))) + norm = np.sqrt(s1**2 + s2**2) + s1n, s2n = s1 / norm, s2 / norm + return s2n / s1n + + +# --------------------------------------------------------------------------- +# Multi-stream DoF sweep (satellite scale) +# --------------------------------------------------------------------------- + +def dof_vs_distance(dist_array, gnd_coords, sat_params, wavelength, threshold=0.1): + """ + Compute singular values and DoF for a 4-element satellite array vs. distance. + + Parameters + ---------- + dist_array : 1-D array of distances (km) + gnd_coords : (M, 3) ground antenna positions (km) + sat_params : dict with keys: + center_theta, center_phi (radians) + subarray_shape [nx, ny] + xdim_subarrays, ydim_subarrays + Dx_km, Dy_km, seed + wavelength : carrier wavelength (km) + threshold : singular-value ratio threshold for DoF count + + Returns + ------- + dof_list : list of int, DoF at each distance + sing_vals_list : list of normalised singular-value arrays + """ + from .array_geometry import place_subarrays, center_dense_positions + + dof_list = [] + sing_vals_list = [] + + for d_km in dist_array: + # satellite centre in Cartesian + ct = sat_params['center_theta'] + cp = sat_params['center_phi'] + sat_center = np.array([ + d_km * np.sin(ct) * np.cos(cp), + d_km * np.sin(ct) * np.sin(cp), + d_km * np.cos(ct), + ]) + + # satellite antenna positions (small 2-D array centred at sat_center) + nx, ny = sat_params['subarray_shape'] + spacing = wavelength / 2.0 + sat_coords = place_subarrays( + np.array([sat_center]), + subarray_shape=(nx * sat_params['xdim_subarrays'], + ny * sat_params['ydim_subarrays']), + element_spacing=spacing, + ) + + H = compute_channel_matrix(sat_coords, gnd_coords, wavelength) + s = compute_singular_values(H, normalise=True) + dof_list.append(int(np.sum(s >= threshold))) + sing_vals_list.append(s) + + return dof_list, sing_vals_list diff --git a/arraylink/utils.py b/arraylink/utils.py new file mode 100644 index 0000000..d1fd6a7 --- /dev/null +++ b/arraylink/utils.py @@ -0,0 +1,77 @@ +""" +Coordinate conversion and decibel helpers. + +All coordinates in kilometres unless otherwise noted in the function signature. +""" +import numpy as np + + +def spherical2cartesian(r, theta, phi, deg=False): + """ + r : radius (km) + theta : inclination from +z axis (0 <= theta <= pi) + phi : azimuth from +x axis (-pi < phi <= pi) + """ + if deg: + theta = np.radians(theta) + phi = np.radians(phi) + x = r * np.sin(theta) * np.cos(phi) + y = r * np.sin(theta) * np.sin(phi) + z = r * np.cos(theta) + return x, y, z + + +def cartesian2spherical(x, y, z, deg=False): + """Returns (r, theta, phi). theta in [0, pi], phi in (-pi, pi].""" + r = np.sqrt(x**2 + y**2 + z**2) + r_safe = np.where(r == 0, np.finfo(float).eps, r) + theta = np.arccos(z / r_safe) + phi = np.arctan2(y, x) + if deg: + theta = np.degrees(theta) + phi = np.degrees(phi) + return r, theta, phi + + +def linear2db(x): + """Convert linear amplitude to dB (20*log10).""" + return 20.0 * np.log10(np.abs(x)) + + +def db2linear(x_db): + """Convert dB to linear amplitude.""" + return 10.0 ** (x_db / 20.0) + + +def generate_grid_points(x, y, z, spherical=False): + """ + Build a 3-D grid of Cartesian points. + + Parameters + ---------- + x, y, z : 1-D arrays + When spherical=False -> Cartesian axes (km). + When spherical=True -> r (km), theta (deg), phi (deg). + + Returns + ------- + grid : (Nx, Ny, Nz, 3) array of Cartesian coordinates + X, Y, Z : (Nx, Ny, Nz) component arrays + """ + x = np.asarray(x, dtype=float) + y = np.asarray(y, dtype=float) + z = np.asarray(z, dtype=float) + + if spherical: + r = x[:, None, None] + theta = np.deg2rad(y)[None, :, None] + phi = np.deg2rad(z)[None, None, :] + R, Theta, Phi = np.broadcast_arrays(r, theta, phi) + X = R * np.sin(Theta) * np.cos(Phi) + Y = R * np.sin(Theta) * np.sin(Phi) + Z = R * np.cos(Theta) + else: + X, Y, Z = np.meshgrid(x, y, z, indexing='ij') + + grid = np.stack((X, Y, Z), axis=-1) + return grid, X, Y, Z diff --git a/configs/arraylink_1km.yaml b/configs/arraylink_1km.yaml new file mode 100644 index 0000000..b7a8698 --- /dev/null +++ b/configs/arraylink_1km.yaml @@ -0,0 +1,56 @@ +# ArrayLink configuration: 16 x 32x32 panels, 1.414km x 1km aperture, 28 GHz +# Matches the satellite-scale simulation in the paper (Sec. IV-C, Fig. 9b) + +# System parameters +frequency_hz: 28.0e9 # 28 GHz (Ka-band) +# wavelength_km derived automatically as c / frequency + +# Ground station (ArrayLink) +ground_station: + mode: arraylink # 'arraylink' (center-dense) or 'uniform' + Nx: 4 # panels in x-direction + Ny: 4 # panels in y-direction + aperture_x_km: 1.414 # sqrt(2) km + aperture_y_km: 1.0 + subarray_shape: [32, 32] # elements per panel + gamma: 3.0 # power-law exponent for center-dense placement + jitter_frac: 0.0 # fractional random jitter (0 = deterministic) + seed: 0 + +# Satellite (transmitter) +satellite: + mode: uniform # uniform grid of tx antennas + Nx: 2 # tx elements in x + Ny: 2 # tx elements in y + aperture_x_m: 1.414 # 1.414 m (note: metres, not km) + aperture_y_m: 1.0 + +# Target satellite location for beamforming weights +target_satellite: + r_km: 500.0 # distance + theta_deg: 0.0 # elevation (0 = zenith) + phi_deg: 0.0 # azimuth + +# Beam pattern evaluation grid (spherical) +eval_grid: + r_km_min: 100.0 + r_km_max: 3000.0 + r_km_step: 50.0 + theta_deg_min: -10.0 + theta_deg_max: 10.0 + theta_deg_step: 0.1 + phi_deg_min: 0.0 + phi_deg_max: 1.0 + phi_deg_step: 1.0 + +# Beam pattern computation settings +computation: + gain_variations: true # include 1/d amplitude weighting + outer_batch_size: 20 # number of r slices per batch + dist_batch_size: 1000 + min_limit_db: -50 + element_gain_dbi: 6.0 # single-element gain added to beam pattern + +# Output +output: + paper_figures_dir: paper_figures diff --git a/configs/upa_baseline.yaml b/configs/upa_baseline.yaml new file mode 100644 index 0000000..1e1c7e2 --- /dev/null +++ b/configs/upa_baseline.yaml @@ -0,0 +1,48 @@ +# UPA baseline: 128x128 monolithic uniform planar array, 28 GHz +# Matches the monolithic comparison in the paper (Sec. IV-C, Fig. 9a) + +frequency_hz: 28.0e9 + +ground_station: + mode: uniform + Nx: 1 # single "panel" (the whole UPA) + Ny: 1 + aperture_x_km: 0.0 # not used in uniform mode when Nx=Ny=1 + aperture_y_km: 0.0 + subarray_shape: [128, 128] # full monolithic UPA + gamma: 3.0 + jitter_frac: 0.0 + seed: 0 + +satellite: + mode: uniform + Nx: 2 + Ny: 2 + aperture_x_m: 1.414 + aperture_y_m: 1.0 + +target_satellite: + r_km: 500.0 + theta_deg: 0.0 + phi_deg: 0.0 + +eval_grid: + r_km_min: 100.0 + r_km_max: 3000.0 + r_km_step: 50.0 + theta_deg_min: -10.0 + theta_deg_max: 10.0 + theta_deg_step: 0.1 + phi_deg_min: 0.0 + phi_deg_max: 1.0 + phi_deg_step: 1.0 + +computation: + gain_variations: true + outer_batch_size: 20 + dist_batch_size: 1000 + min_limit_db: -50 + element_gain_dbi: 6.0 + +output: + paper_figures_dir: paper_figures diff --git a/data/README.md b/data/README.md new file mode 100644 index 0000000..bc7bddb --- /dev/null +++ b/data/README.md @@ -0,0 +1,45 @@ +# Hardware Data + +The hardware experiment data for Fig. 10 will be released as part of the +ArrayLink dataset. + +## Expected file + +Place `hardware_metrics.pkl` in this directory (i.e. `data/hardware_metrics.pkl`). + +## Format + +The pickle file should contain a Python `dict` with the following structure: + +```python +{ + 'case1': { + 'distances': [2.5, 5.0, 10.0, ...], # metres + 'mean_sing_ratio': [0.82, 0.75, 0.61, ...], # mean sigma_2/sigma_1 + 'std_sing_ratio': [0.03, 0.04, 0.05, ...], # std dev across packets + }, + 'case2': { ... }, + 'case3': { ... }, + 'case4': { ... }, +} +``` + +where case1–case4 correspond to the four transmit/receive aperture configurations +in Table I of the paper: + +| Case | d_tx (cm) | d_rx (cm) | +|------|-----------|-----------| +| 1 | 20 | 20 | +| 2 | 50 | 20 | +| 3 | 20 | 50 | +| 4 | 50 | 50 | + +## Running Fig 10 without hardware data + +`fig10_hardware_validation.py` will run with theory and simulation curves only +if the pkl file is absent. It prints a clear message and continues: + +``` +Hardware data not found at: data/hardware_metrics.pkl +Proceeding with theory + simulation curves only. +``` diff --git a/environment.yml b/environment.yml new file mode 100644 index 0000000..19948b7 --- /dev/null +++ b/environment.yml @@ -0,0 +1,19 @@ +name: arraylink +channels: + - conda-forge + - defaults +dependencies: + - python=3.11 + - numpy>=1.24 + - scipy>=1.10 + - matplotlib>=3.7 + - numexpr>=2.8 + - h5py>=3.8 + - pyyaml>=6.0 + - pytest>=7.0 + - pip + - pip: + - scienceplots>=2.0 # not on conda-forge + +# Optional: interactive 2D beam pattern (fig11 --interactive flag) +# conda install -c conda-forge plotly diff --git a/scripts/fig04_gain_vs_arrays.py b/scripts/fig04_gain_vs_arrays.py new file mode 100644 index 0000000..44b7608 --- /dev/null +++ b/scripts/fig04_gain_vs_arrays.py @@ -0,0 +1,104 @@ +""" +Fig 4a/4b — Total and marginal gain vs number of phased-array panels. + +Reproduces Fig. 4 from the ArrayLink paper (INFOCOM 2026). + +Usage +----- + python scripts/fig04_gain_vs_arrays.py [--quick] [--save-dir paper_figures] +""" +import argparse +import os +import sys +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +import numpy as np +import matplotlib.pyplot as plt + +try: + import scienceplots + plt.style.use(['science', 'no-latex']) +except ImportError: + pass + +from arraylink.beamforming import total_array_gain_dbi, marginal_gain_dbi, parabolic_gain_dbi + + +def parse_args(): + ap = argparse.ArgumentParser(description="Fig 4: gain vs number of panels") + ap.add_argument("--quick", action="store_true", help="Skip plot, just print values") + ap.add_argument("--save-dir", default="paper_figures", help="Output directory for PDFs") + return ap.parse_args() + + +def main(): + args = parse_args() + os.makedirs(args.save_dir, exist_ok=True) + + # --- Parameters (paper Sec. III-A) --- + G_PA_DBI = 36.1 # single 32x32 panel gain (6 dBi/element × 1024 elements) + N_MAX = 50 + N = np.arange(1, N_MAX + 1) + + # Reference dish gains (paper Fig. 2) + G_DISH_185 = parabolic_gain_dbi(1.85, 28e9, efficiency=0.6) # ~52.6 dBi + G_DISH_147 = parabolic_gain_dbi(1.47, 28e9, efficiency=0.6) # ~49.5 dBi + ARRAYLINK_N = 16 # number of panels chosen in the paper + THRESHOLD_DB = 0.25 # marginal gain threshold beyond which returns diminish + + G_total = total_array_gain_dbi(N, G_PA_DBI) + G_marginal = marginal_gain_dbi(N, G_PA_DBI) + + if args.quick: + print(f"G(1) = {G_total[0]:.2f} dBi") + print(f"G(16) = {G_total[15]:.2f} dBi") + print(f"G(50) = {G_total[-1]:.2f} dBi") + print(f"Marginal gain at N=16: {G_marginal[15]:.4f} dB") + print(f"1.47m dish: {G_DISH_147:.2f} dBi, 1.85m dish: {G_DISH_185:.2f} dBi") + return + + fig_params = dict(fontsize=16, tick_labelsize=14) + + # --- Fig 4a: total gain vs N --- + fig, ax = plt.subplots(figsize=(6, 5)) + ax.plot(N, G_total, linewidth=2, color='steelblue', label='ArrayLink') + ax.axhline(G_DISH_185, linestyle='--', color='firebrick', linewidth=1.5, + label=f'1.85 m dish ({G_DISH_185:.1f} dBi)') + ax.axhline(G_DISH_147, linestyle=':', color='darkorange', linewidth=1.5, + label=f'1.47 m dish ({G_DISH_147:.1f} dBi)') + ax.axvline(ARRAYLINK_N, linestyle='--', color='gray', linewidth=1, + label=f'ArrayLink (N={ARRAYLINK_N})') + ax.set_xlabel("Number of phased-array panels", fontsize=fig_params['fontsize']) + ax.set_ylabel("Total gain (dBi)", fontsize=fig_params['fontsize']) + ax.tick_params(labelsize=fig_params['tick_labelsize']) + ax.legend(fontsize=13, framealpha=0.5) + ax.grid(True) + fig.tight_layout() + out = os.path.join(args.save_dir, "fig04a_total_gain_vs_arrays.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + + # --- Fig 4b: marginal gain vs N --- + fig, ax = plt.subplots(figsize=(6, 5)) + ax.plot(N, G_marginal, linewidth=2, color='steelblue', label='Marginal gain') + ax.axhline(THRESHOLD_DB, linestyle='--', color='firebrick', linewidth=1.5, + label=f'Threshold ({THRESHOLD_DB} dB)') + ax.axvline(ARRAYLINK_N, linestyle='--', color='gray', linewidth=1, + label=f'N={ARRAYLINK_N}') + ax.set_xlabel("Number of phased-array panels", fontsize=fig_params['fontsize']) + ax.set_ylabel("Marginal gain per added panel (dB)", fontsize=fig_params['fontsize']) + ax.tick_params(labelsize=fig_params['tick_labelsize']) + ax.legend(fontsize=13, framealpha=0.5) + ax.grid(True) + fig.tight_layout() + out = os.path.join(args.save_dir, "fig04b_marginal_gain_vs_arrays.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + + print(f"Saved fig04a and fig04b to {args.save_dir}/") + print(f" ArrayLink (N=16): {total_array_gain_dbi(16, G_PA_DBI):.2f} dBi " + f"vs 1.85m dish {G_DISH_185:.2f} dBi " + f"(gap = {G_DISH_185 - total_array_gain_dbi(16, G_PA_DBI):.1f} dB)") + + +if __name__ == "__main__": + main() diff --git a/scripts/fig06_mimo_boundaries.py b/scripts/fig06_mimo_boundaries.py new file mode 100644 index 0000000..396d565 --- /dev/null +++ b/scripts/fig06_mimo_boundaries.py @@ -0,0 +1,138 @@ +""" +Fig 6 — Singular-value ratio vs distance showing MIMO feasibility boundaries. + +Reproduces Fig. 6 from the ArrayLink paper (INFOCOM 2026). + +Plots sigma_2/sigma_1 vs link distance for small-aperture (hardware) and +satellite-scale aperture configurations, with r_min and r_max annotated. + +Usage +----- + python scripts/fig06_mimo_boundaries.py [--quick] [--save-dir paper_figures] +""" +import argparse +import os +import sys +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.patches as mpatches + +try: + import scienceplots + plt.style.use(['science', 'no-latex']) +except ImportError: + pass + +from arraylink.channel import mimo_region_bounds, theoretical_singular_value_ratio + + +def parse_args(): + ap = argparse.ArgumentParser(description="Fig 6: MIMO feasibility boundaries") + ap.add_argument("--quick", action="store_true") + ap.add_argument("--save-dir", default="paper_figures") + return ap.parse_args() + + +def _plot_boundary(ax, r_km, ratio, r_min, r_max, tau, label_prefix, color, lam): + ax.plot(r_km, ratio, linewidth=2, color=color, label=label_prefix) + ax.axhline(tau, linestyle='--', color='firebrick', linewidth=1.5, label=f'Threshold (τ={tau})') + ax.axvline(r_min, linestyle=':', color='green', linewidth=1.5, + label=rf'$r_{{min}}$ = {r_min:.2f}') + ax.axvline(r_max, linestyle=':', color='navy', linewidth=1.5, + label=rf'$r_{{max}}$ = {r_max:.2f}') + # shade MIMO-feasible region + ylim = ax.get_ylim() + ax.axvspan(r_min, r_max, alpha=0.12, color='green', label='MIMO feasible region') + + +def main(): + args = parse_args() + os.makedirs(args.save_dir, exist_ok=True) + + TAU = 0.1 # feasibility threshold (paper Eq. 7) + F_HZ = 28e9 # 28 GHz + C = 3e8 + LAM_M = C / F_HZ # wavelength in metres + LAM_KM = LAM_M / 1e3 # wavelength in km + + fig_params = dict(fontsize=16, tick_labelsize=14) + + # ---- (a) Small aperture: hardware scale (d_tx = d_rx = 0.2 m) ---- + if args.quick: + for d_tx, d_rx, label in [ + (0.2, 0.2, "d_tx=d_rx=0.2m"), # 20 cm x 20 cm + (np.sqrt(2e3), np.sqrt(2), "satellite scale"), # sqrt(2) km x sqrt(2) m + ]: + r_min, r_max = mimo_region_bounds(d_tx, d_rx, LAM_M, tau=TAU) + print(f"{label}: r_min={r_min:.2f}, r_max={r_max:.2f}") + return + + configs = [ + # (d_tx_m, d_rx_m, r_range_m, label, filename_suffix) + (0.20, 0.20, np.linspace(0.1, 100, 5000), + r"$d_{tx}=d_{rx}=20\,\mathrm{cm}$", "small"), + (0.20, 0.50, np.linspace(0.1, 200, 5000), + r"$d_{tx}=20\,\mathrm{cm},\;d_{rx}=50\,\mathrm{cm}$", "medium"), + ] + + for d_tx, d_rx, r_array, label, suffix in configs: + ratio = theoretical_singular_value_ratio(r_array, d_tx, d_rx, LAM_M) + r_min, r_max = mimo_region_bounds(d_tx, d_rx, LAM_M, tau=TAU) + + fig, ax = plt.subplots(figsize=(6, 5)) + ax.plot(r_array, ratio, linewidth=2, color='steelblue', label='Theory') + ax.axhline(TAU, linestyle='--', color='firebrick', linewidth=1.5, + label=f'Threshold (τ={TAU})') + ax.axvspan(r_min, r_max, alpha=0.15, color='green', label='MIMO feasible region') + ax.axvline(r_min, linestyle=':', color='green', linewidth=1.5, + label=rf'$r_{{min}}={r_min:.1f}\,\mathrm{{m}}$') + ax.axvline(r_max, linestyle=':', color='navy', linewidth=1.5, + label=rf'$r_{{max}}={r_max:.1f}\,\mathrm{{m}}$') + ax.set_xlabel("Distance (m)", fontsize=fig_params['fontsize']) + ax.set_ylabel(r"Singular-value ratio $\sigma_2/\sigma_1$", + fontsize=fig_params['fontsize']) + ax.set_title(label, fontsize=14) + ax.set_ylim(0, 1.05) + ax.tick_params(labelsize=fig_params['tick_labelsize']) + ax.legend(fontsize=12, framealpha=0.5, loc='upper right') + ax.grid(True) + fig.tight_layout() + out = os.path.join(args.save_dir, f"fig06_mimo_boundaries_{suffix}.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out} r_min={r_min:.2f} m, r_max={r_max:.2f} m") + + # ---- Satellite scale (km units) ---- + D_TX_KM = np.sqrt(2) # sqrt(2) km ground aperture diagonal + D_RX_KM = np.sqrt(2) / 1e3 # sqrt(2) m satellite aperture in km + r_km = np.linspace(1, 3200, 10000) + ratio_sat = theoretical_singular_value_ratio(r_km, D_TX_KM, D_RX_KM, LAM_KM) + r_min_km, r_max_km = mimo_region_bounds(D_TX_KM, D_RX_KM, LAM_KM, tau=TAU) + + fig, ax = plt.subplots(figsize=(6, 5)) + ax.plot(r_km, ratio_sat, linewidth=2, color='steelblue', label='Theory') + ax.axhline(TAU, linestyle='--', color='firebrick', linewidth=1.5, + label=f'Threshold (τ={TAU})') + ax.axvspan(r_min_km, r_max_km, alpha=0.15, color='green', label='MIMO feasible region') + ax.axvline(r_min_km, linestyle=':', color='green', linewidth=1.5, + label=rf'$r_{{min}}={r_min_km:.0f}\,\mathrm{{km}}$') + ax.axvline(r_max_km, linestyle=':', color='navy', linewidth=1.5, + label=rf'$r_{{max}}={r_max_km:.0f}\,\mathrm{{km}}$') + ax.set_xlabel("Distance (km)", fontsize=fig_params['fontsize']) + ax.set_ylabel(r"Singular-value ratio $\sigma_2/\sigma_1$", + fontsize=fig_params['fontsize']) + ax.set_ylim(0, 1.05) + ax.set_xlim(0, 3200) + ax.tick_params(labelsize=fig_params['tick_labelsize']) + ax.legend(fontsize=12, framealpha=0.5, loc='upper right') + ax.grid(True) + fig.tight_layout() + out = os.path.join(args.save_dir, "fig06_mimo_boundaries_satellite.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out} r_min={r_min_km:.0f} km, r_max={r_max_km:.0f} km") + + +if __name__ == "__main__": + main() diff --git a/scripts/fig09_beampattern_sim.py b/scripts/fig09_beampattern_sim.py new file mode 100644 index 0000000..6617c8c --- /dev/null +++ b/scripts/fig09_beampattern_sim.py @@ -0,0 +1,339 @@ +""" +Fig 9 — Simulation setup and beam pattern results. + +(a) 2D UPA 128×128 element positions (λ/2 spacing) +(b) ArrayLink: 16 × 32×32 panel center positions in a 2 km × 2 km grid +(c) Beam pattern vs elevation angle θ ∈ [−90°, +90°] at r = 500 km +(d) Beam pattern vs distance r ∈ [0, 2000] km at θ = 0° (boresight) + +Two ArrayLink layouts are shown (seed_0, seed_1 = two different random placements). +--center-dense adds a third ArrayLink curve using a 4×4 power-law grid. + +Panel placement (random layout): + - 4 corner panels fixed at (±1 km, ±1 km) to span the full aperture + - Remaining 12 panels drawn uniformly within the 2 km × 2 km grid + - Minimum separation between any two panel centres enforced + +Usage +----- + python scripts/fig09_beampattern_sim.py [--quick] [--center-dense] + [--save-dir paper_figures] +""" +import argparse +import os +import sys +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.patches import Rectangle + +try: + import scienceplots + plt.style.use(['science', 'no-latex']) +except ImportError: + pass + +from arraylink.array_geometry import upa_positions, place_subarrays +from arraylink.beamforming import dc_weights, compute_beam_pattern_numpy + +# -------------------------------------------------------------------------- +# Constants +# -------------------------------------------------------------------------- +F_HZ = 28e9 +LAM_KM = 3e8 / F_HZ / 1e3 # ~1.0714e-5 km +R_FOCAL = 500.0 # km +GRID_KM = 2.0 # panel centres within ±1 km × ±1 km +MIN_GAP = 0.1 # minimum separation between panel centres [km] +N_PANELS = 16 +SUB_SHAPE = (32, 32) # elements per panel +N_UPA = 128 # UPA dimension (128×128) +ELEM_GAIN = 6.0 # element gain assumption [dBi] + + +# -------------------------------------------------------------------------- +# Panel placement helpers +# -------------------------------------------------------------------------- +def _place_random_2d(n_panels, grid_km, min_gap_km, seed): + """Random panel centres; 4 corners fixed, remainder uniformly drawn.""" + rng = np.random.default_rng(seed) + half = grid_km / 2.0 + panels = [(-half, -half), (-half, half), (half, -half), (half, half)] + attempts = 0 + while len(panels) < n_panels and attempts < 500_000: + p = rng.uniform(-half, half, size=2) + if all(np.hypot(p[0] - q[0], p[1] - q[1]) >= min_gap_km for q in panels): + panels.append(tuple(p)) + attempts += 1 + if len(panels) < n_panels: + raise RuntimeError( + f"Could only place {len(panels)}/{n_panels} panels " + f"with min_gap={min_gap_km} km — reduce MIN_GAP or increase GRID_KM." + ) + return np.array(panels) # (n_panels, 2) + + +def _place_center_dense_2d(nx, ny, grid_km, gamma=3.0): + """4×4 power-law grid — denser at centre (ArrayLink placement).""" + half = grid_km / 2.0 + def warp(n): + t = np.linspace(0, 1, n) + s = 2 * t - 1 + return (half * np.sign(s) * np.abs(s) ** gamma) + xs = warp(nx) + ys = warp(ny) + XX, YY = np.meshgrid(xs, ys) + return np.column_stack([XX.ravel(), YY.ravel()]) # (nx*ny, 2) + + +def _build_arraylink(panel_xy_km, sub_shape, elem_spacing_km): + """All element positions for a distributed array with given panel centres.""" + centres_3d = np.hstack([panel_xy_km, np.zeros((len(panel_xy_km), 1))]) + return place_subarrays(centres_3d, sub_shape, elem_spacing_km) # (N, 3) km + + +# -------------------------------------------------------------------------- +# Beam pattern helpers +# -------------------------------------------------------------------------- +def _gain_vs_theta(ants_km, w, r_km, theta_deg): + thetas = np.deg2rad(theta_deg) + sat_pts = np.column_stack([ + r_km * np.sin(thetas), + np.zeros_like(thetas), + r_km * np.cos(thetas), + ]) + return compute_beam_pattern_numpy(sat_pts, ants_km, w, LAM_KM) # dB + + +def _gain_vs_dist(ants_km, w, r_km_array): + sat_pts = np.column_stack([ + np.zeros(len(r_km_array)), + np.zeros(len(r_km_array)), + r_km_array, + ]) + return compute_beam_pattern_numpy(sat_pts, ants_km, w, LAM_KM) # dB + + +# -------------------------------------------------------------------------- +# Main +# -------------------------------------------------------------------------- +def parse_args(): + ap = argparse.ArgumentParser(description="Fig 9: simulation setup and beam pattern") + ap.add_argument("--quick", action="store_true", + help="Coarse grids for fast smoke test") + ap.add_argument("--center-dense", action="store_true", + help="Add a third ArrayLink curve with center-dense 4×4 layout") + ap.add_argument("--save-dir", default="paper_figures") + return ap.parse_args() + + +def main(): + args = parse_args() + os.makedirs(args.save_dir, exist_ok=True) + + if args.quick: + theta_deg = np.linspace(-90, 90, 19) + r_axis_km = np.linspace(50, 2000, 11) + else: + theta_deg = np.linspace(-90, 90, 1801) + r_axis_km = np.linspace(10, 2000, 500) + + # ------------------------------------------------------------------ # + # Build arrays + # ------------------------------------------------------------------ # + # UPA baseline + upa_ants = upa_positions(N_UPA, N_UPA, LAM_KM / 2) # (16384, 3) km + + # ArrayLink: two random seeds + al_xy_s0 = _place_random_2d(N_PANELS, GRID_KM, MIN_GAP, seed=0) + al_ants_s0 = _build_arraylink(al_xy_s0, SUB_SHAPE, LAM_KM / 2) + + al_xy_s1 = _place_random_2d(N_PANELS, GRID_KM, MIN_GAP, seed=1) + al_ants_s1 = _build_arraylink(al_xy_s1, SUB_SHAPE, LAM_KM / 2) + + # Optional center-dense + if args.center_dense: + al_xy_cd = _place_center_dense_2d(4, 4, GRID_KM, gamma=3.0) + al_ants_cd = _build_arraylink(al_xy_cd, SUB_SHAPE, LAM_KM / 2) + + sat_boresight = np.array([0.0, 0.0, R_FOCAL]) + + # DC weights (focused at 500 km boresight) + w_upa = dc_weights(sat_boresight, upa_ants, LAM_KM) + w_s0 = dc_weights(sat_boresight, al_ants_s0, LAM_KM) + w_s1 = dc_weights(sat_boresight, al_ants_s1, LAM_KM) + if args.center_dense: + w_cd = dc_weights(sat_boresight, al_ants_cd, LAM_KM) + + print(f"UPA elements : {len(upa_ants)}") + print(f"ArrayLink elems: {len(al_ants_s0)}") + + # ------------------------------------------------------------------ # + # Fig 9a — UPA element positions + # ------------------------------------------------------------------ # + fig, ax = plt.subplots(figsize=(4, 4)) + ax.scatter(upa_ants[:, 0] * 1e3, upa_ants[:, 1] * 1e3, + s=0.3, color='royalblue', alpha=0.6) + ax.set_xlabel("X [m]", fontsize=11) + ax.set_ylabel("Y [m]", fontsize=11) + ax.set_title(f"UPA {N_UPA}×{N_UPA} (λ/2 spacing)", fontsize=11) + ax.set_aspect('equal') + ax.grid(True, alpha=0.3) + + # Inset: zoom on top-right 32×32 subarray block + corner_idx = [r * N_UPA + c + for r in range(N_UPA - 32, N_UPA) + for c in range(N_UPA - 32, N_UPA)] + corner_elems = upa_ants[corner_idx] + x_lo = corner_elems[:, 0].min() * 1e3 + x_hi = corner_elems[:, 0].max() * 1e3 + y_lo = corner_elems[:, 1].min() * 1e3 + y_hi = corner_elems[:, 1].max() * 1e3 + pad = max(x_hi - x_lo, y_hi - y_lo) * 0.2 + axins = ax.inset_axes([0.02, 0.58, 0.40, 0.40]) + axins.scatter(corner_elems[:, 0] * 1e3, corner_elems[:, 1] * 1e3, + s=3, color='royalblue', alpha=0.8) + axins.set_xlim(x_lo - pad, x_hi + pad) + axins.set_ylim(y_lo - pad, y_hi + pad) + axins.set_xticklabels([]) + axins.set_yticklabels([]) + axins.grid(True, alpha=0.3) + ax.add_patch(Rectangle( + (x_lo - pad, y_lo - pad), (x_hi - x_lo + 2 * pad), (y_hi - y_lo + 2 * pad), + linewidth=1.5, edgecolor='darkorange', facecolor='none' + )) + + fig.tight_layout() + out = os.path.join(args.save_dir, "fig09a_upa_positions.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out}") + + # ------------------------------------------------------------------ # + # Fig 9b — ArrayLink panel positions (seed 0) + # ------------------------------------------------------------------ # + fig, ax = plt.subplots(figsize=(4, 4)) + # All elements (light) + ax.scatter(al_ants_s0[:, 0], al_ants_s0[:, 1], + s=0.3, color='royalblue', alpha=0.3, label="Subarray Elements") + # Panel centres (prominent) + ax.scatter(al_xy_s0[:, 0], al_xy_s0[:, 1], + s=60, color='crimson', marker='x', linewidths=1.5, + label="Subarray Centers") + ax.set_xlabel("X [km]", fontsize=11) + ax.set_ylabel("Y [km]", fontsize=11) + ax.set_title(f"ArrayLink: {N_PANELS} × {SUB_SHAPE[0]}×{SUB_SHAPE[1]} UPAs", fontsize=11) + ax.set_aspect('equal') + ax.grid(True, alpha=0.3) + + # Inset: zoom on one panel (closest to grid centre) + dists_from_origin = np.linalg.norm(al_xy_s0, axis=1) + panel_idx = np.argmin(dists_from_origin) + cx, cy = al_xy_s0[panel_idx] + start = panel_idx * SUB_SHAPE[0] * SUB_SHAPE[1] + end = start + SUB_SHAPE[0] * SUB_SHAPE[1] + sub_e = al_ants_s0[start:end] + x_lo = sub_e[:, 0].min() + x_hi = sub_e[:, 0].max() + y_lo = sub_e[:, 1].min() + y_hi = sub_e[:, 1].max() + pad = max(x_hi - x_lo, y_hi - y_lo) * 0.25 + axins = ax.inset_axes([0.02, 0.58, 0.40, 0.40]) + axins.scatter(sub_e[:, 0], sub_e[:, 1], s=2, color='royalblue', alpha=0.8) + axins.scatter([cx], [cy], s=40, color='crimson', marker='x', linewidths=1.5) + axins.set_xlim(x_lo - pad, x_hi + pad) + axins.set_ylim(y_lo - pad, y_hi + pad) + axins.set_xticklabels([]) + axins.set_yticklabels([]) + axins.grid(True, alpha=0.3) + + legend_elems = [ + plt.Line2D([0], [0], marker='o', color='w', markerfacecolor='royalblue', + markersize=5, label='Subarray Elements'), + plt.Line2D([0], [0], marker='x', color='crimson', markersize=7, + linewidth=0, label='Subarray Centers'), + ] + ax.legend(handles=legend_elems, fontsize=8, loc='lower right') + fig.tight_layout() + out = os.path.join(args.save_dir, "fig09b_arraylink_positions.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out}") + + if args.quick: + print("Quick mode: running beam pattern on coarse grid.") + + # ------------------------------------------------------------------ # + # Compute beam patterns + # ------------------------------------------------------------------ # + bp_upa_theta = _gain_vs_theta(upa_ants, w_upa, R_FOCAL, theta_deg) + ELEM_GAIN + bp_s0_theta = _gain_vs_theta(al_ants_s0, w_s0, R_FOCAL, theta_deg) + ELEM_GAIN + bp_s1_theta = _gain_vs_theta(al_ants_s1, w_s1, R_FOCAL, theta_deg) + ELEM_GAIN + + bp_upa_dist = _gain_vs_dist(upa_ants, w_upa, r_axis_km) + ELEM_GAIN + bp_s0_dist = _gain_vs_dist(al_ants_s0, w_s0, r_axis_km) + ELEM_GAIN + bp_s1_dist = _gain_vs_dist(al_ants_s1, w_s1, r_axis_km) + ELEM_GAIN + + if args.center_dense: + bp_cd_theta = _gain_vs_theta(al_ants_cd, w_cd, R_FOCAL, theta_deg) + ELEM_GAIN + bp_cd_dist = _gain_vs_dist(al_ants_cd, w_cd, r_axis_km) + ELEM_GAIN + + boresight_upa = bp_upa_dist[np.argmin(np.abs(r_axis_km - R_FOCAL))] + boresight_s0 = bp_s0_dist[np.argmin(np.abs(r_axis_km - R_FOCAL))] + print(f"Boresight gain UPA : {boresight_upa:.1f} dBi") + print(f"Boresight gain ArrayLink: {boresight_s0:.1f} dBi") + + fig_kw = dict(fontsize=14, tick_labelsize=12) + + # ------------------------------------------------------------------ # + # Fig 9c — Gain vs angle + # ------------------------------------------------------------------ # + fig, ax = plt.subplots(figsize=(6, 4.5)) + ax.plot(theta_deg, bp_upa_theta, color='steelblue', linewidth=1.5, label="UPA") + ax.plot(theta_deg, bp_s0_theta, color='darkorange', linewidth=1.2, + linestyle='--', label="ArrayLink") + ax.plot(theta_deg, bp_s1_theta, color='seagreen', linewidth=1.2, + linestyle='-.', label="ArrayLink") + if args.center_dense: + ax.plot(theta_deg, bp_cd_theta, color='purple', linewidth=1.2, + linestyle=':', label="ArrayLink (center-dense)") + ax.set_xlabel(r"$\theta$ (in deg)", fontsize=fig_kw['fontsize']) + ax.set_ylabel("Gain (dB)", fontsize=fig_kw['fontsize']) + ax.set_xlim(-90, 90) + ax.tick_params(labelsize=fig_kw['tick_labelsize']) + ax.legend(fontsize=11, framealpha=0.5) + ax.grid(True, alpha=0.3) + fig.tight_layout() + out = os.path.join(args.save_dir, "fig09c_gain_vs_angle.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out}") + + # ------------------------------------------------------------------ # + # Fig 9d — Gain vs distance + # ------------------------------------------------------------------ # + fig, ax = plt.subplots(figsize=(6, 4.5)) + ax.plot(r_axis_km, bp_upa_dist, color='steelblue', linewidth=1.5, label="UPA") + ax.plot(r_axis_km, bp_s0_dist, color='darkorange', linewidth=1.2, + linestyle='--', label="ArrayLink") + ax.plot(r_axis_km, bp_s1_dist, color='seagreen', linewidth=1.2, + linestyle='-.', label="ArrayLink") + if args.center_dense: + ax.plot(r_axis_km, bp_cd_dist, color='purple', linewidth=1.2, + linestyle=':', label="ArrayLink (center-dense)") + ax.set_xlabel("Distance (km)", fontsize=fig_kw['fontsize']) + ax.set_ylabel("Gain (dB)", fontsize=fig_kw['fontsize']) + ax.set_xlim(0, 2000) + ax.axvline(R_FOCAL, color='gray', linestyle=':', linewidth=1, + label=f"Focal point ({R_FOCAL:.0f} km)") + ax.tick_params(labelsize=fig_kw['tick_labelsize']) + ax.legend(fontsize=11, framealpha=0.5) + ax.grid(True, alpha=0.3) + fig.tight_layout() + out = os.path.join(args.save_dir, "fig09d_gain_vs_distance.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out}") + + +if __name__ == "__main__": + main() diff --git a/scripts/fig10_hardware_validation.py b/scripts/fig10_hardware_validation.py new file mode 100644 index 0000000..87e228f --- /dev/null +++ b/scripts/fig10_hardware_validation.py @@ -0,0 +1,182 @@ +""" +Fig 10 — Hardware validation: singular-value ratio vs distance (4 aperture cases). + +Reproduces Fig. 10 from the ArrayLink paper (INFOCOM 2026). + +Overlays three curves for each aperture case (d_tx, d_rx): + - Hardware measurements (from hardware_metrics.pkl) + - Simulation (computed here from the channel model, Eq. 3) + - Theory (closed-form, Eqs. 4-5) + +Hardware data +------------- +Place `hardware_metrics.pkl` in the `data/` directory (default path). +The pkl file must be a dict keyed by case name with entries: + { + 'case1': {'distances': [...], 'mean_sing_ratio': [...], 'std_sing_ratio': [...]}, + ... + } + +If the file is not found, the script prints instructions and exits. + +Usage +----- + python scripts/fig10_hardware_validation.py [--quick] + [--data-path data/hardware_metrics.pkl] + [--save-dir paper_figures] +""" +import argparse +import os +import sys +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +import pickle +import numpy as np +import matplotlib.pyplot as plt + +try: + import scienceplots + plt.style.use(['science', 'no-latex']) +except ImportError: + pass + +from arraylink.channel import ( + compute_channel_matrix, compute_singular_values, + theoretical_singular_value_ratio, +) +from arraylink.array_geometry import upa_positions + + +# Hardware experiment aperture configurations (paper Sec. IV-B, Table I) +CASES = { + 'case1': {'d_tx_cm': 20, 'd_rx_cm': 20}, + 'case2': {'d_tx_cm': 50, 'd_rx_cm': 20}, + 'case3': {'d_tx_cm': 20, 'd_rx_cm': 50}, + 'case4': {'d_tx_cm': 50, 'd_rx_cm': 50}, +} + + +def parse_args(): + ap = argparse.ArgumentParser(description="Fig 10: hardware validation") + ap.add_argument("--quick", action="store_true", + help="Skip hardware file check; plot theory + sim only") + ap.add_argument("--data-path", default="data/hardware_metrics.pkl", + help="Path to hardware_metrics.pkl (default: data/hardware_metrics.pkl)") + ap.add_argument("--save-dir", default="paper_figures") + return ap.parse_args() + + +def simulate_ratio(d_tx_m, d_rx_m, wavelength_m, dist_m_array): + """Compute sigma_2/sigma_1 from the exact channel model (Eq. 3).""" + ratios = [] + for r_m in dist_m_array: + # 2-element TX and RX arrays, perpendicular to LoS + tx_pos = np.array([[-d_tx_m / 2, 0, 0], [d_tx_m / 2, 0, 0]]) + rx_pos = np.array([[0, -d_rx_m / 2, r_m], [0, d_rx_m / 2, r_m]]) + H = compute_channel_matrix(tx_pos, rx_pos, wavelength_m) + s = compute_singular_values(H, normalise=True) + ratios.append(float(s[1] / s[0]) if s[0] > 0 else 0.0) + return np.array(ratios) + + +def main(): + args = parse_args() + os.makedirs(args.save_dir, exist_ok=True) + + F_HZ = 27e9 # hardware experiment carrier (27 GHz) + LAM_M = 3e8 / F_HZ # ~0.0111 m + TAU = 0.1 + + # Try to load hardware data + hw_data = None + data_path = args.data_path + if not os.path.isabs(data_path): + # resolve relative to repo root (parent of scripts/) + script_dir = os.path.dirname(os.path.abspath(__file__)) + data_path = os.path.join(os.path.dirname(script_dir), data_path) + + if os.path.exists(data_path): + with open(data_path, 'rb') as f: + hw_data = pickle.load(f) + print(f"Loaded hardware data from {data_path}") + else: + if args.quick: + print(f"[quick mode] Hardware data not found at {data_path}. " + "Plotting theory + simulation only.") + else: + print(f"\nHardware data not found at: {data_path}") + print("To include hardware measurements in Fig 10:") + print(" 1. Obtain hardware_metrics.pkl from the data release") + print(" 2. Place it in data/hardware_metrics.pkl (or pass --data-path)") + print("\nProceeding with theory + simulation curves only.\n") + + # Distance range for simulation / theory + if args.quick: + dist_theory_m = np.linspace(1, 100, 50) + dist_sim_m = np.array([2.5, 5.0, 10.0, 20.0, 50.0, 100.0]) + else: + dist_theory_m = np.linspace(0.5, 100, 2000) + dist_sim_m = np.array([2.5, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0, + 50.0, 60.0, 70.0, 80.0, 90.0, 100.0]) + + fig_params = dict(fontsize=16, tick_labelsize=14) + + for case_name, cfg in CASES.items(): + d_tx = cfg['d_tx_cm'] / 100.0 # metres + d_rx = cfg['d_rx_cm'] / 100.0 + + theory_ratio = theoretical_singular_value_ratio(dist_theory_m, d_tx, d_rx, LAM_M) + sim_ratio = simulate_ratio(d_tx, d_rx, LAM_M, dist_sim_m) + + fig, ax = plt.subplots(figsize=(6, 4)) + + # Theory + ax.plot(dist_theory_m, theory_ratio, linestyle=':', linewidth=1.5, + color='black', label='Theory (Eq. 4-5)') + # Simulation + ax.plot(dist_sim_m, sim_ratio, linestyle='--', linewidth=2, + color='steelblue', label='Simulation (Eq. 3)') + # Hardware (if available) + if hw_data is not None and case_name in hw_data: + hw = hw_data[case_name] + ax.errorbar( + hw['distances'], hw['mean_sing_ratio'], + yerr=hw['std_sing_ratio'], + fmt='-o', linewidth=2, markersize=5, capsize=3, + color='darkorange', label='Hardware', + ) + + # MIMO feasible region + from arraylink.channel import mimo_region_bounds + r_min, r_max = mimo_region_bounds(d_tx, d_rx, LAM_M, tau=TAU) + ax.axhline(TAU, linestyle='--', color='firebrick', linewidth=1.5, + label=f'Threshold (τ={TAU})') + r_min_vis = max(r_min, dist_theory_m[0]) + r_max_vis = min(r_max, dist_theory_m[-1]) + if r_max_vis > r_min_vis: + ax.axvspan(r_min_vis, r_max_vis, alpha=0.12, color='green', + label='MIMO feasible') + + ax.set_title( + rf"$d_{{tx}}={cfg['d_tx_cm']}\,\mathrm{{cm}},\;" + rf"d_{{rx}}={cfg['d_rx_cm']}\,\mathrm{{cm}}$", + fontsize=14, + ) + ax.set_xlabel("Distance (m)", fontsize=fig_params['fontsize']) + ax.set_ylabel(r"Singular-value ratio $\sigma_2/\sigma_1$", + fontsize=fig_params['fontsize']) + ax.set_ylim(0, 1.05) + ax.tick_params(labelsize=fig_params['tick_labelsize']) + ax.legend(fontsize=12, framealpha=0.5, loc='upper right') + ax.grid(True) + fig.tight_layout() + out = os.path.join(args.save_dir, f"fig10_{case_name}_validation.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out}") + + if args.quick: + break # one plot is enough for smoke test + + +if __name__ == "__main__": + main() diff --git a/scripts/fig11_2d_beampattern.py b/scripts/fig11_2d_beampattern.py new file mode 100644 index 0000000..92995d7 --- /dev/null +++ b/scripts/fig11_2d_beampattern.py @@ -0,0 +1,295 @@ +""" +Fig 11 — 2D beam patterns (transverse × range heat maps). + +Reproduces Fig. 11 from the ArrayLink paper (INFOCOM 2026). + +The satellite is placed on a Cartesian (X, Y) grid: + X : transverse offset from boresight (km) + Y : range / boresight distance (km) + + satellite position = (X, 0, Y) + +DC weights are focused at (0, 0, 500) km (zenith, 500 km range). + +Shows: + (a) Monolithic UPA 128×128 — far-field gain ridge, constant in range + (b) ArrayLink 16×(32×32) — near-field focal spot at Y=500 km + +ArrayLink panel placement (--placement) +---------------------------------------- + random [default] 16 panels placed uniformly at random within the + aperture (Lx × Ly). Different each run unless + --seed is fixed. + center-dense Power-law concentrated toward boresight (gamma=3), + as used in the paper. + uniform Regular 4×4 grid spanning the full aperture. + +Resolution (default: fast preview, ~30 s) +------------------------------------------ + Default 20 km X steps × 25 km Y steps (301 × 119 ≈ 36 K pts) + --full 5 km X steps × 20 km Y steps (1201 × 148 ≈ 178 K pts, ~5 min) + +Usage +----- + # Fast preview with random panel placement (default): + python scripts/fig11_2d_beampattern.py + + # Full-resolution, paper-exact center-dense layout: + python scripts/fig11_2d_beampattern.py --placement center-dense --full + + # Reproducible random layout, full resolution: + python scripts/fig11_2d_beampattern.py --placement random --seed 42 --full + + # Smoke test (tiny grid, no plot saved): + python scripts/fig11_2d_beampattern.py --quick + + # Interactive Plotly figure (requires plotly): + python scripts/fig11_2d_beampattern.py --interactive +""" +import argparse +import gc +import os +import sys +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +import numpy as np +import matplotlib.pyplot as plt + +try: + import scienceplots + plt.style.use(['science', 'no-latex']) +except ImportError: + pass + +from arraylink.array_geometry import build_ground_station, upa_positions, place_subarrays +from arraylink.beamforming import ( + dc_weights, compute_beam_pattern_numpy, +) + + +def parse_args(): + ap = argparse.ArgumentParser( + description="Fig 11: 2D beam pattern heat map", + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + ap.add_argument("--quick", action="store_true", + help="Tiny grid for smoke test (skips plot save)") + ap.add_argument("--full", action="store_true", + help="Full resolution: 5 km × 20 km grid (~178 K pts, ~5 min). " + "Default is fast preview: 20 km × 25 km grid (~36 K pts, ~30 s).") + ap.add_argument("--placement", default="random", + choices=["random", "center-dense", "uniform"], + help="ArrayLink panel placement strategy (default: random). " + "'random' places 16 panels uniformly at random within the aperture. " + "'center-dense' uses the paper's power-law layout (gamma=3). " + "'uniform' uses a regular 4×4 grid.") + ap.add_argument("--seed", type=int, default=0, + help="RNG seed for random placement (default: 0)") + ap.add_argument("--interactive", action="store_true", + help="Open Plotly interactive figure (requires plotly)") + ap.add_argument("--save-dir", default="paper_figures") + return ap.parse_args() + + +ELEMENT_GAIN_DBI = 6.0 +TARGET_R_KM = 500.0 # DC focus range (km) — matches paper Fig. 11 + +# ArrayLink aperture dimensions (km) +AL_LX = np.sqrt(2) # 1.414 km along X +AL_LY = 1.0 # 1.000 km along Y +AL_NX = 4 # panels along X +AL_NY = 4 # panels along Y + + +def build_arraylink(placement, subarray_shape, element_spacing, seed): + """ + Build the ArrayLink antenna array according to the chosen placement strategy. + + Parameters + ---------- + placement : 'random' | 'center-dense' | 'uniform' + subarray_shape : (M, N) elements per panel + element_spacing: element pitch within a panel (km) + seed : RNG seed (used for 'random' placement) + + Returns + ------- + ants : (Nx*Ny*M*N, 3) all antenna positions + bases : (Nx*Ny, 3) panel base positions + """ + n_panels = AL_NX * AL_NY + + if placement == "random": + rng = np.random.default_rng(seed) + bx = rng.uniform(-AL_LX / 2.0, AL_LX / 2.0, n_panels) + by = rng.uniform(-AL_LY / 2.0, AL_LY / 2.0, n_panels) + bases = np.stack([bx, by, np.zeros(n_panels)], axis=1) # (16, 3) + + elif placement == "center-dense": + _, bases = build_ground_station( + mode='arraylink', subarray_shape=subarray_shape, + element_spacing=element_spacing, + Nx=AL_NX, Ny=AL_NY, Lx=AL_LX, Ly=AL_LY, + gamma=3.0, seed=seed, + ) + + elif placement == "uniform": + _, bases = build_ground_station( + mode='uniform', subarray_shape=subarray_shape, + element_spacing=element_spacing, + Nx=AL_NX, Ny=AL_NY, Lx=AL_LX, Ly=AL_LY, + ) + + else: + raise ValueError(f"Unknown placement '{placement}'") + + ants = place_subarrays(bases, subarray_shape, element_spacing) + return ants, bases + + +def compute_2d_pattern(gnd_ants, weights, lam_km, x_axis, y_axis, quick, + batch_size=2000): + """ + Compute beam pattern on a Cartesian (X, Y) grid. + + Satellite is placed at (X, 0, Y) km. Returns (Ny, Nx) array of gain in dBi. + Uses batched numpy evaluation (~400 MB RAM per batch). + """ + Nx, Ny = len(x_axis), len(y_axis) + XX, YY = np.meshgrid(x_axis, y_axis, indexing='ij') # (Nx, Ny) + flat_pts = np.stack([XX.ravel(), np.zeros(Nx * Ny), YY.ravel()], axis=-1) + + if quick: + bp = compute_beam_pattern_numpy(flat_pts, gnd_ants, weights, lam_km) + bp += ELEMENT_GAIN_DBI + return bp.reshape(Nx, Ny).T # (Ny, Nx) + + # Batched numpy: process batch_size points at a time to keep RAM bounded + N = len(flat_pts) + results = np.empty(N, dtype=np.float64) + for start in range(0, N, batch_size): + end = min(start + batch_size, N) + results[start:end] = compute_beam_pattern_numpy( + flat_pts[start:end], gnd_ants, weights, lam_km + ) + if start % (batch_size * 20) == 0 and start > 0: + pct = 100 * start / N + print(f" {pct:.0f}% ({start:,}/{N:,})", flush=True) + gc.collect() + return (results + ELEMENT_GAIN_DBI).reshape(Nx, Ny).T # (Ny, Nx) + + +def main(): + args = parse_args() + os.makedirs(args.save_dir, exist_ok=True) + + F_HZ = 28e9 + LAM_KM = 3e8 / F_HZ / 1e3 + + if args.quick: + # Minimal grid for smoke test — just verify code paths run + x_axis = np.linspace(-3000, 3000, 13) + y_axis = np.linspace(50, 3000, 7) + elif args.full: + # Full resolution — resolves UPA sidelobes and ArrayLink focal depth + # 1201 × 148 ≈ 178 K pts → ~5 min on a single CPU core + x_axis = np.arange(-3000, 3001, 5, dtype=float) # 1201 pts, 5 km steps + y_axis = np.arange(50, 3001, 20, dtype=float) # 148 pts, 20 km steps + else: + # Fast preview — good enough to see main beam structure and focal spot + # 301 × 119 ≈ 36 K pts → ~30 s on a single CPU core + x_axis = np.arange(-3000, 3001, 20, dtype=float) # 301 pts, 20 km steps + y_axis = np.arange(50, 3001, 25, dtype=float) # 119 pts, 25 km steps + + sat_loc = np.array([0.0, 0.0, TARGET_R_KM]) + + # ---------- Build antenna arrays ---------- + # UPA 128×128 — compact monolithic array, λ/2 spacing + upa_ants = upa_positions(128, 128, element_spacing=LAM_KM / 2.0) + + # ArrayLink 4×4 panels × 32×32 elements — placement set by --placement flag + al_ants, al_bases = build_arraylink( + placement=args.placement, + subarray_shape=(32, 32), + element_spacing=LAM_KM / 2.0, + seed=args.seed, + ) + + w_upa = dc_weights(sat_loc, upa_ants, LAM_KM) + w_al = dc_weights(sat_loc, al_ants, LAM_KM) + + res_tag = "full" if args.full else ("quick" if args.quick else "fast") + print(f"UPA: {len(upa_ants)} elements | ArrayLink: {len(al_ants)} elements " + f"[placement={args.placement}, seed={args.seed}]") + print(f"Grid: {len(x_axis)} × {len(y_axis)} = {len(x_axis)*len(y_axis):,} points " + f"[resolution={res_tag}]") + + print("Computing UPA beam pattern …", flush=True) + bp_upa = compute_2d_pattern(upa_ants, w_upa, LAM_KM, x_axis, y_axis, args.quick) + + print("Computing ArrayLink beam pattern …", flush=True) + bp_al = compute_2d_pattern(al_ants, w_al, LAM_KM, x_axis, y_axis, args.quick) + + if args.quick: + print("Quick mode: patterns computed, skipping plot save.") + print(f" UPA max gain : {bp_upa.max():.1f} dBi") + print(f" ArrayLink max : {bp_al.max():.1f} dBi") + return + + peak = max(bp_upa.max(), bp_al.max()) + vmin = ELEMENT_GAIN_DBI # floor = single-element gain (matches paper colorbar) + vmax = peak + print(f"Peak = {peak:.2f} dBi | colorbar [{vmin:.1f}, {vmax:.1f}] dBi") + + # ---------- Static matplotlib figure ---------- + al_title = { + "random": f"(b) ArrayLink — random (seed={args.seed})", + "center-dense": "(b) ArrayLink — center-dense", + "uniform": "(b) ArrayLink — uniform", + }[args.placement] + + X, Y = np.meshgrid(x_axis, y_axis) # both (Ny, Nx) + fig, axes = plt.subplots(1, 2, figsize=(13, 5.5)) + for ax, bp, title in [ + (axes[0], bp_upa, "(a) Monolithic UPA"), + (axes[1], bp_al, al_title), + ]: + c = ax.pcolormesh(X, Y, bp, shading='auto', + cmap='inferno', vmin=vmin, vmax=vmax) + ax.axhline(TARGET_R_KM, color='white', lw=0.6, ls='--', alpha=0.45) + ax.set_title(title, fontsize=13) + ax.set_xlabel("X (km)", fontsize=12) + ax.set_ylabel("Y (km)", fontsize=12) + fig.colorbar(c, ax=ax, label="Gain (dB)", shrink=0.88) + fig.tight_layout() + out = os.path.join(args.save_dir, "fig11_2d_beampattern.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out}") + + # ---------- Optional Plotly interactive ---------- + if args.interactive: + try: + import plotly.graph_objects as go + from plotly.subplots import make_subplots + fig_p = make_subplots(rows=1, cols=2, + subplot_titles=["Monolithic UPA", al_title]) + for col, (bp, lbl) in enumerate([(bp_upa, "UPA"), (bp_al, "AL")], start=1): + fig_p.add_trace( + go.Heatmap( + x=x_axis, y=y_axis, z=bp, + colorscale='Hot', zmin=vmin, zmax=vmax, + colorbar=dict(title="dB", x=col * 0.5), + showscale=(col == 2), + ), + row=1, col=col, + ) + fig_p.update_xaxes(title_text="X (km)", row=1, col=col) + fig_p.update_yaxes(title_text="Y (km)", row=1, col=col) + fig_p.update_layout(title="2D Beam Patterns: UPA vs ArrayLink") + fig_p.show() + except ImportError: + print("Plotly not installed. Install with: mamba install -c conda-forge plotly") + + +if __name__ == "__main__": + main() diff --git a/scripts/fig12_mimo_dof.py b/scripts/fig12_mimo_dof.py new file mode 100644 index 0000000..f84b5aa --- /dev/null +++ b/scripts/fig12_mimo_dof.py @@ -0,0 +1,158 @@ +""" +Fig 12 — LoS MIMO degrees of freedom vs distance for ArrayLink. + +Reproduces Fig. 12 from the ArrayLink paper (INFOCOM 2026). + +Computes the 4x4 MIMO channel between a satellite 4-element array +(2x2 panels, 1.414m x 1m grid) and the ArrayLink ground station +(16 x 32x32 panels over 1.414km x 1km), then plots: + (a) Normalised singular-value ratios sigma_k/sigma_1 vs distance + (b) Effective number of spatial streams (DoF) vs distance + +Usage +----- + python scripts/fig12_mimo_dof.py [--quick] [--save-dir paper_figures] +""" +import argparse +import os +import sys +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +import numpy as np +import matplotlib.pyplot as plt + +try: + import scienceplots + plt.style.use(['science', 'no-latex']) +except ImportError: + pass + +from arraylink.channel import compute_channel_matrix, compute_singular_values +from arraylink.array_geometry import build_ground_station, upa_positions +from arraylink.utils import spherical2cartesian + + +def parse_args(): + ap = argparse.ArgumentParser(description="Fig 12: MIMO DoF vs distance") + ap.add_argument("--quick", action="store_true", + help="Run on a coarse distance grid (fast, for smoke test)") + ap.add_argument("--save-dir", default="paper_figures") + return ap.parse_args() + + +def main(): + args = parse_args() + os.makedirs(args.save_dir, exist_ok=True) + + # --- System parameters --- + F_HZ = 28e9 + C = 3e8 + LAM_KM = C / F_HZ / 1e3 # wavelength in km (~1.07e-5 km) + TAU = 0.1 # feasibility threshold + + # Ground station: 16 x 32x32 panels over sqrt(2)km x 1km + GND_NX, GND_NY = 4, 4 + GND_LX, GND_LY = np.sqrt(2), np.sqrt(2) # km + PANEL_SHAPE = (32, 32) + + # Satellite: 2x2 element array, 1.414m x 1m aperture + # Assumption: 4-element array arranged as 2x2 UPA on a 1.414m x 1m grid + # in km: 1.414m = 1.414e-3 km, 1m = 1e-3 km + SAT_DX_KM = 1.414e-3 # 1.414 m in km + SAT_DY_KM = 1.0e-3 # 1.0 m in km + + if args.quick: + dist_array_km = np.linspace(100, 3000, 5) + else: + dist_array_km = np.arange(100, 3001, 50, dtype=float) + + # Build ground station once (positions don't change with distance) + gnd_ants, gnd_bases = build_ground_station( + mode='arraylink', + subarray_shape=PANEL_SHAPE, + element_spacing=LAM_KM / 2.0, + Nx=GND_NX, Ny=GND_NY, + Lx=GND_LX, Ly=GND_LY, + gamma=3.0, seed=0, + ) + + # 2x2 satellite array template (centred at origin, will be translated) + sat_template = upa_positions(2, 2, element_spacing=SAT_DX_KM / 2.0) + # The above gives spacing SAT_DX_KM/2 between elements; adjust so the + # max separation is SAT_DX_KM x SAT_DY_KM + span_x = np.ptp(sat_template[:, 0]) + span_y = np.ptp(sat_template[:, 1]) + sat_template[:, 0] *= (SAT_DX_KM / span_x) if span_x > 0 else 1 + sat_template[:, 1] *= (SAT_DY_KM / span_y) if span_y > 0 else 1 + + # Collect singular-value ratios vs distance + # sigma_k / sigma_1 for k=1,2,3 (paper Fig 12a uses 3 ratios beyond sigma_1) + sing_ratios = {1: [], 2: [], 3: []} + dof_list = [] + + for d_km in dist_array_km: + # Place satellite directly above (theta=0, phi=0) at distance d_km + sat_center = np.array([0.0, 0.0, d_km]) + sat_coords = sat_template + sat_center + + H = compute_channel_matrix(sat_coords, gnd_ants, LAM_KM) + s = compute_singular_values(H, normalise=True) + + # sigma_k / sigma_0 for k=1,2,3 + for k in [1, 2, 3]: + ratio = float(s[k] / s[0]) if len(s) > k and s[0] > 0 else 0.0 + sing_ratios[k].append(ratio) + + dof = int(np.sum(s >= TAU)) + dof_list.append(dof) + + if args.quick: + print("DoF at each distance:", list(zip(dist_array_km.astype(int), dof_list))) + return + + fig_params = dict(fontsize=16, tick_labelsize=14) + colors = ['steelblue', 'darkorange', 'green'] + labels = [r'$\sigma_2/\sigma_1$', r'$\sigma_3/\sigma_1$', r'$\sigma_4/\sigma_1$'] + + # --- Fig 12a: singular-value ratios --- + fig, ax = plt.subplots(figsize=(6, 5)) + for idx, k in enumerate([1, 2, 3]): + ax.plot(dist_array_km, sing_ratios[k], linewidth=2, + color=colors[idx], label=labels[idx]) + ax.axhline(TAU, linestyle='--', color='firebrick', linewidth=1.5, + label=f'Threshold (τ={TAU})') + ax.set_xlabel("Distance (km)", fontsize=fig_params['fontsize']) + ax.set_ylabel(r"Singular-value ratio $\sigma_k/\sigma_1$", + fontsize=fig_params['fontsize']) + ax.set_ylim(0, 1.05) + ax.set_xlim(0, 3000) + ax.tick_params(labelsize=fig_params['tick_labelsize']) + ax.legend(fontsize=13, framealpha=0.5) + ax.grid(True) + fig.tight_layout() + out = os.path.join(args.save_dir, "fig12a_singular_value_ratios.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out}") + + # --- Fig 12b: DoF vs distance --- + fig, ax = plt.subplots(figsize=(6, 5)) + ax.plot(dist_array_km, dof_list, linewidth=2, color='steelblue', label='DoFs') + ax.axhline(TAU, linestyle='--', color='firebrick', linewidth=1.5, + label=f'Threshold (τ={TAU})', alpha=0) # hidden, just for legend alignment + ax.set_xlabel("Distance (km)", fontsize=fig_params['fontsize']) + ax.set_ylabel(r"DoF (prominent $\sigma$ values)", fontsize=fig_params['fontsize']) + ax.set_ylim(0.5, 4.5) + ax.set_xlim(0, 3000) + ax.set_yticks([1, 2, 3, 4]) + ax.tick_params(labelsize=fig_params['tick_labelsize']) + ax.legend(fontsize=13, framealpha=0.5) + ax.grid(True) + fig.tight_layout() + out = os.path.join(args.save_dir, "fig12b_dof_vs_distance.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out}") + + +if __name__ == "__main__": + main() diff --git a/scripts/fig13_throughput.py b/scripts/fig13_throughput.py new file mode 100644 index 0000000..83c625d --- /dev/null +++ b/scripts/fig13_throughput.py @@ -0,0 +1,131 @@ +""" +Fig 13 — Per-link and aggregate throughput comparison. + +Reproduces Fig. 13 from the ArrayLink paper (INFOCOM 2026). + +Compares: + - 1 stream @ 52.6 dBi (1.85 m dish baseline) + - 2 streams @ 48.14 dBi (ArrayLink) + - 4 streams @ 48.14 dBi (ArrayLink) + +Two models: + (a) Per-link throughput: each configuration uses total bandwidth B + (b) Aggregate throughput: each stream adds B_single Hz of bandwidth + +Usage +----- + python scripts/fig13_throughput.py [--quick] [--save-dir paper_figures] +""" +import argparse +import os +import sys +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +import numpy as np +import matplotlib.pyplot as plt + +try: + import scienceplots + plt.style.use(['science', 'no-latex']) +except ImportError: + pass + +from arraylink.beamforming import parabolic_gain_dbi, total_array_gain_dbi + + +def parse_args(): + ap = argparse.ArgumentParser(description="Fig 13: throughput comparison") + ap.add_argument("--quick", action="store_true") + ap.add_argument("--save-dir", default="paper_figures") + return ap.parse_args() + + +def spectral_eff(snr_db, gap_db=0.0): + """Shannon spectral efficiency (bps/Hz) with implementation gap.""" + return np.log2(1.0 + 10 ** ((snr_db - gap_db) / 10.0)) + + +def main(): + args = parse_args() + os.makedirs(args.save_dir, exist_ok=True) + + # --- Parameters --- + G_DISH_185 = parabolic_gain_dbi(1.85, 28e9, efficiency=0.6) # ~52.6 dBi (baseline) + G_ARRAYLINK = total_array_gain_dbi(16, 36.1) # ~48.14 dBi + B_SINGLE_HZ = 250e6 # 250 MHz per link + GAP_DB = 2.0 # implementation gap to Shannon + + # SNR baseline: defined for the 1-stream dish case + if args.quick: + SNR0_SWEEP = np.array([10.0, 15.0, 20.0]) + else: + SNR0_SWEEP = np.arange(5, 26, 2.5) + + configs = [ + {"name": "1 stream @ 52.6 dBi (dish)", "N": 1, "gain_dbi": G_DISH_185}, + {"name": "2 streams @ 48.14 dBi (AL)", "N": 2, "gain_dbi": G_ARRAYLINK}, + {"name": "4 streams @ 48.14 dBi (AL)", "N": 4, "gain_dbi": G_ARRAYLINK}, + ] + markers = ['o', 's', '^'] + lines = ['--', '-', '-.'] + colors = ['steelblue', 'darkorange', 'seagreen'] + + # Pre-compute throughput curves + perlink_gbps = {c["name"]: [] for c in configs} + agg_gbps = {c["name"]: [] for c in configs} + + for snr0 in SNR0_SWEEP: + for cfg in configs: + dg = cfg["gain_dbi"] - G_DISH_185 # gain delta vs dish + snr_db = snr0 + dg # per-link SNR + eta = spectral_eff(snr_db, gap_db=GAP_DB) # bps/Hz + # per-link: total BW / N streams shared + perlink_gbps[cfg["name"]].append(eta * B_SINGLE_HZ / 1e9) + # aggregate: each stream gets B_single + agg_gbps[cfg["name"]].append(cfg["N"] * eta * B_SINGLE_HZ / 1e9) + + if args.quick: + for cfg in configs: + print(f"{cfg['name']}: per-link={perlink_gbps[cfg['name']]}, agg={agg_gbps[cfg['name']]}") + return + + fig_params = dict(fontsize=20, tick_labelsize=18) + + # --- Fig 13a: per-link throughput --- + fig, ax = plt.subplots(figsize=(6, 5)) + for i, cfg in enumerate(configs): + ax.plot(SNR0_SWEEP, perlink_gbps[cfg["name"]], marker=markers[i], + linestyle=lines[i], color=colors[i], linewidth=2.5, markersize=6, + label=cfg["name"]) + ax.set_xlabel(f"Baseline SNR₀ for {G_DISH_185:.1f} dBi (dB)", + fontsize=fig_params['fontsize']) + ax.set_ylabel("Per-link throughput (Gbps)", fontsize=fig_params['fontsize']) + ax.tick_params(labelsize=fig_params['tick_labelsize']) + ax.legend(fontsize=14, framealpha=1, loc='upper left') + ax.grid(True) + fig.tight_layout() + out = os.path.join(args.save_dir, "fig13a_perlink_throughput.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out}") + + # --- Fig 13b: aggregate throughput --- + fig, ax = plt.subplots(figsize=(6, 5)) + for i, cfg in enumerate(configs): + ax.plot(SNR0_SWEEP, agg_gbps[cfg["name"]], marker=markers[i], + linestyle=lines[i], color=colors[i], linewidth=2.5, markersize=6, + label=cfg["name"]) + ax.set_xlabel(f"Baseline SNR₀ for {G_DISH_185:.1f} dBi (dB)", + fontsize=fig_params['fontsize']) + ax.set_ylabel("Aggregate throughput (Gbps)", fontsize=fig_params['fontsize']) + ax.tick_params(labelsize=fig_params['tick_labelsize']) + ax.legend(fontsize=14, framealpha=1, loc='upper left') + ax.grid(True) + fig.tight_layout() + out = os.path.join(args.save_dir, "fig13b_aggregate_throughput.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out}") + + +if __name__ == "__main__": + main() diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_array_geometry.py b/tests/test_array_geometry.py new file mode 100644 index 0000000..1d106e9 --- /dev/null +++ b/tests/test_array_geometry.py @@ -0,0 +1,122 @@ +"""Unit tests for arraylink.array_geometry.""" +import numpy as np +import pytest +from arraylink.array_geometry import ( + upa_positions, + center_dense_positions, + uniform_grid_positions, + place_subarrays, + build_ground_station, +) + + +class TestUPAPositions: + def test_shape(self): + p = upa_positions(4, 8, 0.5) + assert p.shape == (32, 3) + + def test_centred_at_origin(self): + p = upa_positions(4, 4, 1.0) + centroid = p.mean(axis=0) + np.testing.assert_allclose(centroid, [0, 0, 0], atol=1e-10) + + def test_spacing(self): + spacing = 0.5 + p = upa_positions(3, 1, spacing) + # x coordinates should be -spacing, 0, +spacing + xs = np.sort(p[:, 0]) + np.testing.assert_allclose(np.diff(xs), spacing, atol=1e-10) + + def test_z_is_zero(self): + p = upa_positions(4, 4, 0.5) + assert np.all(p[:, 2] == 0.0) + + def test_center_offset(self): + center = [1.0, 2.0, 3.0] + p = upa_positions(2, 2, 1.0, center=center) + centroid = p.mean(axis=0) + np.testing.assert_allclose(centroid, center, atol=1e-10) + + +class TestCenterDensePositions: + def test_shape(self): + b = center_dense_positions(4, 4, 1.0, 1.0) + assert b.shape == (16, 3) + + def test_within_aperture(self): + Gx, Gy = 1.414, 1.0 + b = center_dense_positions(4, 4, Gx, Gy) + assert np.all(b[:, 0] >= -Gx / 2 - 1e-9) + assert np.all(b[:, 0] <= Gx / 2 + 1e-9) + assert np.all(b[:, 1] >= -Gy / 2 - 1e-9) + assert np.all(b[:, 1] <= Gy / 2 + 1e-9) + + def test_more_concentrated_than_uniform(self): + """Center-dense should have more panels near origin than uniform grid.""" + Gx, Gy = 1.0, 1.0 + cd = center_dense_positions(4, 4, Gx, Gy, gamma=3.0) + ug = uniform_grid_positions(4, 4, Gx, Gy) + half_r = 0.25 # within central 25% radius + cd_inner = np.sum(np.linalg.norm(cd[:, :2], axis=1) < half_r) + ug_inner = np.sum(np.linalg.norm(ug[:, :2], axis=1) < half_r) + assert cd_inner >= ug_inner, "Center-dense should concentrate more near origin" + + def test_reproducible_with_seed(self): + b1 = center_dense_positions(4, 4, 1.0, 1.0, jitter_frac=0.1, seed=42) + b2 = center_dense_positions(4, 4, 1.0, 1.0, jitter_frac=0.1, seed=42) + np.testing.assert_array_equal(b1, b2) + + def test_different_seeds_differ(self): + b1 = center_dense_positions(4, 4, 1.0, 1.0, jitter_frac=0.2, seed=1) + b2 = center_dense_positions(4, 4, 1.0, 1.0, jitter_frac=0.2, seed=2) + assert not np.allclose(b1, b2) + + +class TestUniformGridPositions: + def test_shape(self): + g = uniform_grid_positions(4, 4, 1.0, 1.0) + assert g.shape == (16, 3) + + def test_centred_at_origin(self): + g = uniform_grid_positions(4, 4, 1.0, 1.0) + centroid = g.mean(axis=0) + np.testing.assert_allclose(centroid[:2], [0, 0], atol=1e-10) + + +class TestPlaceSubarrays: + def test_total_count(self): + bases = np.zeros((4, 3)) + result = place_subarrays(bases, subarray_shape=(8, 8), element_spacing=0.005) + assert result.shape == (4 * 64, 3) + + def test_single_base_centred(self): + base = np.array([[0.0, 0.0, 0.0]]) + result = place_subarrays(base, subarray_shape=(4, 4), element_spacing=1.0) + centroid = result.mean(axis=0) + np.testing.assert_allclose(centroid, [0, 0, 0], atol=1e-10) + + +class TestBuildGroundStation: + def test_arraylink_mode(self): + ants, bases = build_ground_station( + mode='arraylink', + subarray_shape=(32, 32), + element_spacing=0.005, + Nx=4, Ny=4, Lx=1.414, Ly=1.0, + ) + assert bases.shape == (16, 3) + assert ants.shape == (16 * 1024, 3) + + def test_uniform_mode(self): + ants, bases = build_ground_station( + mode='uniform', + subarray_shape=(4, 4), + element_spacing=0.005, + Nx=2, Ny=2, Lx=0.5, Ly=0.5, + ) + assert bases.shape == (4, 3) + assert ants.shape == (64, 3) + + def test_invalid_mode(self): + with pytest.raises(ValueError): + build_ground_station('invalid', (4, 4), 0.005, 2, 2, 0.5, 0.5) diff --git a/tests/test_beamforming.py b/tests/test_beamforming.py new file mode 100644 index 0000000..355a55f --- /dev/null +++ b/tests/test_beamforming.py @@ -0,0 +1,117 @@ +"""Unit tests for arraylink.beamforming.""" +import numpy as np +import pytest +from arraylink.beamforming import ( + parabolic_gain_dbi, + total_array_gain_dbi, + marginal_gain_dbi, + dc_weights, + compute_beam_pattern_numpy, +) + + +class TestParabolicGain: + def test_1p47m_dish_28ghz(self): + """1.47 m dish at 28 GHz, eta=0.6 → ~49.5 dBi (paper Sec. II-A).""" + g = parabolic_gain_dbi(1.47, 28e9, efficiency=0.6) + assert 48.0 < g < 51.0, f"Expected ~49.5 dBi, got {g:.2f}" + + def test_1p85m_dish_28ghz(self): + """1.85 m dish at 28 GHz → ~52.6 dBi (paper Sec. II-A).""" + g = parabolic_gain_dbi(1.85, 28e9, efficiency=0.6) + assert 51.0 < g < 54.0, f"Expected ~52.6 dBi, got {g:.2f}" + + def test_larger_dish_more_gain(self): + g1 = parabolic_gain_dbi(1.0, 28e9) + g2 = parabolic_gain_dbi(2.0, 28e9) + assert g2 > g1 + + +class TestArrayGain: + def test_single_panel_equals_pa_gain(self): + G_pa = 36.1 + g = total_array_gain_dbi(1, G_pa) + assert abs(g - G_pa) < 1e-10 + + def test_16_panels_formula(self): + """Paper example: 16 panels × 36.1 dBi -> ~48.1 dBi.""" + g = total_array_gain_dbi(16, 36.1) + expected = 10 * np.log10(16) + 36.1 + assert abs(g - expected) < 1e-10 + assert 47.0 < g < 49.0 + + def test_marginal_gain_decreasing(self): + N = np.arange(1, 50) + mg = marginal_gain_dbi(N, 36.1) + # Marginal gain should be strictly decreasing + assert np.all(np.diff(mg) < 0) + + def test_marginal_gain_below_0p25_after_16(self): + """ + Paper Fig 4b: marginal gain 'beyond 16 panels' falls below 0.25 dB. + + At N=16, G(17)-G(16) ≈ 0.263 dB (slightly above 0.25). + The threshold is crossed near N=18: G(19)-G(18) ≈ 0.23 dB. + We verify the transition happens within a few panels of N=16. + """ + N_check = np.arange(16, 25) + mg = marginal_gain_dbi(N_check, 36.1) + assert np.any(mg < 0.25), ( + f"Expected marginal gain to drop below 0.25 dB near N=16-20, " + f"but values were {mg.round(4)}" + ) + + +class TestDCWeights: + def test_shape(self): + sat = np.array([0.0, 0.0, 500.0]) + rx = np.random.randn(64, 3) + w = dc_weights(sat, rx, wavelength=0.01 / 1e3) + assert w.shape == (64, 1) + + def test_normalised(self): + sat = np.array([0.0, 0.0, 500.0]) + rx = np.random.randn(32, 3) + w = dc_weights(sat, rx, wavelength=0.01 / 1e3) + assert abs(np.linalg.norm(w) - 1.0) < 1e-10 + + def test_complex(self): + sat = np.array([0.0, 0.0, 1.0]) + rx = np.array([[0.0, 0.0, 0.0], [0.5, 0.0, 0.0]]) + w = dc_weights(sat, rx, wavelength=0.01) + assert np.iscomplexobj(w) + + def test_boresight_peak(self): + """Beam should peak at the steered direction.""" + lam = 0.01 + # Single receive antenna at origin, satellite at (0,0,1) km + rx = np.array([[0.0, 0.0, 0.0]]) + sat_target = np.array([0.0, 0.0, 1.0]) + sat_off = np.array([1.0, 0.0, 1.0]) + + w = dc_weights(sat_target, rx, lam) + + bp_on = compute_beam_pattern_numpy( + sat_target.reshape(1, 3), rx, w, lam + ) + bp_off = compute_beam_pattern_numpy( + sat_off.reshape(1, 3), rx, w, lam + ) + assert bp_on[0] >= bp_off[0] + + +class TestBeamPatternNumpy: + def test_shape(self): + sat_pts = np.random.randn(10, 3) + np.array([0, 0, 5]) + rx = np.random.randn(8, 3) + w = np.ones((8, 1), dtype=complex) / np.sqrt(8) + bp = compute_beam_pattern_numpy(sat_pts, rx, w, wavelength=0.01) + assert bp.shape == (10,) + + def test_values_are_db(self): + sat_pts = np.array([[0.0, 0.0, 1.0]]) + rx = np.zeros((4, 3)) + w = np.ones((4, 1), dtype=complex) / 2.0 + bp = compute_beam_pattern_numpy(sat_pts, rx, w, wavelength=0.01) + # Result should be finite (not NaN, not inf) + assert np.all(np.isfinite(bp)) diff --git a/tests/test_channel.py b/tests/test_channel.py new file mode 100644 index 0000000..8c7732f --- /dev/null +++ b/tests/test_channel.py @@ -0,0 +1,128 @@ +"""Unit tests for arraylink.channel.""" +import numpy as np +import pytest +from arraylink.channel import ( + compute_channel_matrix, + compute_singular_values, + singular_value_ratio, + degrees_of_freedom, + mimo_region_bounds, + theoretical_singular_value_ratio, + spectral_efficiency, +) + + +class TestChannelMatrix: + def test_shape(self): + tx = np.array([[0.0, 0.0, 0.0], [0.1, 0.0, 0.0]]) + rx = np.array([[0.0, 0.0, 1.0], [0.0, 0.2, 1.0], [0.0, -0.2, 1.0]]) + lam = 0.01 + H = compute_channel_matrix(tx, rx, lam) + assert H.shape == (3, 2) + + def test_complex(self): + tx = np.array([[0.0, 0.0, 0.0]]) + rx = np.array([[0.0, 0.0, 1.0]]) + H = compute_channel_matrix(tx, rx, 0.01) + assert np.iscomplexobj(H) + + def test_path_loss_decreases_with_distance(self): + """Amplitude should decrease as distance doubles.""" + tx = np.array([[0.0, 0.0, 0.0]]) + lam = 0.01 + rx1 = np.array([[0.0, 0.0, 1.0]]) + rx2 = np.array([[0.0, 0.0, 2.0]]) + H1 = compute_channel_matrix(tx, rx1, lam) + H2 = compute_channel_matrix(tx, rx2, lam) + assert np.abs(H1[0, 0]) > np.abs(H2[0, 0]) + + def test_2x2_singular_values_analytical(self): + """ + For a symmetric 2×2 LoS geometry at r = d_tx*d_rx/lambda (Delta=2pi), + the two singular values should be equal (well-conditioned channel). + At r_max the ratio should be close to tau=0.1. + """ + lam = 0.01 # 10 mm (28 GHz in metres) + d_tx = 0.2 + d_rx = 0.2 + _, r_max = mimo_region_bounds(d_tx, d_rx, lam, tau=0.1) + # At r_max, ratio should be ~tau + ratio = theoretical_singular_value_ratio(r_max, d_tx, d_rx, lam) + assert abs(ratio - 0.1) < 0.02, f"Expected ~0.1 got {ratio:.4f}" + + +class TestSingularValues: + def test_normalised_norm(self): + tx = np.array([[0.0, 0.0, 0.0], [0.5, 0.0, 0.0]]) + rx = np.array([[0.0, 0.0, 1.0], [0.0, 0.5, 1.0]]) + H = compute_channel_matrix(tx, rx, 0.01) + s = compute_singular_values(H, normalise=True) + assert abs(np.linalg.norm(s) - 1.0) < 1e-10 + + def test_descending_order(self): + tx = np.random.randn(4, 3) + rx = np.random.randn(4, 3) + np.array([0, 0, 10]) + H = compute_channel_matrix(tx, rx, 0.01) + s = compute_singular_values(H) + assert np.all(np.diff(s) <= 0), "Singular values not in descending order" + + +class TestMIMOBounds: + def test_paper_example_metres(self): + """ + Paper Sec. III-C example: d_tx=d_rx=0.2m, lambda=0.01m, tau=0.1 + -> r_max ~ 62 m. + """ + r_min, r_max = mimo_region_bounds(0.2, 0.2, 0.01, tau=0.1) + assert 55 < r_max < 70, f"r_max={r_max:.1f} m, expected ~62 m" + assert r_min < r_max + + def test_satellite_scale(self): + """ + Paper Sec. III-C: d_tx=2km, d_rx=1m, lambda=0.01m, tau=0.1 + -> r_max in the thousands-of-km range. + + Note: The paper quotes "≈2500 km" as a rounded figure; the closed-form + Eq. (11) with these exact parameters yields ~3152 km. We verify the + order-of-magnitude correctness (satellite-scale, not lab-scale). + """ + lam_km = 0.01 / 1e3 + d_tx_km = 2.0 # 2 km ground aperture + d_rx_km = 1e-3 # 1 m satellite aperture in km + _, r_max = mimo_region_bounds(d_tx_km, d_rx_km, lam_km, tau=0.1) + assert 2000 < r_max < 4000, ( + f"r_max={r_max:.0f} km — expected satellite-scale range (2000-4000 km)" + ) + + def test_larger_aperture_extends_range(self): + lam = 0.01 + _, r_max_small = mimo_region_bounds(0.2, 0.2, lam) + _, r_max_large = mimo_region_bounds(1.0, 0.2, lam) + assert r_max_large > r_max_small + + +class TestSpectralEfficiency: + def test_positive(self): + tx = np.array([[0.0, 0.0, 0.0], [0.2, 0.0, 0.0]]) + rx = np.array([[0.0, 0.0, 10.0], [0.0, 0.2, 10.0]]) + H = compute_channel_matrix(tx, rx, 0.01) + C = spectral_efficiency(H, snr_linear=100.0) + assert C > 0 + + def test_increases_with_snr(self): + tx = np.array([[0.0, 0.0, 0.0], [0.2, 0.0, 0.0]]) + rx = np.array([[0.0, 0.0, 10.0], [0.0, 0.2, 10.0]]) + H = compute_channel_matrix(tx, rx, 0.01) + C1 = spectral_efficiency(H, snr_linear=10.0) + C2 = spectral_efficiency(H, snr_linear=1000.0) + assert C2 > C1 + + +class TestDegreesOfFreedom: + def test_returns_int(self): + tx = np.array([[0.0, 0.0, 0.0], [0.2, 0.0, 0.0]]) + rx = np.array([[0.0, 0.0, 10.0], [0.0, 0.2, 10.0]]) + H = compute_channel_matrix(tx, rx, 0.01) + dof = degrees_of_freedom(H) + assert isinstance(dof, int) + assert 1 <= dof <= 2 diff --git a/tests/test_smoke.py b/tests/test_smoke.py new file mode 100644 index 0000000..c5f0f26 --- /dev/null +++ b/tests/test_smoke.py @@ -0,0 +1,53 @@ +""" +Smoke tests: each script runs end-to-end with --quick and exits without error. + +These do NOT check figure correctness — they verify the import chain and +code paths execute without crashing on tiny inputs. +""" +import subprocess +import sys +import os +import pytest + +PYTHON = sys.executable +REPO_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), '..')) +SCRIPTS_DIR = os.path.join(REPO_ROOT, 'scripts') + + +def run_script(name, extra_args=None): + script = os.path.join(SCRIPTS_DIR, name) + cmd = [PYTHON, script, '--quick'] + if extra_args: + cmd.extend(extra_args) + env = os.environ.copy() + env['PYTHONPATH'] = REPO_ROOT + os.pathsep + env.get('PYTHONPATH', '') + result = subprocess.run(cmd, capture_output=True, text=True, timeout=120, env=env) + if result.returncode != 0: + pytest.fail( + f"{name} failed (exit {result.returncode}):\n" + f"STDOUT:\n{result.stdout}\nSTDERR:\n{result.stderr}" + ) + + +class TestSmokeScripts: + def test_fig04_gain_vs_arrays(self): + run_script("fig04_gain_vs_arrays.py") + + def test_fig06_mimo_boundaries(self): + run_script("fig06_mimo_boundaries.py") + + def test_fig09_beampattern_sim(self): + run_script("fig09_beampattern_sim.py") + + def test_fig10_hardware_validation(self): + # Hardware data is not expected in CI — quick mode still runs theory+sim + run_script("fig10_hardware_validation.py") + + def test_fig11_2d_beampattern(self): + run_script("fig11_2d_beampattern.py") + + def test_fig12_mimo_dof(self): + run_script("fig12_mimo_dof.py") + + def test_fig13_throughput(self): + run_script("fig13_throughput.py") From 2a66c2f572e59d129cf7cda8bd13fb03bc5164f2 Mon Sep 17 00:00:00 2001 From: rohithreddyvennam Date: Wed, 27 May 2026 13:13:37 -0700 Subject: [PATCH 02/10] Add pytest-timeout to environment to fix CI smoke test --- environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/environment.yml b/environment.yml index 19948b7..4001f29 100644 --- a/environment.yml +++ b/environment.yml @@ -11,6 +11,7 @@ dependencies: - h5py>=3.8 - pyyaml>=6.0 - pytest>=7.0 + - pytest-timeout>=2.1 - pip - pip: - scienceplots>=2.0 # not on conda-forge From c5aa08e2983e7630362fab669c34ac56c9aa1bcf Mon Sep 17 00:00:00 2001 From: rohithreddyvennam Date: Wed, 27 May 2026 13:13:44 -0700 Subject: [PATCH 03/10] Fix Quick Start units and return value bugs in README --- Readme.md | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/Readme.md b/Readme.md index 71b5a35..db63691 100644 --- a/Readme.md +++ b/Readme.md @@ -35,14 +35,15 @@ r_min, r_max = mimo_region_bounds(d_tx=2000, d_rx=1.0, wavelength=lam) print(f"MIMO region: {r_min/1e3:.0f} km – {r_max/1e3:.0f} km") # Build the ArrayLink ground station (16 panels, center-dense layout) -gnd = build_ground_station( +# Lx/Ly are in km; element_spacing is also in km +gnd, _ = build_ground_station( mode="arraylink", subarray_shape=(32, 32), - element_spacing=lam / 2, + element_spacing=lam / 2 / 1e3, # λ/2 converted to km Nx=4, Ny=4, - Lx=1414.2, Ly=1000.0, + Lx=1.4142, Ly=1.0, # 1.414 km × 1.0 km aperture ) -print(f"Total antenna elements: {len(gnd)}") +print(f"Total antenna elements: {len(gnd)}") # → 16384 ``` ## Reproducing Paper Figures From 30ef76ba2cbd3c226b9915b694dea6d454e16ac9 Mon Sep 17 00:00:00 2001 From: rohithreddyvennam Date: Wed, 27 May 2026 13:49:37 -0700 Subject: [PATCH 04/10] Add hardware metrics dataset and fix fig10 case remapping --- .gitignore | 2 + Readme.md | 9 ++-- data/README.md | 75 +++++++++++++++------------ data/hardware_metrics.pkl | Bin 0 -> 13308 bytes scripts/fig10_hardware_validation.py | 25 ++++++++- 5 files changed, 73 insertions(+), 38 deletions(-) create mode 100644 data/hardware_metrics.pkl diff --git a/.gitignore b/.gitignore index 7c6c8d9..256f18c 100644 --- a/.gitignore +++ b/.gitignore @@ -13,6 +13,8 @@ results/ # Temporary files *.pkl +# Exception: hardware dataset file is tracked intentionally +!data/hardware_metrics.pkl /tmp/ # IDE diff --git a/Readme.md b/Readme.md index db63691..a66d96c 100644 --- a/Readme.md +++ b/Readme.md @@ -113,9 +113,10 @@ environment.yml # Mamba/conda environment If you use this simulator, please cite: ```bibtex -@inproceedings{arraylink2026, - title = {ArrayLink: Distributed Ground Station for Near-Field LoS MIMO Satellite Communications}, - booktitle = {IEEE INFOCOM}, - year = {2026}, +@article{vennam2025satellites, + title = {Satellites are closer than you think: A near field MIMO approach for Ground stations}, + author = {Vennam, Rohith Reddy and Wilson, Luke and Jain, Ish Kumar and Bharadia, Dinesh}, + journal = {arXiv preprint arXiv:2508.09374}, + year = {2025}, } ``` diff --git a/data/README.md b/data/README.md index bc7bddb..b693fd2 100644 --- a/data/README.md +++ b/data/README.md @@ -1,45 +1,56 @@ # Hardware Data -The hardware experiment data for Fig. 10 will be released as part of the -ArrayLink dataset. +Hardware experiment data for Fig. 10 (singular-value ratio vs. distance). -## Expected file +## Included file -Place `hardware_metrics.pkl` in this directory (i.e. `data/hardware_metrics.pkl`). +`hardware_metrics.pkl` is included in this directory. +It contains pre-processed aggregate statistics extracted from the raw hardware +measurements and is all that is needed to reproduce Fig. 10. -## Format - -The pickle file should contain a Python `dict` with the following structure: +## What the file contains ```python -{ - 'case1': { - 'distances': [2.5, 5.0, 10.0, ...], # metres - 'mean_sing_ratio': [0.82, 0.75, 0.61, ...], # mean sigma_2/sigma_1 - 'std_sing_ratio': [0.03, 0.04, 0.05, ...], # std dev across packets - }, - 'case2': { ... }, - 'case3': { ... }, - 'case4': { ... }, -} +import pickle +with open("data/hardware_metrics.pkl", "rb") as f: + data = pickle.load(f) + +# data is a dict keyed by case label (satellite-ground-station convention): +# { +# 'case1': { # Drx=50 cm, Dtx=50 cm +# 'distances': np.array([2.5, 5.0, 10.0, ...]), # metres +# 'mean_sing_ratio': np.array([...]), # mean σ₂/σ₁ across packets +# 'std_sing_ratio': np.array([...]), # std σ₂/σ₁ across packets +# ... (additional keys not used by fig10) +# }, +# 'case2': { ... }, # Drx=20 cm, Dtx=50 cm +# 'case3': { ... }, # Drx=20 cm, Dtx=20 cm +# 'case4': { ... }, # Drx=50 cm, Dtx=20 cm +# } ``` -where case1–case4 correspond to the four transmit/receive aperture configurations -in Table I of the paper: +`fig10_hardware_validation.py` remaps these keys to ArrayLink's case convention +automatically on load — no manual step needed. -| Case | d_tx (cm) | d_rx (cm) | -|------|-----------|-----------| -| 1 | 20 | 20 | -| 2 | 50 | 20 | -| 3 | 20 | 50 | -| 4 | 50 | 50 | +## Running Fig. 10 -## Running Fig 10 without hardware data +```bash +python scripts/fig10_hardware_validation.py --save-dir paper_figures +``` -`fig10_hardware_validation.py` will run with theory and simulation curves only -if the pkl file is absent. It prints a clear message and continues: +## Raw channel data -``` -Hardware data not found at: data/hardware_metrics.pkl -Proceeding with theory + simulation curves only. -``` +The full raw dataset (~31 GB of `.mat` channel matrix files) will be released +on Zenodo. A download link and DOI will be added here before camera-ready. + +Each raw file contains: +- Complex MIMO channel matrices `H[64 freq bins, 2×2, N packets]` +- Measured at 27 GHz with a 2×2 MIMO testbed +- Cases correspond to four (d_tx, d_rx) aperture configurations from Table I + +| Case | d_tx | d_rx | +|------|------|------| +| case1 | 50 cm | 50 cm | +| case2 | 50 cm | 20 cm | +| case3 | 20 cm | 20 cm | +| case4 | 20 cm | 50 cm | diff --git a/data/hardware_metrics.pkl b/data/hardware_metrics.pkl new file mode 100644 index 0000000000000000000000000000000000000000..449b4e92084e985a3fd0576808beeef3fafb3587 GIT binary patch literal 13308 zcmd5?c|29!*H0lCQi(^VN-3F1C1cU7q)F1ijmzcXk|7mkD5P#Bu1Y8gWeRDaLZPUn zIf)7tQF@Xh^yEEtZ_ao7z3+eT=l#6?==5FdyVl-o?Y+-g>zo@b7~3&@$nU?MTvd*c zBg2QO$Nlx+5l*f?z6`b_(}&BEVf%S_ZP&7Q^z>$GdHA{ex-z`I8QZxWaeHs3qbJ+P z*W1t0mm9<72>oVqIm6jbJW34La!450LSQpDg1eHdYAIkT5yoW==Lk9ZZuesTG8yV@ zz>VC?(+=YZ7;~+xtd{gsf4~2+9OBE}9BIMg%~5a;?>d%$oy@;#@vr*)>umnjlz*Ma zzb@uqt$5eq6@0v!e_c0+OAW`7^k6dB_TCI%S5IAg4~7eOBUhCr$PxB&WVkcDsY$a& zQPb9^r_CDspOdx_?AP@$qw1Ildkv*|kjWY8nwusr9uwE$NzjpSO17FL~W6F8c)l_)zuqV-QR&Z zJ`crQl&1zDlzz%ac;O?VW#sHctw#s3oQXYaQQP;y7U5K#(Fh$oHX}5ST8z-XekMZq z8MjS;XAl=`P+af!iefIhWCRe}ndu^wHjzOnrl5#WQu`@p*(ZQluq&;Tk-T$HHIh}0 zg?UruorifsLy)ZFwicnqh2}pCHX4tk7mK<@{XUStm+b0^OU1B(M9g*O0JL=a7qxmx z0VSYi_*5XAk)?^YjZ+N~O1@P@=u;Gop1;J+xelIhEggp1h&pZ5HjEbMxngog!5JGA z*{$$Xm=BL_Nv1G@>|{BF5~czOyL zBPL}bvT>)&hNIKl2OLtOp5!R@9r2eal8{k)XB$?JqOr}-+m9@*LD3Y|&M`X|KSojD z5tf;MaUF_U1K)oUcW6Y>m9uVp3{&5rh+C1n#cAIs6rCEADjlIru-=7+Ryf^hMOKJw z?Zh9w?@-ifoHkX;htFbiM*KRKgVHd&Zj>*LakWZmW~xSD$XEFrp#F|^-&Cp<{wZDe_5Tx6XAByqql=g;=Y6#rBBeTpl$vQ zBhMQsiqg2@FQ#z;MSdyvwaZu!Q6%ZDS@b;S0g4R0FT0Kxy^W$MwG+ktu>`AB!q#T< znsQ`q(0ww!DY^nhm8lW_r<$&#=z@r?^7)sAD7rRg&iOv)TPV8veCZ1#w<|mmYzIg~ z&=T2hE9;|Z@RUp^j>rN;Sp}cysXD$b>DwM z*QAuMh+XoHGOq%8hJOmDU&n%17Sem7c=+ewz>;S;f~9__@D6?Aq~G#?eWCXhk!an+ zwWeF&x1ga*-3^&3AsSOl$`A(R_?pV-7U{%9YI$7zaiSVU0RPwRaUptVsue?e| zU5w9S`tkZk-CunpYZyzIHR4y##}dPiPwM};hhxe9>e+smhz2)R2Ifw+VFE$B5!%EP zHNQ(49XmxnFoI4P1w-6@(Dv3pL`SEP`4$Cwm6d2&BviD?Yk3^Jb!@!u0m+RUYu28jUi9VbTh$N1lcVZ^4SNE zFd_G=FCpZn>mY0tE<)RcaT`#3W|<()_A$9GTlI{PoOQ1bp=~+QTQ2(h9S*e8MT$pe zx;YQTLOnKnjr4jx-fJ943MK03yc(u@oX#;N>O3y``qKOt$l86#xH?&ZU|A&VjIe7V zx}PY2m*^;E!f=0;p@QLH`f-7CyJ@7&D{dgszhOCyVmNoUI|e- z^=|7o;r?vEsaBk`RIFId@1&S+hcJFkU+=t!vxpKTq(el!KwWZRA*lG*Y-yjxSQUXz z<7n}^Z;<u?(OOE;X<8$*UY4G7qBR{vy`QlnX@f6~4&V^k+sCkHvH=!1xt!{_t}8mnb3i z*A00;CGsvh&9@Bj^=1=Q<%xjXSw>DoRau`{ySX!&=###Ob1WaY=Ay4g!L9?sY@%|0 z7%=*^i2+ex2ObnZsO3oPl51M7IvtMp(3f?8uF9D}q9)GSwEq0Gkq?n29?i&j-pmh( z>DGwxZ&+7TBVEEl`bb3)t;WYh<83T!y4+NDE3%yJd_Jh$BMc+V4rs(^cp>Zbs_-9c zxubFY$6m}8cCO1qmf3fQ(MhWO1>%*n>-Aa>ZTT68tS%+Z(M&GE%F8c0I*&`&fYg?d z>g^j-WkiU^O{@))ofgdV#dH%L`scn((N*IP+RT$64(L#^IycsP;+F99URPT1FrrUB zSoCz8&<*0nL)~l7i#szaaLzlX1PhnQ5}REy>QP+uIAX-;`I=|Kue?V?PQt6K=*}h- z)veqc{>XqB@k>~v{N%t|WL3_&FOW5oIRAIAXvqBdtA@v7x=qLUJ6mH<%j**-du5Zc z*1z3n&}r@r#cK=?y11qi=Wpr0`1xkyoGq<9zhR=%7Bnnz;4ouv5X)#Co+I7&mB^Ce z;by&ly8~!g`bK_hx6(-zdAMbSPE-j&QQqT8>}Ee=UPofsjnKJQZ z$@d!vwGn~DM~>^G&Y|a|s}e|V+A5DwO)?b85}y|Ev`px5az_y@&CPYh*Sv>b+#B%$MGJ8vFLp{ zkGXo5EQWdB*Qbo|laL}#KTup$?y|BJCCZu6p6j*qb`q6D+7zpg=0-#d@*pTQeP$t1^dx_r#@%B> zh%A);J^z<9A$mFVi*uVNQI(8!Qdz;3Aj(-Yrtaw!`PPB8qiPC++vZscUi7h{M0Y3e zFuC}ZC`Qvgt-iZBHKCz|=)3Rbg?!ea4(G+sbJejm2}qY)9aX!i?i7m5qb4>jSV%uj zsiBXBh4&3@NkrDCb<+$q_S27GigmrwZbrmCesRgGUFAW>nxWe3sE6Uu$qsD;d|bozx3CF;(W%`;LXzV(>zYD-VCy}t|*Gj&KRx~ zPVCFq8rOp)8N?(;bQb4jxTf+92X#0v=H;pvcz-~rR}F8P{?O((1iU@F(IIT_^oZNY zS|R&1#9i+ristO_EFFH4XjCd>1jfZ168Dm)e%_q6aH0i?k8a<6Y!q>ODQ#_<>|si5 zFZVs1J(cf>EpgxS;SQNL{>wq$R3ptax!U9>p0!QQYr6LES%W$p&H6|YA3vf68Tqwq zt&%@oiczc4w)wlR#6e<#hGhjdTzOf5@e35~8F${DxRBFE%^Pc!Plq%8f9i0=sSi5THyno{zyJQ1euwuz{QeHU-MBG{pU@>W)3H;J;z&h0A_}Ih z@E*tv1>5cDwvvL|jr~x|*hAQs*${Q~9g~?zR^G3Pu=iR&K2(c3`7agmV~wBalm6Lk|Re>M{1*iep{8Luxj{+LobcdobIf1g25#dJ(l=BKABjwj+US(lFSf1H?zP;ninQ zd(edVO;c+3A=C!xu<=-JI&=WFqupu|_MCi;FuAo0Vbm+4l~D**;_nvn3rqKxkys9| zV-c{(#t*X*Dogx6a7Hctit^X)9Yb^plBzjqD_$>!X`f>~4-;!3uZ!d-TVJF0Vd1^v z--|NnlS=xgVhp)vT@GfiJc?MCkO^0j{5drS;l`w?IO0k-BJXy^@ZAgY3| zc84r}EG9(l;op=-r4j{fNxjLw!zbV3h=%{TC{}gQ$)+apnQiUp=hlI&*evVv%+N*@ zg*i>S#WWy9=9<|#wR`zB&7gL|Hg@0B(A(+g!64yl;6FOV^S;a{tZGdP{Zyi6Q}!g~ z=KEWOs9~6c5o;pxuJZZu<7+hKqS1Z(gyTASJCC6#Xv(TdV(?mua({_^;E+ngS!7AD zQwyjz8$~IzU6Nu;312D3OXntApX9Ly^$yEtO&RKQlPJivos1GE2@$VvSG#ObtKNDS z$B|XOe!@77DBE4*m-(-%zlyBJZ2fJQo)PuFY>2_8JM|^Vs=hE(Bh8nn_eIO)+U%BQ zAxpeX`lMkSeKVy}lCo?`!?)x#WGUClojqBd$rBA~6qc2SZ@<0dG}4dC)D|C|M*mKP za{A`+Rk>N}?#Pl3x)gBkCIdxnmS5FQ#yF#BWbhkj6E;z{r=G8|TBl5xL6m#1;MCRT z-1TUvQJE5NC`iAKrC62ECo{8_*duGUQ15ZWxpdJ)u~fc{d$1~E8;>=pE6|s=X`fd? z)YY9Cnd@xVyg;Xla(aE6A2y-L{N&T*;R|a~#N4G{U}g6RMK!Z03vs^_8)r$FEN4bR zC9*bdKY8;-U@3}DO^(^&+D{y>sKz`fR`8d6}-gD8q28vn>m+&{Yy8;`HzDHb^J7mq+JXr<8DF<5(X{_FU=eizALM QEDKDx=5HLPpP$x$0WqP)c>n+a literal 0 HcmV?d00001 diff --git a/scripts/fig10_hardware_validation.py b/scripts/fig10_hardware_validation.py index 87e228f..97c0752 100644 --- a/scripts/fig10_hardware_validation.py +++ b/scripts/fig10_hardware_validation.py @@ -54,6 +54,23 @@ 'case4': {'d_tx_cm': 50, 'd_rx_cm': 50}, } +# The hardware_metrics.pkl (from satellite-ground-station) uses a different +# case numbering convention. This map translates pkl keys → ArrayLink keys. +# +# pkl key Drx Dtx ArrayLink key +# case1 50cm 50cm → case4 +# case2 20cm 50cm → case2 (same apertures, same label — coincidence) +# case3 20cm 20cm → case1 +# case4 50cm 20cm → case3 +# case0 40cm 40cm → (not used in ArrayLink) +# +PKL_TO_ARRAYLINK = { + 'case1': 'case4', + 'case2': 'case2', + 'case3': 'case1', + 'case4': 'case3', +} + def parse_args(): ap = argparse.ArgumentParser(description="Fig 10: hardware validation") @@ -96,8 +113,12 @@ def main(): if os.path.exists(data_path): with open(data_path, 'rb') as f: - hw_data = pickle.load(f) - print(f"Loaded hardware data from {data_path}") + raw = pickle.load(f) + # Remap pkl case keys to ArrayLink's case convention (apertures differ) + hw_data = {PKL_TO_ARRAYLINK[k]: v for k, v in raw.items() + if k in PKL_TO_ARRAYLINK} + print(f"Loaded hardware data from {data_path} " + f"({len(hw_data)} cases remapped)") else: if args.quick: print(f"[quick mode] Hardware data not found at {data_path}. " From cb755f7cd3bf894d992882640a5789d48a8397bf Mon Sep 17 00:00:00 2001 From: rohithreddyvennam Date: Wed, 27 May 2026 13:53:22 -0700 Subject: [PATCH 05/10] Document hardware_metrics.pkl provenance and add dataset link placeholder --- data/README.md | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/data/README.md b/data/README.md index b693fd2..ade5a47 100644 --- a/data/README.md +++ b/data/README.md @@ -38,12 +38,31 @@ automatically on load — no manual step needed. python scripts/fig10_hardware_validation.py --save-dir paper_figures ``` +## How `hardware_metrics.pkl` was generated + +`hardware_metrics.pkl` was produced from the raw `.mat` channel files using: + +``` +python_simulator/targeted_scripts/hardware_experiments_analysis.py +``` + +At a high level, the script: +1. Reads a logbook CSV that maps each `.mat` file to its distance, case, and aperture config +2. Loads each `.mat` file (complex MIMO channel matrices, shape `[64 freq bins, 2×2, N packets]`) +3. At a fixed frequency bin, computes σ₂/σ₁ (singular-value ratio) for every captured packet +4. Aggregates per distance point → `mean_sing_ratio`, `std_sing_ratio` +5. Saves the result as `hardware_metrics.pkl` + +> ⚠️ This script is part of the raw-data repository (link below) and requires +> the full `.mat` dataset to run. It is **not** needed to reproduce Fig. 10 — +> `hardware_metrics.pkl` is already included here. + ## Raw channel data -The full raw dataset (~31 GB of `.mat` channel matrix files) will be released -on Zenodo. A download link and DOI will be added here before camera-ready. +> 📌 **Dataset link coming soon** — the full raw dataset will be released on Zenodo +> before camera-ready. This section will be updated with the DOI and download instructions. -Each raw file contains: +Each raw `.mat` file contains: - Complex MIMO channel matrices `H[64 freq bins, 2×2, N packets]` - Measured at 27 GHz with a 2×2 MIMO testbed - Cases correspond to four (d_tx, d_rx) aperture configurations from Table I From 4f63ae0e62bb048f963582c4d5af87e39e726a89 Mon Sep 17 00:00:00 2001 From: rohithreddyvennam Date: Wed, 27 May 2026 15:47:59 -0700 Subject: [PATCH 06/10] Add Apache 2.0 license, fig02, YAML config loading, and bug fixes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit License - Add LICENSE (Apache 2.0, copyright 2025 UCSD WCSNG) - Add SPDX-License-Identifier headers to all 18 Python source files - Add license badge and License section to README New figure - scripts/fig02_parabolic_gain.py: Airy-disk beam pattern for 1.85 m and 1.47 m dishes (η=0.6) at 28 GHz; outputs full ±90° and zoomed ±1° PDFs - arraylink/beamforming.py: add parabolic_beam_pattern_dbi() using scipy j1; stable near θ=0 via masked denominator YAML config loading (fig09 / fig11 / fig12) - Load frequency, aperture, subarray shape, and element gain from configs/arraylink_1km.yaml via --config CLI arg; no more scattered module-level magic numbers - Fix PyYAML float parsing: 28.0e9 → 28.0e+9 in both config files Bug fixes - channel.degrees_of_freedom: use σ_k/σ_1 ≥ τ (paper Eq. 7), not L2-normalised absolute threshold; fixes undercounting of spatial streams - fig12_mimo_dof.py: same DoF criterion fix applied inline - fig06 quick mode: √2e3 → √2×1e3 (was ~45 m instead of ~1414 m) - fig13 config labels: dynamic f-strings instead of hardcoded "52.6 dBi" - fig13 --quick: print completion message instead of silent exit - arraylink/__init__.py: export parabolic_beam_pattern_dbi and load_config Tests - test_channel.py: add test_uses_sigma1_ratio_not_l2_norm (guards DoF fix) - test_channel.py: add test_rank1_channel_dof_is_1 --- LICENSE | 181 +++++++++++++++++++++++++++ Readme.md | 11 ++ arraylink/__init__.py | 5 + arraylink/array_geometry.py | 3 + arraylink/beamforming.py | 43 ++++++- arraylink/channel.py | 17 ++- arraylink/utils.py | 39 +++++- configs/arraylink_1km.yaml | 2 +- configs/upa_baseline.yaml | 2 +- data/README.md | 3 +- scripts/fig02_parabolic_gain.py | 146 +++++++++++++++++++++ scripts/fig04_gain_vs_arrays.py | 3 + scripts/fig06_mimo_boundaries.py | 5 +- scripts/fig09_beampattern_sim.py | 56 +++++---- scripts/fig10_hardware_validation.py | 3 + scripts/fig11_2d_beampattern.py | 69 +++++----- scripts/fig12_mimo_dof.py | 39 +++--- scripts/fig13_throughput.py | 12 +- tests/__init__.py | 3 + tests/test_array_geometry.py | 3 + tests/test_beamforming.py | 3 + tests/test_channel.py | 30 +++++ tests/test_smoke.py | 3 + 23 files changed, 599 insertions(+), 82 deletions(-) create mode 100644 LICENSE create mode 100644 scripts/fig02_parabolic_gain.py diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..cc0d5b0 --- /dev/null +++ b/LICENSE @@ -0,0 +1,181 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship made available under + the License, as indicated by a copyright notice that is included in + or attached to the work (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean, as submitted to the Licensor for inclusion + in the Work by the copyright owner or by an individual or Legal Entity + authorized to submit on behalf of the copyright owner. For the purposes + of this definition, "submitted" means any form of electronic, verbal, + or written communication sent to the Licensor or its representatives, + including but not limited to communication on electronic mailing lists, + source code control systems, and issue tracking systems that are managed + by, or on behalf of, the Licensor for the purpose of discussing and + improving the Work, but excluding communication that is conspicuously + marked or designated in writing by the copyright owner as "Not a + Contribution." + + "Contributor" shall mean Licensor and any Legal Entity on behalf of + whom a Contribution has been received by the Licensor and incorporated + within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contributions + with the Work to which such Contributions were submitted. If You + institute patent litigation against any entity (including a cross-claim + or counterclaim in a lawsuit) alleging that the Work or any + Contribution embodied within the Work constitutes direct or + contributory patent infringement, then any patent licenses granted + to You under this License for that Work shall terminate as of the + date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or Derivative + Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, You must include a readable copy of the + attribution notices contained within such NOTICE file, in + at least one of the following places: within a NOTICE text file + distributed as part of the Derivative Works; within the Source + form or documentation, if provided along with the Derivative + Works; or, within a display generated by the Derivative Works, + if and wherever such third-party notices normally appear. The + contents of the NOTICE file are for informational purposes only + and do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own license statement for Your modifications and + may provide additional grant of rights to use, copy, modify, merge, + publish, distribute, sublicense, and/or sell copies of the + Derivative Works, and to permit persons to whom the Derivative + Works is furnished to do so, subject to Your terms and conditions. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or reproducing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or exemplary damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or all other + commercial damages or losses), even if such Contributor has been + advised of the possibility of such damages. + + 9. Accepting Warranty or Liability. While redistributing the Work or + Derivative Works thereof, You may choose to offer, and charge a fee + for, acceptance of support, warranty, indemnity, or other liability + obligations and/or rights consistent with this License. However, in + accepting such obligations, You may offer only conditions that You + can solely fulfill, and not on behalf of any other Contributor. + + END OF TERMS AND CONDITIONS + + Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia + UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/Readme.md b/Readme.md index a66d96c..26f0c31 100644 --- a/Readme.md +++ b/Readme.md @@ -1,5 +1,7 @@ # ArrayLink Simulator +[![License](https://img.shields.io/badge/License-Apache%202.0-blue.svg)](LICENSE) + Near-field LoS MIMO simulator for the ArrayLink distributed ground station (INFOCOM 2026). ArrayLink uses 16 phased-array panels spread across a km-scale aperture to form a @@ -52,6 +54,7 @@ Run any figure script from the repo root: | Script | Figure | Description | |--------|--------|-------------| +| `scripts/fig02_parabolic_gain.py` | Fig. 2 | Parabolic dish beam pattern vs scan angle (1.85 m and 1.47 m dishes) | | `scripts/fig04_gain_vs_arrays.py` | Fig. 4 | Array gain vs number of panels | | `scripts/fig06_mimo_boundaries.py` | Fig. 6 | Theoretical MIMO region boundaries | | `scripts/fig09_beampattern_sim.py` | Fig. 9 | Simulation setup: UPA/ArrayLink positions and gain vs θ / distance | @@ -62,6 +65,7 @@ Run any figure script from the repo root: ```bash # Generate all figures (outputs to paper_figures/) +# Runtime: ~2-3 min total (fig11 ~30 s, fig12 ~1 min, others seconds each) for s in scripts/fig*.py; do python $s; done # Interactive 2D beam pattern @@ -108,6 +112,13 @@ environment.yml # Mamba/conda environment | Elements per panel | 32×32 = 1024 | Half-wavelength spacing | | Center-dense exponent γ | 3.0 | Power-law placement transform | +## License + +This project is licensed under the **Apache License 2.0** — see the [LICENSE](LICENSE) file for details. + +Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) + ## Citation If you use this simulator, please cite: diff --git a/arraylink/__init__.py b/arraylink/__init__.py index 38b6285..e12df4d 100644 --- a/arraylink/__init__.py +++ b/arraylink/__init__.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """ArrayLink near-field LoS MIMO simulator.""" from .channel import ( compute_channel_matrix, @@ -17,6 +20,7 @@ ) from .beamforming import ( parabolic_gain_dbi, + parabolic_beam_pattern_dbi, total_array_gain_dbi, marginal_gain_dbi, dc_weights, @@ -25,6 +29,7 @@ compute_beam_pattern_numpy, ) from .utils import ( + load_config, spherical2cartesian, cartesian2spherical, linear2db, diff --git a/arraylink/array_geometry.py b/arraylink/array_geometry.py index 50d91b5..41d563d 100644 --- a/arraylink/array_geometry.py +++ b/arraylink/array_geometry.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """ Antenna array placement functions. diff --git a/arraylink/beamforming.py b/arraylink/beamforming.py index 525c60f..af8eeb8 100644 --- a/arraylink/beamforming.py +++ b/arraylink/beamforming.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """ Beamforming weight generation and beam pattern computation. @@ -33,7 +36,7 @@ def parabolic_gain_dbi(D, frequency, efficiency=0.6): """ - Directional gain of a parabolic dish antenna (dBi). + On-axis (peak) gain of a circular parabolic dish antenna (dBi). G = eta * (pi*D/lambda)^2 @@ -55,6 +58,44 @@ def parabolic_gain_dbi(D, frequency, efficiency=0.6): return 10.0 * np.log10(gain_linear) +def parabolic_beam_pattern_dbi(D_m, frequency_hz, theta_deg, efficiency=0.6): + """ + Full beam pattern of a circular parabolic dish vs scan angle (dBi). + + Uses the Airy-disk (uniformly illuminated circular aperture) pattern: + + G(θ) = G0 · [2 J1(x) / x]², x = π·D·sin(θ) / λ + + where G0 = eta · (π·D/λ)² is the on-axis gain. + + Paper reference: Fig. 2 — parabolic dish gain pattern. + + Parameters + ---------- + D_m : dish diameter (metres) + frequency_hz: carrier frequency (Hz) + theta_deg : scan angles (degrees), scalar or array + efficiency : aperture efficiency eta (default 0.6) + + Returns + ------- + gain_dbi : ndarray, same shape as theta_deg, in dBi + (clipped at −100 dBi to avoid −∞ at nulls) + """ + from scipy.special import j1 + theta = np.deg2rad(np.asarray(theta_deg, dtype=float)) + lam = 3e8 / frequency_hz + k = 2 * np.pi / lam + a = D_m / 2.0 + G0 = efficiency * (np.pi * D_m / lam) ** 2 + x = k * a * np.sin(theta) + # Stable evaluation: replace zeros with 1 before dividing, restore limit = 1 + x_safe = np.where(np.abs(x) < 1e-12, 1.0, x) + pattern = np.where(np.abs(x) < 1e-12, 1.0, (2.0 * j1(x_safe) / x_safe) ** 2) + gain_linear = np.maximum(G0 * pattern, 1e-10) # floor avoids log(0) + return 10.0 * np.log10(gain_linear) + + # --------------------------------------------------------------------------- # Gain aggregation model # --------------------------------------------------------------------------- diff --git a/arraylink/channel.py b/arraylink/channel.py index 6dfc98b..27565db 100644 --- a/arraylink/channel.py +++ b/arraylink/channel.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """ Free-space LoS channel matrix and near-field MIMO analysis. @@ -69,12 +72,18 @@ def singular_value_ratio(H): def degrees_of_freedom(H, threshold=0.1): """ - Number of spatial streams whose normalised singular value exceeds `threshold`. + Number of spatial streams satisfying the paper's MIMO feasibility criterion. - threshold=0.1 matches the feasibility criterion in Eq. (7) of the paper. + Counts k where σ_k / σ_1 ≥ threshold (paper Eq. 7, default threshold=0.1). + + Note: uses the ratio σ_k/σ_1, NOT the L2-normalised absolute value. + The two are equivalent only for rank-1 channels; for multi-stream channels + the L2 norm suppresses all values and would under-count DoF. """ - s = compute_singular_values(H, normalise=True) - return int(np.sum(s >= threshold)) + s = compute_singular_values(H, normalise=False) + if s[0] == 0: + return 0 + return int(np.sum(s / s[0] >= threshold)) def spectral_efficiency(H, snr_linear): diff --git a/arraylink/utils.py b/arraylink/utils.py index d1fd6a7..2f8487f 100644 --- a/arraylink/utils.py +++ b/arraylink/utils.py @@ -1,11 +1,48 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """ -Coordinate conversion and decibel helpers. +Coordinate conversion, decibel helpers, and config loading. All coordinates in kilometres unless otherwise noted in the function signature. """ +import os import numpy as np +def load_config(path): + """ + Load a YAML configuration file. + + The path may be absolute or relative to the repository root. + The repository root is inferred as the parent directory of the + ``arraylink/`` package (i.e. two levels up from this file). + + Parameters + ---------- + path : str + Path to the YAML config file, e.g. ``'configs/arraylink_1km.yaml'``. + + Returns + ------- + cfg : dict + Parsed YAML contents. + + Examples + -------- + >>> from arraylink.utils import load_config + >>> cfg = load_config('configs/arraylink_1km.yaml') + >>> cfg['frequency_hz'] + 28000000000.0 # stored as 28.0e+9 in YAML → parsed as float by PyYAML + """ + import yaml + if not os.path.isabs(path): + repo_root = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + path = os.path.join(repo_root, path) + with open(path) as f: + return yaml.safe_load(f) + + def spherical2cartesian(r, theta, phi, deg=False): """ r : radius (km) diff --git a/configs/arraylink_1km.yaml b/configs/arraylink_1km.yaml index b7a8698..9036b81 100644 --- a/configs/arraylink_1km.yaml +++ b/configs/arraylink_1km.yaml @@ -2,7 +2,7 @@ # Matches the satellite-scale simulation in the paper (Sec. IV-C, Fig. 9b) # System parameters -frequency_hz: 28.0e9 # 28 GHz (Ka-band) +frequency_hz: 28.0e+9 # 28 GHz (Ka-band) # wavelength_km derived automatically as c / frequency # Ground station (ArrayLink) diff --git a/configs/upa_baseline.yaml b/configs/upa_baseline.yaml index 1e1c7e2..af30380 100644 --- a/configs/upa_baseline.yaml +++ b/configs/upa_baseline.yaml @@ -1,7 +1,7 @@ # UPA baseline: 128x128 monolithic uniform planar array, 28 GHz # Matches the monolithic comparison in the paper (Sec. IV-C, Fig. 9a) -frequency_hz: 28.0e9 +frequency_hz: 28.0e+9 ground_station: mode: uniform diff --git a/data/README.md b/data/README.md index ade5a47..1a195be 100644 --- a/data/README.md +++ b/data/README.md @@ -59,8 +59,7 @@ At a high level, the script: ## Raw channel data -> 📌 **Dataset link coming soon** — the full raw dataset will be released on Zenodo -> before camera-ready. This section will be updated with the DOI and download instructions. +> 📌 **Dataset link coming soon** — the full raw dataset will be released soon. This section will be updated with the DOI and download instructions. Each raw `.mat` file contains: - Complex MIMO channel matrices `H[64 freq bins, 2×2, N packets]` diff --git a/scripts/fig02_parabolic_gain.py b/scripts/fig02_parabolic_gain.py new file mode 100644 index 0000000..8345ac8 --- /dev/null +++ b/scripts/fig02_parabolic_gain.py @@ -0,0 +1,146 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) +""" +Fig 2 — Parabolic dish beam pattern vs scan angle. + +Reproduces Fig. 2 from the ArrayLink paper (INFOCOM 2026). + +Shows the far-field gain (dBi) of two circular parabolic dishes as a function +of scan angle θ ∈ [−90°, +90°] at 28 GHz: + - D = 1.85 m (η = 0.577) → ~52.3 dBi peak [baseline dish] + - D = 1.47 m (η = 0.50) → ~49.7 dBi peak [smaller reference] + +The pattern uses the Airy-disk model for a uniformly illuminated circular +aperture: G(θ) = G0 · [2 J1(x)/x]², x = π·D·sin(θ)/λ + +Two output files +---------------- + fig02_parabolic_gain.pdf — full ±90° scan range + fig02_parabolic_gain_zoomed.pdf — zoomed ±1°, 40–52.5 dB (main-lobe detail) + +Usage +----- + python scripts/fig02_parabolic_gain.py [--quick] [--save-dir paper_figures] + [--config configs/arraylink_1km.yaml] +""" +import argparse +import os +import sys +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +import numpy as np +import matplotlib.pyplot as plt + +try: + import scienceplots + plt.style.use(['science', 'no-latex']) +except ImportError: + pass + +from arraylink.beamforming import parabolic_beam_pattern_dbi +from arraylink.utils import load_config + + +# --------------------------------------------------------------------------- +# Fig-2-specific dish configurations (paper Sec. II-A / Table I) +# η=0.6 (60% aperture efficiency) used consistently across all figures. +# --------------------------------------------------------------------------- +DISHES = [ + {"label": "Diameter=1.85 m", "D_m": 1.85, "eta": 0.6, + "color": "steelblue", "ls": "-"}, + {"label": "Diameter=1.47 m", "D_m": 1.47, "eta": 0.6, + "color": "darkorange", "ls": "--"}, +] + + +def parse_args(): + ap = argparse.ArgumentParser( + description="Fig 2: parabolic dish beam pattern vs scan angle" + ) + ap.add_argument("--config", default="configs/arraylink_1km.yaml", + help="Path to YAML config (default: configs/arraylink_1km.yaml)") + ap.add_argument("--quick", action="store_true", + help="Print peak gains and exit, skip plot save") + ap.add_argument("--save-dir", default="paper_figures") + return ap.parse_args() + + +def main(): + args = parse_args() + os.makedirs(args.save_dir, exist_ok=True) + + cfg = load_config(args.config) + F_HZ = cfg['frequency_hz'] + + # Angular grid — high resolution to resolve narrow main lobe at 28 GHz + if args.quick: + theta_deg = np.linspace(-90, 90, 181) + else: + theta_deg = np.linspace(-90, 90, 18001) + + # Compute patterns + patterns = [] + for dish in DISHES: + gain = parabolic_beam_pattern_dbi( + dish["D_m"], F_HZ, theta_deg, efficiency=dish["eta"] + ) + patterns.append(gain) + peak = gain.max() + print(f" {dish['label']:20s} peak = {peak:.2f} dBi (η={dish['eta']})") + + if args.quick: + print("Quick mode: computations done, skipping plot save.") + return + + fig_kw = dict(fontsize=18, tick_labelsize=15) + + # ------------------------------------------------------------------ # + # Fig 2 — Full beam pattern (±90°) + # ------------------------------------------------------------------ # + fig, ax = plt.subplots(figsize=(10, 6)) + for dish, gain in zip(DISHES, patterns): + ax.plot(theta_deg, gain, + label=dish["label"], + color=dish["color"], + linestyle=dish["ls"], + linewidth=1.5) + + ax.set_xlabel("Scan Angle θ (degrees)", fontsize=fig_kw['fontsize']) + ax.set_ylabel("Gain (dB)", fontsize=fig_kw['fontsize']) + ax.set_xlim(-90, 90) + ax.set_ylim(-60, 60) + ax.tick_params(labelsize=fig_kw['tick_labelsize']) + ax.legend(fontsize=16, loc="upper right", framealpha=0.5) + ax.grid(True, which='major', alpha=0.6) + ax.minorticks_on() + ax.grid(True, which='minor', alpha=0.3, linestyle=':') + fig.tight_layout() + out = os.path.join(args.save_dir, "fig02_parabolic_gain.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out}") + + # ------------------------------------------------------------------ # + # Fig 2 (zoomed) — Main-lobe detail (±1°, 40–52.5 dB) + # ------------------------------------------------------------------ # + fig, ax = plt.subplots(figsize=(6, 6)) + for dish, gain in zip(DISHES, patterns): + ax.plot(theta_deg, gain, + label=dish["label"], + color=dish["color"], + linestyle=dish["ls"], + linewidth=3) + + ax.set_xlim(-1, 1) + ax.set_ylim(40, 52.5) + ax.tick_params(labelsize=fig_kw['tick_labelsize']) + ax.grid(True, which='major', alpha=0.6) + fig.tight_layout() + out = os.path.join(args.save_dir, "fig02_parabolic_gain_zoomed.pdf") + fig.savefig(out, bbox_inches='tight') + plt.close(fig) + print(f"Saved {out}") + + +if __name__ == "__main__": + main() diff --git a/scripts/fig04_gain_vs_arrays.py b/scripts/fig04_gain_vs_arrays.py index 44b7608..3991101 100644 --- a/scripts/fig04_gain_vs_arrays.py +++ b/scripts/fig04_gain_vs_arrays.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """ Fig 4a/4b — Total and marginal gain vs number of phased-array panels. diff --git a/scripts/fig06_mimo_boundaries.py b/scripts/fig06_mimo_boundaries.py index 396d565..8324a95 100644 --- a/scripts/fig06_mimo_boundaries.py +++ b/scripts/fig06_mimo_boundaries.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """ Fig 6 — Singular-value ratio vs distance showing MIMO feasibility boundaries. @@ -62,7 +65,7 @@ def main(): if args.quick: for d_tx, d_rx, label in [ (0.2, 0.2, "d_tx=d_rx=0.2m"), # 20 cm x 20 cm - (np.sqrt(2e3), np.sqrt(2), "satellite scale"), # sqrt(2) km x sqrt(2) m + (np.sqrt(2) * 1e3, np.sqrt(2), "satellite scale"), # sqrt(2) km x sqrt(2) m (metres, matching LAM_M) ]: r_min, r_max = mimo_region_bounds(d_tx, d_rx, LAM_M, tau=TAU) print(f"{label}: r_min={r_min:.2f}, r_max={r_max:.2f}") diff --git a/scripts/fig09_beampattern_sim.py b/scripts/fig09_beampattern_sim.py index 6617c8c..ec9b501 100644 --- a/scripts/fig09_beampattern_sim.py +++ b/scripts/fig09_beampattern_sim.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """ Fig 9 — Simulation setup and beam pattern results. @@ -35,19 +38,14 @@ from arraylink.array_geometry import upa_positions, place_subarrays from arraylink.beamforming import dc_weights, compute_beam_pattern_numpy +from arraylink.utils import load_config # -------------------------------------------------------------------------- -# Constants +# Fig-09-specific constants (not in YAML — visualization layout choices) # -------------------------------------------------------------------------- -F_HZ = 28e9 -LAM_KM = 3e8 / F_HZ / 1e3 # ~1.0714e-5 km -R_FOCAL = 500.0 # km -GRID_KM = 2.0 # panel centres within ±1 km × ±1 km -MIN_GAP = 0.1 # minimum separation between panel centres [km] -N_PANELS = 16 -SUB_SHAPE = (32, 32) # elements per panel -N_UPA = 128 # UPA dimension (128×128) -ELEM_GAIN = 6.0 # element gain assumption [dBi] +GRID_KM = 2.0 # panel centres drawn within ±1 km × ±1 km (display aperture) +MIN_GAP = 0.1 # minimum separation between panel centres [km] +N_UPA = 128 # UPA dimension (128×128) used only in this figure # -------------------------------------------------------------------------- @@ -94,23 +92,23 @@ def _build_arraylink(panel_xy_km, sub_shape, elem_spacing_km): # -------------------------------------------------------------------------- # Beam pattern helpers # -------------------------------------------------------------------------- -def _gain_vs_theta(ants_km, w, r_km, theta_deg): +def _gain_vs_theta(ants_km, w, r_km, theta_deg, lam_km): thetas = np.deg2rad(theta_deg) sat_pts = np.column_stack([ r_km * np.sin(thetas), np.zeros_like(thetas), r_km * np.cos(thetas), ]) - return compute_beam_pattern_numpy(sat_pts, ants_km, w, LAM_KM) # dB + return compute_beam_pattern_numpy(sat_pts, ants_km, w, lam_km) # dB -def _gain_vs_dist(ants_km, w, r_km_array): +def _gain_vs_dist(ants_km, w, r_km_array, lam_km): sat_pts = np.column_stack([ np.zeros(len(r_km_array)), np.zeros(len(r_km_array)), r_km_array, ]) - return compute_beam_pattern_numpy(sat_pts, ants_km, w, LAM_KM) # dB + return compute_beam_pattern_numpy(sat_pts, ants_km, w, lam_km) # dB # -------------------------------------------------------------------------- @@ -118,6 +116,8 @@ def _gain_vs_dist(ants_km, w, r_km_array): # -------------------------------------------------------------------------- def parse_args(): ap = argparse.ArgumentParser(description="Fig 9: simulation setup and beam pattern") + ap.add_argument("--config", default="configs/arraylink_1km.yaml", + help="Path to YAML config (default: configs/arraylink_1km.yaml)") ap.add_argument("--quick", action="store_true", help="Coarse grids for fast smoke test") ap.add_argument("--center-dense", action="store_true", @@ -130,6 +130,18 @@ def main(): args = parse_args() os.makedirs(args.save_dir, exist_ok=True) + # ------------------------------------------------------------------ # + # Load config — students can change these in configs/arraylink_1km.yaml + # and see the effect on all figures that use this config. + # ------------------------------------------------------------------ # + cfg = load_config(args.config) + F_HZ = cfg['frequency_hz'] + LAM_KM = 3e8 / F_HZ / 1e3 + R_FOCAL = cfg['target_satellite']['r_km'] + SUB_SHAPE = tuple(cfg['ground_station']['subarray_shape']) + N_PANELS = cfg['ground_station']['Nx'] * cfg['ground_station']['Ny'] + ELEM_GAIN = cfg['computation']['element_gain_dbi'] + if args.quick: theta_deg = np.linspace(-90, 90, 19) r_axis_km = np.linspace(50, 2000, 11) @@ -265,17 +277,17 @@ def main(): # ------------------------------------------------------------------ # # Compute beam patterns # ------------------------------------------------------------------ # - bp_upa_theta = _gain_vs_theta(upa_ants, w_upa, R_FOCAL, theta_deg) + ELEM_GAIN - bp_s0_theta = _gain_vs_theta(al_ants_s0, w_s0, R_FOCAL, theta_deg) + ELEM_GAIN - bp_s1_theta = _gain_vs_theta(al_ants_s1, w_s1, R_FOCAL, theta_deg) + ELEM_GAIN + bp_upa_theta = _gain_vs_theta(upa_ants, w_upa, R_FOCAL, theta_deg, LAM_KM) + ELEM_GAIN + bp_s0_theta = _gain_vs_theta(al_ants_s0, w_s0, R_FOCAL, theta_deg, LAM_KM) + ELEM_GAIN + bp_s1_theta = _gain_vs_theta(al_ants_s1, w_s1, R_FOCAL, theta_deg, LAM_KM) + ELEM_GAIN - bp_upa_dist = _gain_vs_dist(upa_ants, w_upa, r_axis_km) + ELEM_GAIN - bp_s0_dist = _gain_vs_dist(al_ants_s0, w_s0, r_axis_km) + ELEM_GAIN - bp_s1_dist = _gain_vs_dist(al_ants_s1, w_s1, r_axis_km) + ELEM_GAIN + bp_upa_dist = _gain_vs_dist(upa_ants, w_upa, r_axis_km, LAM_KM) + ELEM_GAIN + bp_s0_dist = _gain_vs_dist(al_ants_s0, w_s0, r_axis_km, LAM_KM) + ELEM_GAIN + bp_s1_dist = _gain_vs_dist(al_ants_s1, w_s1, r_axis_km, LAM_KM) + ELEM_GAIN if args.center_dense: - bp_cd_theta = _gain_vs_theta(al_ants_cd, w_cd, R_FOCAL, theta_deg) + ELEM_GAIN - bp_cd_dist = _gain_vs_dist(al_ants_cd, w_cd, r_axis_km) + ELEM_GAIN + bp_cd_theta = _gain_vs_theta(al_ants_cd, w_cd, R_FOCAL, theta_deg, LAM_KM) + ELEM_GAIN + bp_cd_dist = _gain_vs_dist(al_ants_cd, w_cd, r_axis_km, LAM_KM) + ELEM_GAIN boresight_upa = bp_upa_dist[np.argmin(np.abs(r_axis_km - R_FOCAL))] boresight_s0 = bp_s0_dist[np.argmin(np.abs(r_axis_km - R_FOCAL))] diff --git a/scripts/fig10_hardware_validation.py b/scripts/fig10_hardware_validation.py index 97c0752..fc4c9b1 100644 --- a/scripts/fig10_hardware_validation.py +++ b/scripts/fig10_hardware_validation.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """ Fig 10 — Hardware validation: singular-value ratio vs distance (4 aperture cases). diff --git a/scripts/fig11_2d_beampattern.py b/scripts/fig11_2d_beampattern.py index 92995d7..15b1669 100644 --- a/scripts/fig11_2d_beampattern.py +++ b/scripts/fig11_2d_beampattern.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """ Fig 11 — 2D beam patterns (transverse × range heat maps). @@ -64,6 +67,7 @@ from arraylink.beamforming import ( dc_weights, compute_beam_pattern_numpy, ) +from arraylink.utils import load_config def parse_args(): @@ -86,21 +90,14 @@ def parse_args(): help="RNG seed for random placement (default: 0)") ap.add_argument("--interactive", action="store_true", help="Open Plotly interactive figure (requires plotly)") + ap.add_argument("--config", default="configs/arraylink_1km.yaml", + help="Path to YAML config (default: configs/arraylink_1km.yaml)") ap.add_argument("--save-dir", default="paper_figures") return ap.parse_args() -ELEMENT_GAIN_DBI = 6.0 -TARGET_R_KM = 500.0 # DC focus range (km) — matches paper Fig. 11 - -# ArrayLink aperture dimensions (km) -AL_LX = np.sqrt(2) # 1.414 km along X -AL_LY = 1.0 # 1.000 km along Y -AL_NX = 4 # panels along X -AL_NY = 4 # panels along Y - - -def build_arraylink(placement, subarray_shape, element_spacing, seed): +def build_arraylink(placement, subarray_shape, element_spacing, seed, + nx, ny, lx, ly): """ Build the ArrayLink antenna array according to the chosen placement strategy. @@ -110,25 +107,27 @@ def build_arraylink(placement, subarray_shape, element_spacing, seed): subarray_shape : (M, N) elements per panel element_spacing: element pitch within a panel (km) seed : RNG seed (used for 'random' placement) + nx, ny : panels along X and Y (e.g. 4, 4) + lx, ly : aperture extents in km (e.g. sqrt(2), 1.0) Returns ------- - ants : (Nx*Ny*M*N, 3) all antenna positions - bases : (Nx*Ny, 3) panel base positions + ants : (nx*ny*M*N, 3) all antenna positions + bases : (nx*ny, 3) panel base positions """ - n_panels = AL_NX * AL_NY + n_panels = nx * ny if placement == "random": rng = np.random.default_rng(seed) - bx = rng.uniform(-AL_LX / 2.0, AL_LX / 2.0, n_panels) - by = rng.uniform(-AL_LY / 2.0, AL_LY / 2.0, n_panels) - bases = np.stack([bx, by, np.zeros(n_panels)], axis=1) # (16, 3) + bx = rng.uniform(-lx / 2.0, lx / 2.0, n_panels) + by = rng.uniform(-ly / 2.0, ly / 2.0, n_panels) + bases = np.stack([bx, by, np.zeros(n_panels)], axis=1) # (n_panels, 3) elif placement == "center-dense": _, bases = build_ground_station( mode='arraylink', subarray_shape=subarray_shape, element_spacing=element_spacing, - Nx=AL_NX, Ny=AL_NY, Lx=AL_LX, Ly=AL_LY, + Nx=nx, Ny=ny, Lx=lx, Ly=ly, gamma=3.0, seed=seed, ) @@ -136,7 +135,7 @@ def build_arraylink(placement, subarray_shape, element_spacing, seed): _, bases = build_ground_station( mode='uniform', subarray_shape=subarray_shape, element_spacing=element_spacing, - Nx=AL_NX, Ny=AL_NY, Lx=AL_LX, Ly=AL_LY, + Nx=nx, Ny=ny, Lx=lx, Ly=ly, ) else: @@ -147,7 +146,7 @@ def build_arraylink(placement, subarray_shape, element_spacing, seed): def compute_2d_pattern(gnd_ants, weights, lam_km, x_axis, y_axis, quick, - batch_size=2000): + elem_gain_dbi, batch_size=2000): """ Compute beam pattern on a Cartesian (X, Y) grid. @@ -160,7 +159,7 @@ def compute_2d_pattern(gnd_ants, weights, lam_km, x_axis, y_axis, quick, if quick: bp = compute_beam_pattern_numpy(flat_pts, gnd_ants, weights, lam_km) - bp += ELEMENT_GAIN_DBI + bp += elem_gain_dbi return bp.reshape(Nx, Ny).T # (Ny, Nx) # Batched numpy: process batch_size points at a time to keep RAM bounded @@ -175,15 +174,23 @@ def compute_2d_pattern(gnd_ants, weights, lam_km, x_axis, y_axis, quick, pct = 100 * start / N print(f" {pct:.0f}% ({start:,}/{N:,})", flush=True) gc.collect() - return (results + ELEMENT_GAIN_DBI).reshape(Nx, Ny).T # (Ny, Nx) + return (results + elem_gain_dbi).reshape(Nx, Ny).T # (Ny, Nx) def main(): args = parse_args() os.makedirs(args.save_dir, exist_ok=True) - F_HZ = 28e9 - LAM_KM = 3e8 / F_HZ / 1e3 + # Load config — students can change these in configs/arraylink_1km.yaml + cfg = load_config(args.config) + F_HZ = cfg['frequency_hz'] + LAM_KM = 3e8 / F_HZ / 1e3 + TARGET_R_KM = cfg['target_satellite']['r_km'] + ELEMENT_GAIN_DBI = cfg['computation']['element_gain_dbi'] + AL_NX = cfg['ground_station']['Nx'] + AL_NY = cfg['ground_station']['Ny'] + AL_LX = cfg['ground_station']['aperture_x_km'] + AL_LY = cfg['ground_station']['aperture_y_km'] if args.quick: # Minimal grid for smoke test — just verify code paths run @@ -206,12 +213,14 @@ def main(): # UPA 128×128 — compact monolithic array, λ/2 spacing upa_ants = upa_positions(128, 128, element_spacing=LAM_KM / 2.0) - # ArrayLink 4×4 panels × 32×32 elements — placement set by --placement flag + # ArrayLink panels × elements — placement set by --placement flag al_ants, al_bases = build_arraylink( placement=args.placement, - subarray_shape=(32, 32), + subarray_shape=tuple(cfg['ground_station']['subarray_shape']), element_spacing=LAM_KM / 2.0, seed=args.seed, + nx=AL_NX, ny=AL_NY, + lx=AL_LX, ly=AL_LY, ) w_upa = dc_weights(sat_loc, upa_ants, LAM_KM) @@ -224,10 +233,12 @@ def main(): f"[resolution={res_tag}]") print("Computing UPA beam pattern …", flush=True) - bp_upa = compute_2d_pattern(upa_ants, w_upa, LAM_KM, x_axis, y_axis, args.quick) + bp_upa = compute_2d_pattern(upa_ants, w_upa, LAM_KM, x_axis, y_axis, args.quick, + elem_gain_dbi=ELEMENT_GAIN_DBI) print("Computing ArrayLink beam pattern …", flush=True) - bp_al = compute_2d_pattern(al_ants, w_al, LAM_KM, x_axis, y_axis, args.quick) + bp_al = compute_2d_pattern(al_ants, w_al, LAM_KM, x_axis, y_axis, args.quick, + elem_gain_dbi=ELEMENT_GAIN_DBI) if args.quick: print("Quick mode: patterns computed, skipping plot save.") @@ -236,7 +247,7 @@ def main(): return peak = max(bp_upa.max(), bp_al.max()) - vmin = ELEMENT_GAIN_DBI # floor = single-element gain (matches paper colorbar) + vmin = ELEMENT_GAIN_DBI # floor = single-element gain (matches paper colorbar) vmax = peak print(f"Peak = {peak:.2f} dBi | colorbar [{vmin:.1f}, {vmax:.1f}] dBi") diff --git a/scripts/fig12_mimo_dof.py b/scripts/fig12_mimo_dof.py index f84b5aa..ad8d0c3 100644 --- a/scripts/fig12_mimo_dof.py +++ b/scripts/fig12_mimo_dof.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """ Fig 12 — LoS MIMO degrees of freedom vs distance for ArrayLink. @@ -28,11 +31,13 @@ from arraylink.channel import compute_channel_matrix, compute_singular_values from arraylink.array_geometry import build_ground_station, upa_positions -from arraylink.utils import spherical2cartesian +from arraylink.utils import spherical2cartesian, load_config def parse_args(): ap = argparse.ArgumentParser(description="Fig 12: MIMO DoF vs distance") + ap.add_argument("--config", default="configs/arraylink_1km.yaml", + help="Path to YAML config (default: configs/arraylink_1km.yaml)") ap.add_argument("--quick", action="store_true", help="Run on a coarse distance grid (fast, for smoke test)") ap.add_argument("--save-dir", default="paper_figures") @@ -43,22 +48,22 @@ def main(): args = parse_args() os.makedirs(args.save_dir, exist_ok=True) - # --- System parameters --- - F_HZ = 28e9 - C = 3e8 - LAM_KM = C / F_HZ / 1e3 # wavelength in km (~1.07e-5 km) - TAU = 0.1 # feasibility threshold + # --- Load config — students can change these in configs/arraylink_1km.yaml --- + cfg = load_config(args.config) + F_HZ = cfg['frequency_hz'] + LAM_KM = 3e8 / F_HZ / 1e3 # wavelength in km (~1.07e-5 km) + TAU = 0.1 # feasibility threshold - # Ground station: 16 x 32x32 panels over sqrt(2)km x 1km - GND_NX, GND_NY = 4, 4 - GND_LX, GND_LY = np.sqrt(2), np.sqrt(2) # km - PANEL_SHAPE = (32, 32) + # Ground station geometry from YAML + GND_NX = cfg['ground_station']['Nx'] + GND_NY = cfg['ground_station']['Ny'] + GND_LX = cfg['ground_station']['aperture_x_km'] # sqrt(2) km + GND_LY = cfg['ground_station']['aperture_y_km'] # 1.0 km + PANEL_SHAPE = tuple(cfg['ground_station']['subarray_shape']) - # Satellite: 2x2 element array, 1.414m x 1m aperture - # Assumption: 4-element array arranged as 2x2 UPA on a 1.414m x 1m grid - # in km: 1.414m = 1.414e-3 km, 1m = 1e-3 km - SAT_DX_KM = 1.414e-3 # 1.414 m in km - SAT_DY_KM = 1.0e-3 # 1.0 m in km + # Satellite: 2x2 element array (aperture from YAML, in metres → convert to km) + SAT_DX_KM = cfg['satellite']['aperture_x_m'] / 1e3 + SAT_DY_KM = cfg['satellite']['aperture_y_m'] / 1e3 if args.quick: dist_array_km = np.linspace(100, 3000, 5) @@ -102,7 +107,9 @@ def main(): ratio = float(s[k] / s[0]) if len(s) > k and s[0] > 0 else 0.0 sing_ratios[k].append(ratio) - dof = int(np.sum(s >= TAU)) + # Paper criterion: σ_k/σ_1 ≥ τ (Eq. 7). s[k]/s[0] = σ_k/σ_1 + # (L2-normalisation factors cancel in the ratio). + dof = int(np.sum(s / s[0] >= TAU)) if s[0] > 0 else 0 dof_list.append(dof) if args.quick: diff --git a/scripts/fig13_throughput.py b/scripts/fig13_throughput.py index 83c625d..606ce02 100644 --- a/scripts/fig13_throughput.py +++ b/scripts/fig13_throughput.py @@ -1,10 +1,13 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """ Fig 13 — Per-link and aggregate throughput comparison. Reproduces Fig. 13 from the ArrayLink paper (INFOCOM 2026). Compares: - - 1 stream @ 52.6 dBi (1.85 m dish baseline) + - 1 stream @ ~52.5 dBi (1.85 m dish baseline, η=0.6) - 2 streams @ 48.14 dBi (ArrayLink) - 4 streams @ 48.14 dBi (ArrayLink) @@ -61,9 +64,9 @@ def main(): SNR0_SWEEP = np.arange(5, 26, 2.5) configs = [ - {"name": "1 stream @ 52.6 dBi (dish)", "N": 1, "gain_dbi": G_DISH_185}, - {"name": "2 streams @ 48.14 dBi (AL)", "N": 2, "gain_dbi": G_ARRAYLINK}, - {"name": "4 streams @ 48.14 dBi (AL)", "N": 4, "gain_dbi": G_ARRAYLINK}, + {"name": f"1 stream @ {G_DISH_185:.1f} dBi (dish)", "N": 1, "gain_dbi": G_DISH_185}, + {"name": f"2 streams @ {G_ARRAYLINK:.2f} dBi (AL)", "N": 2, "gain_dbi": G_ARRAYLINK}, + {"name": f"4 streams @ {G_ARRAYLINK:.2f} dBi (AL)", "N": 4, "gain_dbi": G_ARRAYLINK}, ] markers = ['o', 's', '^'] lines = ['--', '-', '-.'] @@ -86,6 +89,7 @@ def main(): if args.quick: for cfg in configs: print(f"{cfg['name']}: per-link={perlink_gbps[cfg['name']]}, agg={agg_gbps[cfg['name']]}") + print("Quick mode: computations done, skipping plot save.") return fig_params = dict(fontsize=20, tick_labelsize=18) diff --git a/tests/__init__.py b/tests/__init__.py index e69de29..1922f34 100644 --- a/tests/__init__.py +++ b/tests/__init__.py @@ -0,0 +1,3 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) diff --git a/tests/test_array_geometry.py b/tests/test_array_geometry.py index 1d106e9..6c495e5 100644 --- a/tests/test_array_geometry.py +++ b/tests/test_array_geometry.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """Unit tests for arraylink.array_geometry.""" import numpy as np import pytest diff --git a/tests/test_beamforming.py b/tests/test_beamforming.py index 355a55f..44f0639 100644 --- a/tests/test_beamforming.py +++ b/tests/test_beamforming.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """Unit tests for arraylink.beamforming.""" import numpy as np import pytest diff --git a/tests/test_channel.py b/tests/test_channel.py index 8c7732f..5576660 100644 --- a/tests/test_channel.py +++ b/tests/test_channel.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """Unit tests for arraylink.channel.""" import numpy as np import pytest @@ -126,3 +129,30 @@ def test_returns_int(self): dof = degrees_of_freedom(H) assert isinstance(dof, int) assert 1 <= dof <= 2 + + def test_uses_sigma1_ratio_not_l2_norm(self): + """ + Confirm the criterion is σ_k/σ_1 ≥ threshold (paper Eq. 7), NOT + the L2-normalised absolute value. + + Construct a synthetic H whose SVD is [1, 1, 0.12, 0.03]: + - paper criterion (σ_k/σ_1 ≥ 0.1): DoF = 3 (σ_1, σ_2, σ_3 pass) + - L2-norm criterion (s_k ≥ 0.1) : DoF = 2 (only σ_1, σ_2 pass + because s_2 = 0.12/‖s‖ ≈ 0.084 < 0.1 after L2 normalisation) + """ + # Build a 4×4 diagonal matrix with the desired singular values + s_target = np.array([1.0, 1.0, 0.12, 0.03]) + H_synthetic = np.diag(s_target.astype(complex)) + dof = degrees_of_freedom(H_synthetic, threshold=0.1) + assert dof == 3, ( + f"Expected DoF=3 (paper σ_k/σ_1 criterion), got {dof}. " + "Likely using L2-norm threshold which under-counts streams." + ) + + def test_rank1_channel_dof_is_1(self): + """A rank-1 channel should always have DoF=1.""" + tx = np.array([[0.0, 0.0, 0.0]]) + rx = np.array([[0.0, 0.0, 10.0], [0.1, 0.0, 10.0], + [0.0, 0.1, 10.0], [0.1, 0.1, 10.0]]) + H = compute_channel_matrix(tx, rx, 0.01) # 4×1 → rank 1 + assert degrees_of_freedom(H) == 1 diff --git a/tests/test_smoke.py b/tests/test_smoke.py index c5f0f26..dac76c6 100644 --- a/tests/test_smoke.py +++ b/tests/test_smoke.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia +# UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) """ Smoke tests: each script runs end-to-end with --quick and exits without error. From 5b5234cf9bbd9a8255008bf62cdc2d408135e29f Mon Sep 17 00:00:00 2001 From: rohithreddyvennam Date: Wed, 27 May 2026 15:58:28 -0700 Subject: [PATCH 07/10] Simplify README for clarity and conciseness MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Move optional plotly install out of Installation into a callout under Reproducing Figures, next to the fig11 interactive flag - Merge Hardware data subsection into a matching callout, same pattern - Trim Quick Start to a single clean code block (remove verbose comments) - Rename Key Assumptions → Key Parameters; trim from 8 to 5 rows, drop Notes column and internal params already covered by YAML configs - Condense License section to one line + copyright --- Readme.md | 76 +++++++++++++++++++++---------------------------------- 1 file changed, 29 insertions(+), 47 deletions(-) diff --git a/Readme.md b/Readme.md index 26f0c31..c382832 100644 --- a/Readme.md +++ b/Readme.md @@ -4,9 +4,7 @@ Near-field LoS MIMO simulator for the ArrayLink distributed ground station (INFOCOM 2026). -ArrayLink uses 16 phased-array panels spread across a km-scale aperture to form a -distributed ground station that exploits spatial multiplexing in the radiative near-field -of LEO/MEO satellites. +ArrayLink uses 16 phased-array panels spread across a km-scale aperture to exploit spatial multiplexing in the radiative near-field of LEO/MEO satellites. ## Installation @@ -19,32 +17,21 @@ mamba env create -f environment.yml mamba activate arraylink ``` -For the interactive 2D beam pattern (Fig. 11 `--interactive` flag): -```bash -mamba install -c conda-forge plotly -``` - ## Quick Start ```python from arraylink.channel import compute_channel_matrix, mimo_region_bounds from arraylink.array_geometry import build_ground_station -from arraylink.beamforming import total_array_gain_dbi, parabolic_gain_dbi # MIMO feasibility range for the ArrayLink aperture at 28 GHz lam = 3e8 / 28e9 r_min, r_max = mimo_region_bounds(d_tx=2000, d_rx=1.0, wavelength=lam) -print(f"MIMO region: {r_min/1e3:.0f} km – {r_max/1e3:.0f} km") - -# Build the ArrayLink ground station (16 panels, center-dense layout) -# Lx/Ly are in km; element_spacing is also in km -gnd, _ = build_ground_station( - mode="arraylink", - subarray_shape=(32, 32), - element_spacing=lam / 2 / 1e3, # λ/2 converted to km - Nx=4, Ny=4, - Lx=1.4142, Ly=1.0, # 1.414 km × 1.0 km aperture -) +print(f"MIMO region: {r_min/1e3:.0f} – {r_max/1e3:.0f} km") + +# Build the ArrayLink ground station (16 panels, center-dense, √2 km × 1 km aperture) +gnd, _ = build_ground_station(mode="arraylink", subarray_shape=(32, 32), + element_spacing=lam/2/1e3, Nx=4, Ny=4, + Lx=1.4142, Ly=1.0) print(f"Total antenna elements: {len(gnd)}") # → 16384 ``` @@ -54,29 +41,28 @@ Run any figure script from the repo root: | Script | Figure | Description | |--------|--------|-------------| -| `scripts/fig02_parabolic_gain.py` | Fig. 2 | Parabolic dish beam pattern vs scan angle (1.85 m and 1.47 m dishes) | +| `scripts/fig02_parabolic_gain.py` | Fig. 2 | Parabolic dish beam pattern vs scan angle | | `scripts/fig04_gain_vs_arrays.py` | Fig. 4 | Array gain vs number of panels | | `scripts/fig06_mimo_boundaries.py` | Fig. 6 | Theoretical MIMO region boundaries | -| `scripts/fig09_beampattern_sim.py` | Fig. 9 | Simulation setup: UPA/ArrayLink positions and gain vs θ / distance | +| `scripts/fig09_beampattern_sim.py` | Fig. 9 | UPA/ArrayLink positions and beam pattern | | `scripts/fig10_hardware_validation.py` | Fig. 10 | Hardware experiment validation | | `scripts/fig11_2d_beampattern.py` | Fig. 11 | 2D beam pattern heatmap | | `scripts/fig12_mimo_dof.py` | Fig. 12 | MIMO degrees of freedom vs distance | | `scripts/fig13_throughput.py` | Fig. 13 | Throughput comparison | ```bash -# Generate all figures (outputs to paper_figures/) -# Runtime: ~2-3 min total (fig11 ~30 s, fig12 ~1 min, others seconds each) +# Generate all figures (outputs to paper_figures/; ~2–3 min total) for s in scripts/fig*.py; do python $s; done - -# Interactive 2D beam pattern -python scripts/fig11_2d_beampattern.py --interactive ``` -### Hardware data (Fig. 10) +> **Fig. 11 interactive mode** — install `plotly` first, then run with `--interactive`: +> ```bash +> mamba install -c conda-forge plotly +> python scripts/fig11_2d_beampattern.py --interactive +> ``` -Fig. 10 shows theory and simulation curves even without hardware data. To add the -measured curves, place `hardware_metrics.pkl` at `data/hardware_metrics.pkl`. See -[data/README.md](data/README.md) for the expected format. +> **Fig. 10 hardware curves** — place `hardware_metrics.pkl` at `data/hardware_metrics.pkl` +> (theory and simulation curves render without it; see [data/README.md](data/README.md)). ## Running Tests @@ -91,30 +77,26 @@ arraylink/ channel.py # Channel matrix, singular values, MIMO bounds, spectral efficiency array_geometry.py # UPA, center-dense, and uniform array placement beamforming.py # DC weights, beam pattern, array/dish gain formulas - utils.py # Coordinate transforms and helpers + utils.py # Coordinate transforms, decibel helpers, config loading scripts/ # One script per paper figure -configs/ # YAML parameter files for ArrayLink and UPA baseline +configs/ # YAML parameter files (arraylink_1km.yaml, upa_baseline.yaml) tests/ # Unit tests (42) + smoke tests (8) -data/ # Hardware experiment data (not included; see data/README.md) -environment.yml # Mamba/conda environment +data/ # Hardware experiment data (see data/README.md) ``` -## Key Assumptions +## Key Parameters -| Parameter | Value | Notes | -|-----------|-------|-------| -| Carrier frequency | 28 GHz | Ka-band (λ ≈ 10.7 mm) | -| Hardware experiment frequency | 27 GHz | Matches lab setup | -| Panel element gain | 6 dBi | Isotropic radiator baseline | -| Dish efficiency | 0.6 (60%) | Standard aperture efficiency | -| MIMO threshold τ | 0.1 | σ₂/σ₁ > τ for spatial multiplexing | -| Panels | 16 (4×4 layout) | Center-dense over √2 km × 1 km | -| Elements per panel | 32×32 = 1024 | Half-wavelength spacing | -| Center-dense exponent γ | 3.0 | Power-law placement transform | +| Parameter | Value | +|-----------|-------| +| Carrier frequency | 28 GHz | +| Panels | 16 (4×4, center-dense over √2 km × 1 km) | +| Elements per panel | 32×32 = 1024 (λ/2 spacing) | +| Dish aperture efficiency | 60% | +| MIMO threshold τ | 0.1 (σ₂/σ₁ > τ for spatial multiplexing) | ## License -This project is licensed under the **Apache License 2.0** — see the [LICENSE](LICENSE) file for details. +Apache License 2.0 — see [LICENSE](LICENSE). Copyright 2025 Rohith Reddy Vennam, Luke Wilson, Ish Kumar Jain, Dinesh Bharadia UC San Diego Wireless Communications Sensing and Networking Group (WCSNG) From 8c9e168a9d2ef2f78602092fdfa36951e99ae4d8 Mon Sep 17 00:00:00 2001 From: rohithreddyvennam Date: Wed, 27 May 2026 16:01:59 -0700 Subject: [PATCH 08/10] Improve README opening description for clarity and impact --- Readme.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Readme.md b/Readme.md index c382832..389806a 100644 --- a/Readme.md +++ b/Readme.md @@ -2,9 +2,9 @@ [![License](https://img.shields.io/badge/License-Apache%202.0-blue.svg)](LICENSE) -Near-field LoS MIMO simulator for the ArrayLink distributed ground station (INFOCOM 2026). +Conventional satellite ground stations rely on a single large parabolic dish. ArrayLink replaces it with 16 distributed phased-array panels spanning a km-scale aperture — placing LEO/MEO satellites in the radiative near-field where the line-of-sight channel supports spatial multiplexing across independent streams, delivering throughput that no single dish can match. -ArrayLink uses 16 phased-array panels spread across a km-scale aperture to exploit spatial multiplexing in the radiative near-field of LEO/MEO satellites. +This repository is the full simulator for the INFOCOM 2026 paper. ## Installation From 14c2cbe3740ed1258d6e923e152d972800e285e0 Mon Sep 17 00:00:00 2001 From: rohithreddyvennam Date: Wed, 27 May 2026 16:03:48 -0700 Subject: [PATCH 09/10] Fix README description to mention hardware experiment results --- Readme.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Readme.md b/Readme.md index 389806a..b142b1c 100644 --- a/Readme.md +++ b/Readme.md @@ -4,7 +4,7 @@ Conventional satellite ground stations rely on a single large parabolic dish. ArrayLink replaces it with 16 distributed phased-array panels spanning a km-scale aperture — placing LEO/MEO satellites in the radiative near-field where the line-of-sight channel supports spatial multiplexing across independent streams, delivering throughput that no single dish can match. -This repository is the full simulator for the INFOCOM 2026 paper. +This repository contains the simulator and hardware experiment results for the INFOCOM 2026 paper. ## Installation From d3b582611d93eeebdc7b8f529c01ed4a0e967f9b Mon Sep 17 00:00:00 2001 From: Rohith Reddy Vennam <48006337+rohithreddyvennam@users.noreply.github.com> Date: Wed, 27 May 2026 16:19:12 -0700 Subject: [PATCH 10/10] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- tests/test_smoke.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/test_smoke.py b/tests/test_smoke.py index dac76c6..ab02ee3 100644 --- a/tests/test_smoke.py +++ b/tests/test_smoke.py @@ -33,6 +33,9 @@ def run_script(name, extra_args=None): class TestSmokeScripts: + def test_fig02_parabolic_gain(self): + run_script("fig02_parabolic_gain.py") + def test_fig04_gain_vs_arrays(self): run_script("fig04_gain_vs_arrays.py")