diff --git a/CHANGELOG.md b/CHANGELOG.md
index d82d6a2..ad858c8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,29 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [2.0.1] - 2026-06-03
+
+### Added
+
+- Concavity penalty for the L2 annotation cost function. The optimizer can occasionally
+  converge to a physically wrong, inverted arc — a failure mode on off-apex images or with
+  loose parameter bounds. The penalty detects this by checking whether the predicted arc
+  curves in the correct direction relative to its chord, and adds a large cost when it does
+  not. Enabled by default; configurable via `concavity_penalty` (bool) and
+  `concavity_penalty_scale` (float) kwargs in `fit_arc()`. The correct curvature direction
+  is auto-detected from `theta_z`, so non-standard viewing geometries (e.g. camera pointing
+  away from nadir) are handled correctly without any user configuration.
+
+### Changed
+
+- Sagitta pre-fitter (`fit_sagitta`) upgraded from a 2-D to a 3-D optimisation: the arc
+  apex x-position is now a free parameter alongside curvature and camera angle. This removes
+  a systematic radius bias on asymmetric arcs whose true apex falls outside the annotated
+  region.
+- Sagitta pre-fitter now seeds the differential-evolution warm-start with a geometrically
+  computed camera angle (`arccos(r / (r + h))`) rather than the raw optimiser estimate,
+  avoiding a near-degeneracy that could misguide the final fit.
+
 ## [2.0.0] - 2026-05-07
 
 ### Added
diff --git a/docs/api.rst b/docs/api.rst
index b70b2df..ce196be 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -418,6 +418,10 @@ Complete workflow class for limb-based planetary radius determination. Default d
 * Multiple uncertainty estimation methods (population spread, Hessian, profile likelihood)
 * Flexible cost functions (gradient-field flux, L1, L2, log-L1)
 * Support for multiple minimizers (differential-evolution, dual-annealing, basinhopping)
+* Concavity penalty on the L2 cost that prevents convergence to wrong-curvature (inverted) arc solutions;
+  configurable via ``concavity_penalty`` and ``concavity_penalty_scale`` in :meth:`fit_arc`
+* Sagitta pre-stage warm-start with 3-D optimizer (``fit_sagitta``) that seeds the differential-evolution
+  search with a geometrically correct ``theta_x`` estimate via :func:`~planet_ruler.geometry.limb_camera_angle`
 
 MaskSegmenter
 ~~~~~~~~~~~~~
diff --git a/docs/science/methodology.rst b/docs/science/methodology.rst
index c049ffa..ae3bb97 100644
--- a/docs/science/methodology.rst
+++ b/docs/science/methodology.rst
@@ -57,16 +57,41 @@ then search parameter space to minimize this cost.
 Parameter Bounds and Optimization
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-Physical constraints provide bounds on free parameters. Earth's radius is approximately 6371 km 
-(±50 km for various reference models). Altitude for aircraft photography ranges from 5-20 km. 
-Camera pitch typically ranges from -90° to +90°, roll from -180° to +180°. Focal length can 
+Physical constraints provide bounds on free parameters. Earth's radius is approximately 6371 km
+(±50 km for various reference models). Altitude for aircraft photography ranges from 5-20 km.
+Camera pitch typically ranges from -90° to +90°, roll from -180° to +180°. Focal length can
 often be constrained from EXIF metadata to within ±10%.
 
-Global optimization algorithms (differential evolution, basin-hopping, dual annealing) search this 
-bounded parameter space to find the best-fit solution. Multi-resolution strategies for gradient-field 
-optimization start with coarse image resolution (fast evaluation) and progressively refine to full 
+Global optimization algorithms (differential evolution, basin-hopping, dual annealing) search this
+bounded parameter space to find the best-fit solution. Multi-resolution strategies for gradient-field
+optimization start with coarse image resolution (fast evaluation) and progressively refine to full
 resolution, helping avoid local minima in the cost landscape.
 
