Research Recipes#

Recipes for researchers exploring advanced array concepts.

Custom Optimization Objectives#

Problem: Optimize array thinning for minimum peak sidelobe.

import phased_array as pa
import numpy as np

# Full array to thin
full_geom = pa.create_rectangular_array(24, 24, dx=0.5, dy=0.5)
k = pa.wavelength_to_k(1.0)

def peak_sidelobe_objective(geom):
    """Objective function: minimize peak sidelobe level."""
    weights = np.ones(geom.n_elements)
    theta_deg, E_plane, _ = pa.compute_pattern_cuts(
        geom.x, geom.y, weights, k, n_points=721
    )

    # Find main beam region
    peak_idx = np.argmax(E_plane)
    half_bw = 10  # samples to exclude

    # Peak sidelobe outside main beam
    sidelobes = np.concatenate([
        E_plane[:max(0, peak_idx-half_bw)],
        E_plane[min(len(E_plane), peak_idx+half_bw):]
    ])

    return np.max(sidelobes) if len(sidelobes) > 0 else 0

# Optimize using genetic algorithm
optimized_geom = pa.thin_array_genetic_algorithm(
    full_geom,
    n_target=300,  # Keep 300 of 576 elements
    objective_func=peak_sidelobe_objective,
    population_size=50,
    n_generations=100,
    mutation_rate=0.1,
    seed=42
)

print(f"Original elements: {full_geom.n_elements}")
print(f"Optimized elements: {optimized_geom.n_elements}")

# Compare results
w_full = np.ones(full_geom.n_elements)
w_opt = np.ones(optimized_geom.n_elements)

_, E_full, _ = pa.compute_pattern_cuts(full_geom.x, full_geom.y, w_full, k)
_, E_opt, _ = pa.compute_pattern_cuts(optimized_geom.x, optimized_geom.y, w_opt, k)

print(f"Full array SLL: {np.max(E_full[100:]):.1f} dB")
print(f"Optimized SLL: {np.max(E_opt[100:]):.1f} dB")

Adaptive Beamforming Simulation#

Problem: Simulate MVDR (Capon) beamformer with interference.

def mvdr_beamformer(geom, k, theta_desired, phi_desired, interference_dirs,
                    snr_dB=20, inr_dB=30):
    """
    Minimum Variance Distortionless Response beamformer.

    Parameters
    ----------
    interference_dirs : list of (theta, phi) tuples
        Interference directions
    snr_dB : float
        Signal-to-noise ratio
    inr_dB : float
        Interference-to-noise ratio
    """
    n = geom.n_elements

    # Steering vector for desired signal
    a_d = pa.steering_vector(k, geom.x, geom.y, theta_desired, phi_desired)

    # Build interference-plus-noise covariance matrix
    sigma_n = 1.0  # Noise power (normalized)
    sigma_s = sigma_n * 10 ** (snr_dB / 10)
    sigma_i = sigma_n * 10 ** (inr_dB / 10)

    R = sigma_n * np.eye(n, dtype=complex)

    for theta_i, phi_i in interference_dirs:
        a_i = pa.steering_vector(k, geom.x, geom.y, theta_i, phi_i)
        R += sigma_i * np.outer(a_i, a_i.conj())

    # MVDR weights: w = R^-1 * a / (a^H * R^-1 * a)
    R_inv = np.linalg.inv(R)
    w = R_inv @ a_d
    w = w / (a_d.conj() @ R_inv @ a_d)

    return w

# Example
geom = pa.create_rectangular_array(16, 16, dx=0.5, dy=0.5)
k = pa.wavelength_to_k(1.0)

# Signal at 20 deg, interference at 45 and 60 deg
w_mvdr = mvdr_beamformer(
    geom, k,
    theta_desired=20, phi_desired=0,
    interference_dirs=[(45, 0), (60, 0)],
    snr_dB=20, inr_dB=40
)

# Compare with conventional
w_conv = pa.steering_vector(k, geom.x, geom.y, 20, 0)

theta_deg, E_conv, _ = pa.compute_pattern_cuts(geom.x, geom.y, w_conv, k)
_, E_mvdr, _ = pa.compute_pattern_cuts(geom.x, geom.y, w_mvdr, k)

# MVDR will have deep nulls at interference locations

Custom Array Geometries#

Problem: Create a non-standard array geometry.

from phased_array.geometry import ArrayGeometry

def create_log_periodic_array(n_elements, min_spacing, ratio):
    """
    Create a log-periodic linear array.

    Spacing increases geometrically: d_n = min_spacing * ratio^n
    """
    x = [0]
    for i in range(n_elements - 1):
        spacing = min_spacing * (ratio ** i)
        x.append(x[-1] + spacing)

    x = np.array(x)
    x -= x.mean()  # Center at origin

    return ArrayGeometry(
        x=x,
        y=np.zeros_like(x),
        z=np.zeros_like(x),
        nx=np.zeros_like(x),
        ny=np.zeros_like(x),
        nz=np.ones_like(x),
        element_indices=np.arange(len(x))
    )

# Create and test
log_geom = create_log_periodic_array(20, min_spacing=0.3, ratio=1.1)
print(f"Elements: {log_geom.n_elements}")
print(f"Total aperture: {log_geom.x.max() - log_geom.x.min():.2f} wavelengths")

