Grow E-ELT subaperture packing

Author: Ben Mazin

nullsim/physics/pupil_geometry.py

@@ -788,6 +788,106 @@ def _evaluate_placement(
     return float(-log_dist_sum)
 
 
+def _max_feasible_diameter_for_positions(
+    positions: NDArray[np.float64],
+    pupil: PupilGeometry,
+    diameter_m: float,
+) -> float:
+    """Largest sub-aperture diameter allowed by fixed center positions."""
+
+    positions = np.asarray(positions, dtype=np.float64)
+    max_diameter = np.inf
+    r_obs = pupil.central_obstruction_diameter_m / 2
+    r_outer = diameter_m / 2
+
+    for x, y in positions:
+        radial = float(np.hypot(x, y))
+        if pupil.boundary_vertices is not None:
+            edge_clearance = _signed_distance_to_polygon_edge(
+                float(x), float(y), pupil.boundary_vertices
+            )
+            max_diameter = min(max_diameter, 2.0 * edge_clearance)
+        else:
+            max_diameter = min(max_diameter, 2.0 * (r_outer - radial))
+
+        if r_obs > 0.0:
+            max_diameter = min(max_diameter, 2.0 * (radial - r_obs))
+
+        for vane in pupil.spider_vanes:
+            vane_clearance = point_to_segment_distance(
+                float(x),
+                float(y),
+                vane.x_start,
+                vane.y_start,
+                vane.x_end,
+                vane.y_end,
+            ) - vane.width_m / 2
+            max_diameter = min(max_diameter, 2.0 * vane_clearance)
+
+    if positions.shape[0] > 1:
+        diffs = positions[:, None, :] - positions[None, :, :]
+        distances = np.sqrt((diffs ** 2).sum(axis=-1))
+        distances = distances[np.triu_indices(positions.shape[0], k=1)]
+        max_diameter = min(max_diameter, float(np.min(distances)))
+
+    if not np.isfinite(max_diameter):
+        return 0.0
+    return max(0.0, float(max_diameter))
+
+
+def _grow_refined_placement(
+    positions: NDArray[np.float64],
+    d_sub: float,
+    d_sub_hi: float,
+    n_sub: int,
+    pupil: PupilGeometry,
+    diameter_m: float,
+) -> tuple[NDArray[np.float64], float]:
+    """Grow a feasible placement using hard-constraint refinement.
+
+    The seed search is intentionally conservative. For crowded pupils with
+    spider vanes, such as E-ELT N=48, a feasible seed can often be enlarged by
+    jointly relaxing all sub-aperture centers. This pass keeps the binary
+    search's proven feasible layout as the starting point and only accepts a
+    larger diameter when the existing hard-constraint score remains feasible.
+    """
+
+    if pupil.boundary_vertices is not None or not pupil.spider_vanes or n_sub < 24:
+        return np.asarray(positions, dtype=np.float64), float(d_sub)
+
+    best_positions = np.asarray(positions, dtype=np.float64)
+    best_d_sub = float(d_sub)
+    low = best_d_sub
+    high = max(float(d_sub_hi), low)
+    if high <= low:
+        return best_positions, best_d_sub
+
+    for _ in range(12):
+        trial_d_sub = 0.5 * (low + high)
+        trial_r_sub = 0.5 * trial_d_sub
+        refined = _refine_positions(
+            best_positions,
+            n_sub,
+            trial_r_sub,
+            pupil,
+            diameter_m,
+        )
+        max_diameter = _max_feasible_diameter_for_positions(
+            refined,
+            pupil,
+            diameter_m,
+        )
+        accepted_d_sub = min(trial_d_sub, max_diameter)
+        if accepted_d_sub > best_d_sub:
+            best_positions = refined
+            best_d_sub = accepted_d_sub
+            low = accepted_d_sub
+        else:
+            high = trial_d_sub
+
+    return best_positions, best_d_sub
+
+
 def _place_single_subaperture(
     diameter_m: float, r_sub: float, pupil: PupilGeometry,
 ) -> NDArray[np.float64]:
@@ -852,6 +952,7 @@ def optimize_subaperture_positions(
 
     ring_d_max = _max_ring_d_sub(n_sub, pupil, diameter_m)
     d_sub_hi = max(diameter_m / np.sqrt(n_sub), ring_d_max) * 1.1
+    d_sub_upper_bound = d_sub_hi
     d_sub_lo = 0.1
     best_d_sub = d_sub_lo
 
@@ -919,6 +1020,15 @@ def optimize_subaperture_positions(
             f"diameter={diameter_m:.1f}m"
         )
 
+    best_result, best_d_sub = _grow_refined_placement(
+        best_result,
+        best_d_sub,
+        d_sub_upper_bound,
+        n_sub,
+        pupil,
+        diameter_m,
+    )
+
     return best_result, best_d_sub
 
 

tests/physics/test_pupil_geometry.py

@@ -274,3 +274,19 @@ def test_optimize_eelt_n18_returns_feasible_subapertures() -> None:
     r_sub = d_sub / 2.0
     for x, y in positions:
         assert is_valid_position(x, y, r_sub, pupil, EELT_DIAMETER_M)
+
+
+def test_optimize_eelt_n48_grows_packed_layout() -> None:
+    pupil = eelt_pupil()
+    positions, d_sub = optimize_subaperture_positions(48, EELT_DIAMETER_M, pupil)
+    assert positions.shape == (48, 2)
+    assert d_sub > 4.5
+
+    r_sub = d_sub / 2.0
+    for x, y in positions:
+        assert is_valid_position(x, y, r_sub, pupil, EELT_DIAMETER_M)
+
+    diffs = positions[:, np.newaxis, :] - positions[np.newaxis, :, :]
+    distances = np.sqrt(np.sum(diffs * diffs, axis=2))
+    distances = distances[np.triu_indices(positions.shape[0], k=1)]
+    assert float(np.min(distances)) >= d_sub - 1.0e-8