+Cost Function Structure and Curvature Disambiguation
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The L2 annotation cost function measures the mean squared pixel distance between predicted
+and observed horizon points. However, because a planetary limb arc and its mirror image
+(with the curvature flipped) can produce nearly identical pixel residuals, the optimizer
+can sometimes converge to a physically wrong solution: a predicted arc that curves in the
+opposite direction from the true limb.
+
+Planet Ruler resolves this ambiguity with an additive *concavity penalty* that fires when
+the predicted arc curves the wrong way. The penalty evaluates the predicted arc at the
+annotated x-coordinates, draws a chord between the leftmost and rightmost predictions,
+and checks whether the interior of the arc lies above or below that chord:
+
+* **∩-shaped arc** (interior above the chord): the limb peaks toward the top of the image —
+  the correct shape for standard high-altitude photography (e.g. an airplane window).
+  Penalty = 0.
+* **∪-shaped arc** (interior below the chord): the limb sags toward the bottom of the image
+  — the wrong-curvature mirror solution. Penalty fires, adding a large cost that steers
+  differential evolution out of that basin.
+
+The penalty is enabled by default. It can be disabled with ``concavity_penalty=False`` in
+:meth:`~planet_ruler.observation.LimbObservation.fit_arc` if needed (for example, when
+the viewing geometry is inverted and the ∪ shape is genuinely correct).
+
 Uncertainty Quantification
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/docs/tutorial_airplane.rst b/docs/tutorial_airplane.rst
index d4d1180..1190793 100644
--- a/docs/tutorial_airplane.rst
+++ b/docs/tutorial_airplane.rst
@@ -610,6 +610,31 @@ Use the crop tool to remove obstructions:
    # Continue with normal workflow
    obs = pr.LimbObservation("airplane_photo_cropped.jpg", config)
 
+**"The fitted arc curves the wrong way (arc sags downward instead of peaking upward)"**
+
+The optimizer occasionally locks onto a mirror-image solution where the predicted horizon
+curves in the opposite direction from the true limb. Planet Ruler includes a *concavity penalty*
+that prevents this by default, but if you see a clearly wrong fit you can try:
+
+.. code-block:: python
+
+   # Increase the penalty weight (default is 100) to push harder out of the wrong basin
+   obs.fit_arc(
+       minimizer='differential-evolution',
+       max_iter=1000,
+       concavity_penalty_scale=500,
+   )
+
+   # Or re-run with a different random seed
+   obs.fit_arc(minimizer='differential-evolution', max_iter=1000, seed=99)
+
+If you are photographing from an unusual angle where the horizon genuinely curves downward
+(rare outside of ISS or near-vertical imagery), disable the penalty:
+
+.. code-block:: python
+
+   obs.fit_arc(concavity_penalty=False)
+
 **"Result varies between photos"**
 
 Normal! Try:
diff --git a/notebooks/tutorials/measure_your_planet.ipynb b/notebooks/tutorials/measure_your_planet.ipynb
index 7b4d142..35331ba 100644
--- a/notebooks/tutorials/measure_your_planet.ipynb
+++ b/notebooks/tutorials/measure_your_planet.ipynb
@@ -149,12 +149,14 @@
     "gps_altitude = get_gps_altitude(image_path)\n",
     "if gps_altitude:\n",
     "    print(f\"\\n✓ GPS Altitude Found: {gps_altitude:.1f} meters\")\n",
-    "    print(f\"   Using this instead of your specified altitude: {altitude_m} meters\")\n",
     "    altitude_m = gps_altitude\n",
     "else:\n",
-    "    print(\n",
-    "        f\"\\nGPS Altitude: Not found - using your specified altitude of {altitude_m} meters\"\n",
-    "    )\n",
+    "    try:\n",
+    "        print(\n",
+    "            f\"\\nGPS Altitude: Not found - using your specified altitude of {altitude_m} meters\"\n",
+    "        )\n",
+    "    except:\n",
+    "        print(\"No GPS Altitude found - specify altitude_m or try a different image\")\n",
     "\n",
     "print(\"\\n\" + \"=\" * 60)"
    ]
diff --git a/planet_ruler/fit.py b/planet_ruler/fit.py
index fa59028..9fbdd68 100644
--- a/planet_ruler/fit.py
+++ b/planet_ruler/fit.py
@@ -25,7 +25,12 @@
     bilinear_interpolate,
     gradient_field,
 )
