@@ -865,7 +865,7 @@ def calc_foclen(self):
 
         for surf in self.surfaces:
             n2 = surf.mat2.ior(wvln).double()
-            omega = omega - y * surf.paraxial_power(n, n2).double()
+            omega = omega - y * surf.paraxial_power(n, n2, wvln).double()
             y_last = y
             y = y + omega * (surf._get_effective_d_next().double() / n2)
             n = n2
@@ -1210,7 +1210,7 @@ def calc_pupil_paraxial(self, reverse=False):
                 y = y - omega * (d / n)
                 z = z - d
                 n1 = surfaces[i - 1].mat2.ior(wvln).double() if i > 0 else air
-                omega = omega + y * surfaces[i].paraxial_power(n1, n).double()
+                omega = omega + y * surfaces[i].paraxial_power(n1, n, wvln).double()
                 n = n1
         else:
             n = aper_surf.mat2.ior(wvln).double()
@@ -1220,7 +1220,7 @@ def calc_pupil_paraxial(self, reverse=False):
                 z = z + d
                 surf = surfaces[i + 1]
                 n2 = surf.mat2.ior(wvln).double()
-                omega = omega - y * surf.paraxial_power(n, n2).double()
+                omega = omega - y * surf.paraxial_power(n, n2, wvln).double()
                 n = n2
 
         # A stop conjugate to infinity (telecentric side) has no finite pupil.
 
@@ -51,7 +51,8 @@ def _resolve_zmx_glass(glass_record):
             # Preserve catalog dispersion whenever the named glass is known.
             Material(normalized_name)
             return normalized_name
-        except NotImplementedError:
+        except (NotImplementedError, ValueError):
+            # ValueError: a name with "/" that is not a numeric n/V pair.
             pass
 
         # Zemax model-glass records carry nd/Vd in fields four and five. These
@@ -472,7 +473,8 @@ def read_lens_seq(self, filename="./test.seq"):
         print("Beginning to parse surface data...\n")
 
         for line_num, line in enumerate(lines, 1):