Mutual Coupling Compensation#

Problem: Pre-compensate weights for known mutual coupling.

def compensate_mutual_coupling(weights, coupling_matrix):
    """
    Pre-distort weights to compensate for mutual coupling.

    The coupled weights are: w_coupled = C @ w
    We want: C @ w_predistorted = w_desired
    So: w_predistorted = C^-1 @ w_desired
    """
    C_inv = np.linalg.inv(coupling_matrix)
    return C_inv @ weights

# Create coupling matrix
geom = pa.create_rectangular_array(8, 8, dx=0.5, dy=0.5)
coupling = pa.mutual_coupling_matrix_theoretical(
    geom.x, geom.y,
    coupling_coefficient=0.3,
    coupling_exponent=2.0
)

# Desired steering
k = pa.wavelength_to_k(1.0)
w_desired = pa.steering_vector(k, geom.x, geom.y, 30, 0)
w_desired *= pa.taylor_taper_2d(8, 8, sidelobe_dB=-30)

# Compensated weights
w_comp = compensate_mutual_coupling(w_desired, coupling)

# Verify: coupling @ w_comp should give pattern close to desired
w_result = coupling @ w_comp

_, E_desired, _ = pa.compute_pattern_cuts(geom.x, geom.y, w_desired, k)
_, E_result, _ = pa.compute_pattern_cuts(geom.x, geom.y, w_result, k)

# E_result should be close to E_desired

Statistical Analysis with Monte Carlo#

Problem: Characterize performance statistics with manufacturing tolerances.

def monte_carlo_analysis(geom, k, weights, n_trials=1000,
                         pos_error_std=0.01,  # wavelengths
                         amp_error_std=0.5,   # dB
                         phase_error_std=3.0  # degrees
                         ):
    """
    Monte Carlo analysis of array performance with random errors.
    """
    results = {
        'sll': [],
        'hpbw': [],
        'beam_pointing': [],
        'directivity_loss': []
    }

    for _ in range(n_trials):
        # Position errors
        x_err = geom.x + pos_error_std * np.random.randn(geom.n_elements)
        y_err = geom.y + pos_error_std * np.random.randn(geom.n_elements)

        # Amplitude and phase errors
        amp_err = 10 ** (amp_error_std * np.random.randn(geom.n_elements) / 20)
        phase_err = np.deg2rad(phase_error_std * np.random.randn(geom.n_elements))
        w_err = weights * amp_err * np.exp(1j * phase_err)

        # Compute pattern
        theta_deg, E_plane, _ = pa.compute_pattern_cuts(
            x_err, y_err, w_err, k, n_points=721
        )

        # Extract metrics
        peak_idx = np.argmax(E_plane)
        results['beam_pointing'].append(theta_deg[peak_idx])
        results['hpbw'].append(pa.compute_half_power_beamwidth(theta_deg, E_plane))

        sidelobes = E_plane[np.abs(theta_deg) > 15]
        results['sll'].append(np.max(sidelobes))

    return {k: np.array(v) for k, v in results.items()}

# Run analysis
geom = pa.create_rectangular_array(16, 16, dx=0.5, dy=0.5)
k = pa.wavelength_to_k(1.0)
weights = pa.steering_vector(k, geom.x, geom.y, 20, 0)
weights *= pa.taylor_taper_2d(16, 16, sidelobe_dB=-35)

stats = monte_carlo_analysis(geom, k, weights, n_trials=500)

print("Performance Statistics (500 trials)")
print("-" * 40)
print(f"SLL: {np.mean(stats['sll']):.1f} dB (std: {np.std(stats['sll']):.1f})")
print(f"HPBW: {np.mean(stats['hpbw']):.2f} deg (std: {np.std(stats['hpbw']):.2f})")
print(f"Pointing: {np.mean(stats['beam_pointing']):.2f} deg (std: {np.std(stats['beam_pointing']):.2f})")

Frequency-Selective Surface Integration#

Problem: Model pattern with frequency-selective surface (FSS) radome.

def apply_fss_radome(pattern_dB, theta, fss_transmission_dB_func):
    """
    Apply FSS radome transmission function to pattern.

    Parameters
    ----------
    fss_transmission_dB_func : callable
        Function(theta_deg) -> transmission in dB
    """
    theta_deg = np.rad2deg(theta)
    transmission = fss_transmission_dB_func(theta_deg)
    return pattern_dB + transmission

# Example FSS model: bandpass with angle-dependent cutoff
def fss_model(theta_deg):
    """Simple FSS model: transmission degrades with scan angle."""
    # Normal incidence: 0 dB loss
    # Grazing incidence: -10 dB loss
    return -0.002 * theta_deg ** 2

geom = pa.create_rectangular_array(16, 16, dx=0.5, dy=0.5)
k = pa.wavelength_to_k(1.0)
weights = pa.steering_vector(k, geom.x, geom.y, 45, 0)

theta, phi, pattern_dB = pa.compute_full_pattern(geom.x, geom.y, weights, k)

# Apply FSS
for i, t in enumerate(theta):
    pattern_dB[i, :] = apply_fss_radome(
        pattern_dB[i, :], np.array([t] * len(phi)), fss_model
    )