-from planet_ruler.geometry import _r_from_K, limb_arc_sagitta, limb_arc
+from planet_ruler.geometry import (
+    _r_from_K,
+    limb_arc_sagitta,
+    limb_arc,
+    limb_camera_angle,
+)
 from typing import Callable
 import logging
 
@@ -133,12 +138,23 @@ def __init__(
         free_parameters: list,
         init_parameter_values: dict,
         loss_function: str = "l2",
+        concavity_penalty: bool = True,
+        concavity_penalty_scale: float = 100.0,
     ):
         _VALID = {"l2", "l1", "log-l1"}
         if loss_function not in _VALID:
             raise ValueError(f"Unrecognized loss function: {loss_function!r}")
         super().__init__(target, function, free_parameters, init_parameter_values)
         self.loss_function = loss_function
+        self.concavity_penalty = concavity_penalty
+        self.concavity_penalty_scale = concavity_penalty_scale
+        # Auto-detect expected arc curvature from theta_z.
+        # When theta_z ≈ π (standard airplane/smartphone orientation) the arc
+        # is ∩-shaped; penalise ∪.  When theta_z ≈ 0 (some ISS / telescope
+        # configurations) the arc is ∪-shaped; penalise ∩ instead (invert).
+        # cos(theta_z) > 0 ↔ theta_z closer to 0 than to π ↔ invert.
+        theta_z = float(init_parameter_values.get("theta_z", np.pi))
+        self._invert_concavity = bool(np.cos(theta_z) > 0)
         valid_mask = ~np.isnan(target)
         self.x = np.where(valid_mask)[0]
         self.target = target[valid_mask]
@@ -154,7 +170,18 @@ def cost(self, params: np.ndarray | dict) -> float:
         if self.is_sparse and y is not None and len(y) != len(self.target):
             y = y[self.x]
         if self.loss_function == "l2":
-            return np.mean((y - self.target) ** 2)
+            l2 = float(np.mean((y - self.target) ** 2))
+            penalty = (
+                _chord_sagitta_penalty(
+                    y,
+                    self.x,
+                    self.concavity_penalty_scale,
+                    invert=self._invert_concavity,
+                )
+                if self.concavity_penalty
+                else 0.0
+            )
+            return l2 + penalty
         elif self.loss_function == "l1":
             return np.mean(np.abs(y - self.target))
         elif self.loss_function == "log-l1":
@@ -163,6 +190,96 @@ def cost(self, params: np.ndarray | dict) -> float:
         raise ValueError(f"Unrecognized loss function: {self.loss_function}")
 
 