-            line = line.strip()
+            # "!" starts a Code V comment, also after a command.
+            line = line.split("!", 1)[0].strip()
 
             # Skip irrelevant lines
             if not line or line.startswith(
@@ -535,7 +537,8 @@ def read_lens_seq(self, filename="./test.seq"):
                     surfaces.append(current_surface)
                     surface_index += 1
 
-                parts = line.split()
+                # ";" ends the S command; another command may follow it.
+                parts = line.split(";", 1)[0].split()
                 radius_value = float(parts[1]) if len(parts) > 1 else 0.0
                 thickness = float(parts[2]) if len(parts) > 2 else 0.0
                 material = parts[3].upper() if len(parts) > 3 else "AIR"
 
@@ -40,11 +40,12 @@ def prune_surf(self, mounting_margin=None):
         """Prune surface radii so all valid rays pass through, then enforce manufacturability.
 
         Traces 16 meridional fields from 0 to the full FoV to find the maximum
-        ray height [mm] on each surface, expands it by a mounting margin, then
-        caps the proposed radii to satisfy an edge-sag limit and edge-clearance
-        (air-gap / edge-thickness) constraints with neighbouring surfaces.
+        sensor-reaching ray height [mm] on each surface, expands it by a mounting
+        margin, then caps the proposed radii to satisfy an edge-sag limit and
+        edge-clearance (air-gap / edge-thickness) constraints with neighbouring surfaces.
         Aperture surfaces are not resized. The capped radii are committed via
-        each surface's `update_r`.
+        each surface's `update_r`. With no sensor-reaching ray height, the saved
+        radius is used without a mounting margin; geometry caps still apply.
 
         Args:
             mounting_margin (float or None, optional): Absolute mounting margin
@@ -75,8 +76,11 @@ def prune_surf(self, mounting_margin=None):
 
         # Ray record, shape [num_rays, num_surfaces + 2, 3]
         ray_o_record = torch.stack(ray_o_record, dim=-2)
-        ray_o_record = torch.nan_to_num(ray_o_record, nan=0.0, posinf=0.0, neginf=0.0)
         ray_o_record = ray_o_record.reshape(-1, ray_o_record.shape[-2], 3)
+        # A ray lost later must not inflate an earlier surface's aperture.
+        reaches_sensor = torch.isfinite(ray_o_record[:, -1, :]).all(-1)
+        ray_o_record = torch.where(reaches_sensor[:, None, None], ray_o_record, 0.0)
+        ray_o_record = torch.nan_to_num(ray_o_record, nan=0.0, posinf=0.0, neginf=0.0)
 
         # Compute the maximum ray height for each surface
         ray_r_record = (ray_o_record[..., :2] ** 2).sum(-1).sqrt()
@@ -91,7 +95,9 @@ def prune_surf(self, mounting_margin=None):
             if surf_r_max[i] > 0:
                 base = float(surf_r_max[i].item())
             else:
-                base = float(self.surfaces[i].r)
+                # No ray evidence: repeated pruning must not add another margin.
+                proposed_r[i] = float(saved_radii[i])
+                continue
 
             # Expand the ray-traced radius by a mounting margin
             if mounting_margin is None:
 
@@ -81,27 +81,33 @@ def psf(self, points, wvln=None, ks=PSF_KS, **kwargs):
                 recenter (bool): If True (default), center the PSF on the chief
                 ray; otherwise on the pinhole projection.
                 model (str): One of 'geometric' (default), 'coherent', 'huygens'.
+                return_field (bool): For the coherent and Huygens models, return
+                    the energy-normalized complex sensor-plane field instead of
+                    intensity. Defaults to False.
 
         Returns:
-            psf (torch.Tensor): PSF normalized to sum to 1. Shape [ks, ks] for a
-                single point, or [N, ks, ks] for the geometric model with N points.
+            torch.Tensor: Intensity PSF normalized to sum to 1, or an
+                energy-normalized complex field when ``return_field=True``.
+                Shape [ks, ks] for a single point, or [N, ks, ks] for the
+                geometric model with N points.
 
         Raises:
             ValueError: If `model` is not one of the supported names.
         """
         wvln = self.primary_wvln if wvln is None else wvln
         spp = kwargs.get("spp", None)
         recenter = kwargs.get("recenter", True)
+        return_field = kwargs.get("return_field", False)
         model = kwargs.get("model", "geometric")
         if model == "geometric":
             spp = SPP_PSF if spp is None else spp
             return self.psf_geometric(points, ks, wvln, spp, recenter)
         elif model == "coherent":
             spp = SPP_COHERENT if spp is None else spp
-            return self.psf_coherent(points, ks, wvln, spp, recenter)
+            return self.psf_coherent(points, ks, wvln, spp, recenter, return_field)
         elif model == "huygens":
             spp = SPP_COHERENT if spp is None else spp
-            return self.psf_huygens(points, ks, wvln, spp, recenter)
+            return self.psf_huygens(points, ks, wvln, spp, recenter, return_field)
         else:
             raise ValueError(f"Unknown PSF model: {model}")
 
@@ -182,7 +188,13 @@ def psf_geometric(self, points, ks=PSF_KS, wvln=None, spp=SPP_PSF, recenter=True
         return diff_float(psf)
 
     def psf_coherent(
-        self, points, ks=PSF_KS, wvln=None, spp=SPP_COHERENT, recenter=True
+        self,
+        points,
+        ks=PSF_KS,
+        wvln=None,
+        spp=SPP_COHERENT,
+        recenter=True,
+        return_field=False,
     ):
         """Compute the coherent exit-pupil PSF (alias for `psf_pupil_prop`).
 
@@ -198,22 +210,43 @@ def psf_coherent(
                 back to `self.primary_wvln`.
             spp (int, optional): Rays sampled. Defaults to SPP_COHERENT.
             recenter (bool, optional): If True (default), center on the chief ray.
+            return_field (bool, optional): Return the energy-normalized complex
+                sensor-plane field instead of intensity. Defaults to False.
 
         Returns:
-            psf (torch.Tensor): PSF normalized to sum to 1. Shape [ks, ks].
+            torch.Tensor: Intensity PSF normalized to sum to 1, or an
+                energy-normalized complex field when ``return_field=True``.
+                Shape [ks, ks].
+
+        Note:
+            The complex field has an arbitrary global phase because optical path
+            length is measured relative to a reference; relative phase is meaningful.
         """
         wvln = self.primary_wvln if wvln is None else wvln
-        return self.psf_pupil_prop(points, ks=ks, wvln=wvln, spp=spp, recenter=recenter)
+        return self.psf_pupil_prop(
+            points,
+            ks=ks,
+            wvln=wvln,
+            spp=spp,
+            recenter=recenter,
+            return_field=return_field,
+        )
 
     def psf_pupil_prop(
-        self, points, ks=PSF_KS, wvln=None, spp=SPP_COHERENT, recenter=True
+        self,
+        points,
+        ks=PSF_KS,
+        wvln=None,
+        spp=SPP_COHERENT,
+        recenter=True,
+        return_field=False,
     ):
         """Compute the single-point monochromatic PSF via the exit-pupil diffraction model.
 
         Steps:
             1. Compute the complex wavefront at the exit-pupil plane by coherent ray tracing.
             2. Propagate to the sensor plane with the Angular Spectrum Method (ASM)
-               and take the intensity as the PSF. This function is differentiable.
+               and return either its field or intensity. This function is differentiable.
 
         Args:
             points (torch.Tensor or list): Single normalized point source [3] or
@@ -226,12 +259,14 @@ def psf_pupil_prop(
             spp (int, optional): Number of rays to sample. Defaults to SPP_COHERENT.
             recenter (bool, optional): If True (default), center on the chief ray;
                 otherwise on the pinhole projection.
+            return_field (bool, optional): Return the energy-normalized complex
+                sensor-plane field instead of intensity. Defaults to False.
 
         Returns:
-            psf (torch.Tensor): PSF normalized to sum to 1. Shape [ks, ks] when
-                `ks` is given, or the full uncropped intensity field of shape
-                [1, 1, 2H, 2H] (twice the exit-pupil grid after zero-padding)
-                when `ks` is None.
+            torch.Tensor: Intensity PSF normalized to sum to 1, or an
+                energy-normalized complex field when ``return_field=True``.
+                Shape [ks, ks] when `ks` is given, or [1, 1, 2H, 2H] when
+                `ks` is None (twice the exit-pupil grid after zero-padding).
 
         Reference:
             [1] "End-to-End Hybrid Refractive-Diffractive Lens Design with Differentiable Ray-Wave Model", SIGGRAPH Asia 2024.
@@ -241,6 +276,8 @@ def psf_pupil_prop(
             Angular Spectrum Method (ASM) instead of a single FFT. ASM is more
             accurate because the FFT approach assumes a far-field condition
             (e.g., chief ray perpendicular to the image plane).
+            The complex field has an arbitrary global phase because optical path
+            length is measured relative to a reference; relative phase is meaningful.
         """
         wvln = self.primary_wvln if wvln is None else wvln
         # Pupil field by coherent ray tracing
@@ -275,6 +312,26 @@ def psf_pupil_prop(
         psfc_idx_i = ((2 - psfc[1]) * h / 4).round().long()
         psfc_idx_j = ((2 + psfc[0]) * w / 4).round().long()
 
+        if return_field:
+            field = sensor_field
+            if ks is not None:
+                field = (
+                    F.pad(
+                        field,
+                        [ks // 2, ks // 2, ks // 2, ks // 2],
+                        mode="constant",
+                        value=0,
+                    )
+                    .squeeze(0)
+                    .squeeze(0)
+                )
+                field = field[
+                    psfc_idx_i : psfc_idx_i + ks,
+                    psfc_idx_j : psfc_idx_j + ks,
+                ]
+            energy = field.abs().square().sum(dim=(-2, -1), keepdim=True)
+            return field / torch.sqrt(energy + EPSILON)
+
         # Crop valid PSF region and normalize
         if ks is not None:
             psf_inten_pad = (
@@ -392,7 +449,13 @@ def pupil_field(self, points, wvln=None, spp=SPP_COHERENT, recenter=True):
         return wavefront, psf_center
 
     def psf_huygens(
-        self, points, ks=PSF_KS, wvln=None, spp=SPP_COHERENT, recenter=True
+        self,
+        points,
+        ks=PSF_KS,
+        wvln=None,
+        spp=SPP_COHERENT,
+        recenter=True,
+        return_field=False,
     ):
         """Compute the single-wavelength Huygens PSF by spherical-wave integration.
 
@@ -414,16 +477,22 @@ def psf_huygens(
             spp (int, optional): Rays sampled. Defaults to SPP_COHERENT.
             recenter (bool, optional): If True (default), center on the chief ray;
                 otherwise on the pinhole projection.
+            return_field (bool, optional): Return the energy-normalized complex
+                sensor-plane field instead of intensity. Defaults to False.
 
         Returns:
-            psf (torch.Tensor): PSF normalized to sum to 1. Shape [ks, ks].
+            torch.Tensor: Intensity PSF normalized to sum to 1, or an
+                energy-normalized complex field when ``return_field=True``.
+                Shape [ks, ks].
 
         Reference:
             [1] "Optical Aberrations Correction in Postprocessing Using Imaging Simulation", TOG 2021.
 
         Note:
             Different from the ZEMAX Huygens PSF, which traces rays to the image
             plane and performs plane-wave integration.
+            The complex field has an arbitrary global phase because optical path
+            length is measured relative to a reference; relative phase is meaningful.
         """
         wvln = self.primary_wvln if wvln is None else wvln
         if self.dtype != torch.float64:
@@ -541,6 +610,11 @@ def psf_huygens(
             # Sum contributions from this batch
             psf_complex += complex_amp.sum(dim=-1)  # [ks, ks]
 
+        if return_field:
+            field = torch.flip(psf_complex, [-2, -1])
+            energy = field.abs().square().sum(dim=(-2, -1), keepdim=True)
+            return field / torch.sqrt(energy + EPSILON)
+
         # Convert complex field to intensity
         psf = psf_complex.abs() ** 2
 
 
@@ -178,7 +178,69 @@ def _get_curvature_params(self):
         """
         return self.c, self.k
 
-    def paraxial_power(self, n1, n2):
+    def newtons_method(self, ray):
+        """Find an intersection, rejecting rays entering behind the aperture rim."""
+        t, valid = super().newtons_method(ray)
+        return t, valid & ~self._enters_through_edge(ray)
+
+    @torch.no_grad()
+    def _enters_through_edge(self, ray):
+        """Identify rays entering the aperture or conic domain behind the face.
+
+        The entry point is pulled radially inside the boundary before evaluating
+        sag, so a grazing ray is not sampled past its intersection. Rays starting
+        inside the valid region retain the usual Newton validity checks.
+
+        Args:
+            ray (Ray): Input ray bundle in local surface coordinates.
+
+        Returns:
+            mask (torch.Tensor): Boolean edge-entry mask, shape [...].
+        """
+        o, d = ray.o, ray.d
+        c, k = self._get_curvature_params()
+        aperture_sq = (self.w**2 + self.h**2) / 4 if self.is_square else self.r**2
+        conic_power = (1 + k) * c**2
+        domain_sq = torch.where(
+            conic_power > 0,
+            1 / torch.where(conic_power > 0, conic_power, 1.0),
+            torch.full_like(conic_power, float("inf")),
+        )
+        radius_sq = domain_sq.clamp(max=aperture_sq).to(o)
+        a = d[..., :2].square().sum(-1)
+        b = (o[..., :2] * d[..., :2]).sum(-1)
+        outside = o[..., :2].square().sum(-1) - radius_sq
+        disc = b.square() - a * outside
+        axial = a == 0
+        denominator = torch.where(axial, 1.0, a)
+        root = disc.clamp(min=0).sqrt()
+        t_in = torch.where(axial, -float("inf"), (-b - root) / denominator)
+        t_out = torch.where(axial, float("inf"), (-b + root) / denominator)
+        reaches = torch.where(axial, outside <= 0, disc > 0)
+        starts_inside = outside <= 0
+
+        if self.is_square:
+            # Square apertures use half-side lengths, not a circular substitute.
+            half_size = o.new_tensor([self.w / 2, self.h / 2])
+            parallel = d[..., :2] == 0
+            safe_d = torch.where(parallel, 1.0, d[..., :2])
+            t1 = (-half_size - o[..., :2]) / safe_d
+            t2 = (half_size - o[..., :2]) / safe_d
+            near = torch.where(parallel, -float("inf"), torch.minimum(t1, t2))
+            far = torch.where(parallel, float("inf"), torch.maximum(t1, t2))
+            t_in = torch.maximum(t_in, near.amax(-1))
+            t_out = torch.minimum(t_out, far.amin(-1))
+            inside_square = o[..., :2].abs() <= half_size
+            reaches = reaches & (~parallel | inside_square).all(-1)
+            starts_inside = starts_inside & inside_square.all(-1)
+
+        enters = reaches & ~starts_inside & (t_in > 0) & (t_out > t_in)
+        entry = o + d * torch.where(enters, t_in, 0.0).unsqueeze(-1)
+        xy = torch.where(enters.unsqueeze(-1), entry[..., :2], 0.0) * (1 - 1e-6)
+        behind = (self._sag(xy[..., 0], xy[..., 1]) - entry[..., 2]) * d[..., 2] < 0
+        return enters & behind
+
+    def paraxial_power(self, n1, n2, wvln=None):
         """Return the paraxial optical power of this surface [1/mm].
 
         Only the vertex curvature contributes, so the conic constant and the
@@ -190,6 +252,8 @@ def paraxial_power(self, n1, n2):
         Args:
             n1 (torch.Tensor): Refractive index of the incident medium.
             n2 (torch.Tensor): Refractive index of the transmission medium.
+            wvln (float, optional): Wavelength [µm]. Unused by refractive
+                surfaces, which see it only through `n1`/`n2`.
 
         Returns:
             power (torch.Tensor): Surface power `(n2 - n1) * (c + 2 * a2)`