+def _chord_sagitta_penalty(
+    y: np.ndarray,
+    x: np.ndarray,
+    scale: float = 100.0,
+    invert: bool = False,
+) -> float:
+    """Penalise an incorrectly curved arc prediction in the L2 cost function.
+
+    For a correctly curved planetary limb in the standard orientation
+    (theta_z ≈ π), the arc peaks toward the top of the image (∩ shape in
+    image coordinates where y increases downward): the arc centre has a
+    *smaller* y-value than the chord connecting the leftmost and rightmost
+    annotation points.
+
+    When the optimiser converges to the mirror basin the L2 residuals can be
+    arbitrarily small yet the geometry is physically wrong.  This penalty
+    detects that case and adds a large additive cost so differential evolution
+    escapes.
+
+    The ``invert`` flag handles viewing geometries where the correct arc is
+    ∪-shaped (e.g. theta_z ≈ 0 for some ISS or telescope configurations).
+    L2CostFunction sets this automatically from theta_z in
+    init_parameter_values.
+
+    Algorithm:
+        1. Guard: fewer than 3 points or any non-finite prediction → 0.
+        2. Guard: max(y) − min(y) < 0.5 px (flat proxy-mode arc, entirely
+           outside the FOV) → 1e6 (heavy penalty for degenerate basin).
+        3. Sort annotation points by x; define chord between leftmost and
+           rightmost (x, y_predicted) pairs.
+        4. For each point compute d = y_predicted − chord_y.
+           d > 0: arc is below chord (∪ shape).
+           d < 0: arc is above chord (∩ shape).
+        5. sum_raw = Σ d.
+           invert=False: penalise if sum_raw > 0  (∪ shape, wrong for θz≈π).
+           invert=True:  penalise if sum_raw < 0  (∩ shape, wrong for θz≈0).
+           Penalty = scale × mean(d²).
+
+    Args:
+        y: Predicted y-pixel values at the annotated x-coordinates.
+        x: Annotated x-pixel coordinates (parallel to y).
+        scale: Multiplicative weight applied to mean(d²) when the arc is in
+            the wrong basin.  Defaults to 100.0, which keeps the penalty
+            roughly two orders of magnitude above the typical L2 cost (O(1)
+            px²) so the optimiser reliably avoids the wrong basin.
+        invert: If False (default), penalise ∪-shaped arcs (correct when
+            theta_z ≈ π).  If True, penalise ∩-shaped arcs (correct when
+            theta_z ≈ 0).
+
+    Returns:
+        Non-negative penalty term to be added to the L2 cost.
+    """
+    n = len(x)
+    if n < 3:
+        return 0.0
+    if not np.all(np.isfinite(y)):
+        return 0.0
+
+    # Flat arc: limbArc returns a near-constant when the arc is entirely
+    # outside the FOV, creating a false basin with finite but meaningless cost.
+    if float(np.max(y) - np.min(y)) < 0.5:
+        return 1e6
+
+    order = np.argsort(x)
+    xs = x[order]
+    ys = y[order]
+
+    x_l, y_l = float(xs[0]), float(ys[0])
+    x_r, y_r = float(xs[-1]), float(ys[-1])
+    dx = x_r - x_l
+    if dx < 1e-6:
+        return 0.0
+
+    # d > 0: arc is below the chord in image space (∪ shape).
+    # d < 0: arc is above the chord (∩ shape).
+    chord_y = y_l + (y_r - y_l) * (xs - x_l) / dx
+    d = ys - chord_y
+    sum_raw = float(np.sum(d))
+
+    if invert:
+        # Expected arc is ∪-shaped (theta_z ≈ 0): penalise ∩ (sum_raw < 0).
+        if sum_raw >= 0.0:
+            return 0.0
+    else:
+        # Expected arc is ∩-shaped (theta_z ≈ π): penalise ∪ (sum_raw > 0).
+        if sum_raw <= 0.0:
+            return 0.0
+    return scale * float(np.mean(d**2))
+
+
 class GradientFieldCostFunction(BaseCostFunction):
     """Cost function using gradient field flux alignment."""
 
@@ -383,6 +500,8 @@ def __init__(
         directional_decay_rate: float = 0.15,
         prefer_direction=None,
         n_jobs: int = 1,
+        concavity_penalty: bool = True,
+        concavity_penalty_scale: float = 100.0,
     ):
         self.target = target
         self.free_parameters = free_parameters
@@ -399,6 +518,8 @@ def __init__(
         self.directional_decay_rate = directional_decay_rate
         self.prefer_direction = prefer_direction
         self.n_jobs = n_jobs
+        self.concavity_penalty = concavity_penalty
+        self.concavity_penalty_scale = concavity_penalty_scale
 
     def fit(self) -> dict:
         working_parameters = dict(self.init_parameter_values)
@@ -421,6 +542,8 @@ def fit(self) -> dict:
                 free_parameters=self.free_parameters,
                 init_parameter_values=working_parameters,
                 loss_function=self.loss_function,
+                concavity_penalty=self.concavity_penalty,
+                concavity_penalty_scale=self.concavity_penalty_scale,
             )
         else:
             cost_fn = GradientFieldCostFunction(
@@ -846,9 +969,8 @@ def fit(self) -> dict:
 
         xs = np.where(~np.isnan(self.limb))[0].astype(float)
         ys = self.limb[~np.isnan(self.limb)].astype(float)
-        u = xs - float(est["x_apex"])
-        s = ys - float(est["y_apex"])
 
+        x_apex0 = float(est["x_apex"])
         K0 = float(est["K"]) if np.isfinite(est["K"]) else 100.0
         tx0 = float(est["theta_x_est"]) if np.isfinite(est["theta_x_est"]) else 0.0
         log_kappa0 = np.log(1.0 / max(K0, 1e-12))
@@ -857,33 +979,60 @@ def fit(self) -> dict:
         log_kappa_bounds = (-3.0, 3.0)
         tx_bounds = (float(tx_lo), float(tx_hi))
 
-        # Step 2: 2-D L-BFGS-B over [log_kappa, theta_x]
-        def cost_2d(params: np.ndarray) -> float:
-            log_kap, tx = params
+        x_min, x_max = float(xs.min()), float(xs.max())
+        x_range = x_max - x_min
+        # Allow apex to range 3× the annotation width beyond either edge so that
+        # asymmetric arcs with an off-frame apex are handled correctly.
+        x_apex_bounds = (x_min - 3.0 * x_range, x_max + 3.0 * x_range)
+
+        # Step 2: 3-D L-BFGS-B over [log_kappa, theta_x, x_apex].
+        # Freeing x_apex avoids a systematic bias when the true arc apex lies
+        # outside the annotated region.  The y-apex offset is absorbed by the
+        # floating intercept s0 computed analytically inside the cost.
+        def cost_3d(params: np.ndarray) -> float:
+            log_kap, tx, x_apex_free = params
             r_t = _r_from_K(1.0 / np.exp(log_kap), self.h)
+            u_free = xs - x_apex_free
             s_pred = np.array(
-                [limb_arc_sagitta(float(ui), tx, self.f_px, r_t, self.h) for ui in u]
+                [
+                    limb_arc_sagitta(float(ui), tx, self.f_px, r_t, self.h)
+                    for ui in u_free
+                ]
             )
-            s0 = float(np.mean(s - s_pred))
-            return float(np.sum((s - s0 - s_pred) ** 2))
+            s0 = float(np.mean(ys - s_pred))
+            return float(np.sum((ys - s0 - s_pred) ** 2))
 
         try:
-            res2d = minimize(
-                cost_2d,
-                x0=[log_kappa0, tx0],
+            res3d = minimize(
+                cost_3d,
+                x0=[log_kappa0, tx0, x_apex0],
                 method="L-BFGS-B",
-                bounds=[log_kappa_bounds, tx_bounds],
+                bounds=[log_kappa_bounds, tx_bounds, x_apex_bounds],
                 options={"maxiter": 500, "ftol": 1e-12},
             )
-            K_opt = 1.0 / float(np.exp(res2d.x[0]))
-            theta_x_opt = float(res2d.x[1])
+            K_opt = 1.0 / float(np.exp(res3d.x[0]))
+            theta_x_opt = float(res3d.x[1])
+            x_apex_opt = float(res3d.x[2])
         except Exception as exc:
-            warn.append(f"SagittaFitter 2-D minimisation failed: {exc}")
+            warn.append(f"SagittaFitter 3-D minimisation failed: {exc}")
             K_opt = K0
             theta_x_opt = tx0
+            x_apex_opt = x_apex0
 
         r_opt = float(_r_from_K(K_opt, self.h))
 
+        # Recompute u and s in the optimised x_apex frame for the jackknife loop.
+        # The jackknife re-derives its own floating intercept s0 internally, so
+        # only the coordinate frame of u must be consistent.
+        u = xs - x_apex_opt
+        s_pred_opt = np.array(
+            [
+                limb_arc_sagitta(float(ui), theta_x_opt, self.f_px, r_opt, self.h)
+                for ui in u
+            ]
+        )
+        s = ys - float(np.mean(ys - s_pred_opt))
+
         # Step 3: uncertainty bounds — jackknife (theta_x fixed) or OLS fallback
         K_sigma_ols = float(est.get("K_sigma", 0.0))
         K_sigma_jack = float("nan")
@@ -953,7 +1102,16 @@ def cost_1d(log_kap: float) -> float:
             "theta_x_est": theta_x_opt,
             "theta_x_sigma": float("nan"),
             "K_sigma_jack": K_sigma_jack,
-            "updated_init": {"r": r_opt, "theta_x": theta_x_opt},
+            # Use the geometric dip angle as the theta_x warm-start rather than
+            # theta_x_opt from the sagitta fit.  limb_arc_sagitta is nearly
+            # degenerate between alpha and pi-alpha for small kappa*sin(theta_x)
+            # values typical at cruising altitude, so the 3-D minimiser may
+            # converge to the wrong basin.  limb_camera_angle(r, h) = arccos(r/(r+h))
+            # is always the correct geometric estimate given a reliable r.
+            "updated_init": {
+                "r": r_opt,
+                "theta_x": limb_camera_angle(r_opt, self.h),
+            },
             "updated_limits": {"r": [r_low, r_high]},
             "status": "ok",
             "warnings": warn,
diff --git a/planet_ruler/observation.py b/planet_ruler/observation.py
index 4dd3b91..1366b4a 100644
--- a/planet_ruler/observation.py
+++ b/planet_ruler/observation.py
@@ -1025,6 +1025,8 @@ def fit_arc(
         seed: int = 0,
         verbose: bool = False,
         n_jobs: int = 1,
+        concavity_penalty: bool = True,
+        concavity_penalty_scale: float = 100.0,
         _dashboard=None,
     ) -> "LimbObservation":
         """
@@ -1044,6 +1046,16 @@ def fit_arc(
             verbose: Print progress.
             n_jobs: Parallel workers. Effective only for differential-evolution
                 and shgo; emits a UserWarning for other minimizers.
+            concavity_penalty: If True (default), add a chord-relative sagitta
+                penalty to the l2 cost that discourages ∪-shaped (inverted) arc
+                solutions.  The penalty fires only when the predicted arc sags
+                below the chord connecting the leftmost and rightmost annotation
+                points, which is geometrically wrong for a planetary limb.
+                Disable if the optimizer is reliably finding the correct ∩-shaped
+                basin and the penalty overhead is unwanted.
+            concavity_penalty_scale: Multiplicative weight applied to the
+                mean-squared chord deviation when the arc is inverted.  Defaults
+                to 100.0.  Has no effect when concavity_penalty is False.
             _dashboard: Internal — FitDashboard instance passed by fit_limb.
 
         Returns:
@@ -1102,6 +1114,8 @@ def _arc_db_callback(xk, convergence=None, *args, **kwargs):
             seed=seed,
             verbose=verbose,
             n_jobs=n_jobs,
+            concavity_penalty=concavity_penalty,
+            concavity_penalty_scale=concavity_penalty_scale,
         )
 
         result = fitter.fit()
@@ -1125,6 +1139,8 @@ def _arc_db_callback(xk, convergence=None, *args, **kwargs):
                 "loss_function": loss_function,
                 "minimizer": effective_minimizer,
                 "minimizer_preset": minimizer_preset,
+                "concavity_penalty": concavity_penalty,
+                "concavity_penalty_scale": concavity_penalty_scale,
                 "best_parameters": self.best_parameters,
                 "fit_results": self.fit_results,
                 "status": result.get("status", "ok"),
diff --git a/tests/test_fit.py b/tests/test_fit.py
index 97d7956..a0af4f7 100644
--- a/tests/test_fit.py
+++ b/tests/test_fit.py
@@ -1562,5 +1562,131 @@ def test_scipy_minimize_n_restarts(self):
         assert result["best_parameters"] is not None
 
 
+# ============================================================================
+# CHORD-RELATIVE SAGITTA PENALTY
+# ============================================================================
+
+
+@pytest.mark.unit
+class TestChordSagittaPenalty:
+    """Unit tests for _chord_sagitta_penalty and its integration into L2CostFunction.
+
+    Image-coordinate convention (y increases downward):
+      ∩-shaped arc  (centre has smaller y than chord) = correct  = no penalty.
+      ∪-shaped arc  (centre has larger  y than chord) = inverted = penalty fires.
+    """
+
+    def _penalty(self, y, x, scale=100.0):
+        from planet_ruler.fit import _chord_sagitta_penalty
+
+        return _chord_sagitta_penalty(
+            np.asarray(y, float), np.asarray(x, float), scale=scale
+        )
+
+    def test_correct_arc_zero_penalty(self):
+        """∩-shaped arc (centre has smaller y than chord) must produce zero penalty."""
+        x = np.arange(11, dtype=float)
+        # Arc peaks toward the top of the image: centre (x=5) has y=0 (smallest);
+        # endpoints have y=2.5.  Correct ∩ shape — no penalty.
+        y = 0.1 * (x - 5.0) ** 2
+        assert self._penalty(y, x) == 0.0
+
+    def test_inverted_arc_positive_penalty(self):
+        """∪-shaped arc (centre has larger y than chord) must produce a positive penalty."""
+        x = np.arange(11, dtype=float)
+        # Arc sags toward the bottom of the image: centre (x=5) has y=2.5 (largest);
+        # endpoints have y=0.  Inverted ∪ shape — penalty fires.
+        y = -0.1 * (x - 5.0) ** 2 + 2.5
+        assert self._penalty(y, x) > 0.0
+
+    def test_proxy_mode_flat_arc(self):
+        """A near-constant arc (out-of-FOV proxy mode) must return 1e6."""
+        x = np.arange(10, dtype=float)
+        y = np.full(10, 50.0)
+        assert self._penalty(y, x) == 1e6
+
+    def test_fewer_than_three_points(self):
+        """Guard: fewer than 3 points must return 0 without raising."""
+        assert self._penalty([1.0, 2.0], [0.0, 5.0]) == 0.0
+
+    def test_penalty_scales_with_scale_param(self):
+        """Penalty magnitude must scale linearly with the scale parameter."""
+        x = np.arange(11, dtype=float)
+        y = -0.1 * (x - 5.0) ** 2 + 2.5  # inverted ∪ arc
+        p1 = self._penalty(y, x, scale=100.0)
+        p2 = self._penalty(y, x, scale=200.0)
+        assert p2 > 0.0
+        assert np.isclose(p2 / p1, 2.0, rtol=1e-9)
+
+    def test_l2_cost_adds_penalty_for_inverted_arc(self):
+        """L2CostFunction.cost() must be > 0 when the prediction is a perfect ∪ arc.
+
+        Target == prediction (zero L2 residuals) but ∪-shaped, so only the
+        penalty drives the total cost above zero.
+        """
+        x_full = np.arange(11, dtype=float)
+        y_target = -0.1 * (x_full - 5.0) ** 2 + 2.5  # ∪ shape (inverted)
+
+        def inverted_arc(**kwargs):
+            return y_target.copy()
+
+        cost_fn = L2CostFunction(
+            target=y_target,
+            function=inverted_arc,
+            free_parameters=[],
+            init_parameter_values={},
+            loss_function="l2",
+            concavity_penalty=True,
+        )
+        assert cost_fn.cost({}) > 0.0
+
+    def test_l2_cost_penalty_disabled(self):
+        """When concavity_penalty=False, a perfect-fit ∪ arc must have cost == 0."""
+        x_full = np.arange(11, dtype=float)
+        y_target = -0.1 * (x_full - 5.0) ** 2 + 2.5  # ∪ shape (inverted)
+
+        def inverted_arc(**kwargs):
+            return y_target.copy()
+
+        cost_fn = L2CostFunction(
+            target=y_target,
+            function=inverted_arc,
+            free_parameters=[],
+            init_parameter_values={},
+            loss_function="l2",
+            concavity_penalty=False,
+        )
+        assert np.isclose(cost_fn.cost({}), 0.0, atol=1e-10)
+
+    def test_l2_cost_custom_scale(self):
+        """concavity_penalty_scale is reflected in the cost value."""
+        x_full = np.arange(11, dtype=float)
+        y_target = -0.1 * (x_full - 5.0) ** 2 + 2.5  # ∪ shape (inverted)
+
+        def inverted_arc(**kwargs):
+            return y_target.copy()
+
+        cost_fn_100 = L2CostFunction(
+            target=y_target,
+            function=inverted_arc,
+            free_parameters=[],
+            init_parameter_values={},
+            loss_function="l2",
+            concavity_penalty=True,
+            concavity_penalty_scale=100.0,
+        )
+        cost_fn_200 = L2CostFunction(
+            target=y_target,
+            function=inverted_arc,
+            free_parameters=[],
+            init_parameter_values={},
+            loss_function="l2",
+            concavity_penalty=True,
+            concavity_penalty_scale=200.0,
+        )
+        # Both have zero raw L2 (perfect fit), so cost ratio == scale ratio.
+        assert cost_fn_200.cost({}) > cost_fn_100.cost({}) > 0.0
+
+
 if __name__ == "__main__":
     pytest.main([__file__, "-v"])