docs(example-04): fixed meridional plane reduction and updated text and visualizations

Author: Littie28

docs/examples/index.rst

@@ -9,3 +9,4 @@ Interactive examples demonstrating coordinate frame transformations.
    01-surface_intersections
    02-frame-hierarchies
    03-simple-robot-arm
+   04-simple-raytracing

examples/04-simple-raytracing.ipynb

@@ -14,13 +14,13 @@
     "\n",
     "**Key demonstrations**:\n",
     "1. **Frame transformations simplify optics**: Position/rotate lenses by modifying frames - the raytracing algorithm stays unchanged\n",
-    "2. **Meridional plane reduction**: Transform 3D intersection problem to 2D by exploiting rotational symmetry\n",
-    "3. **Interchangeable profiles**: Swap optical surfaces (spherical, hyperbolic) without changing the tracing code\n",
+    "2. **Interchangeable profiles**: Swap optical surfaces (spherical, hyperbolic) without changing the tracing code\n",
     "\n",
     "**Note**: This is a simplified raytracer for demonstration. Production implementations would handle:\n",
     "  - Rays missing the surface (returns error instead of graceful fallback)\n",
     "  - Rays starting inside the surface\n",
-    "  - Numerical stability near grazing incidence"
+    "  - Numerical stability near grazing incidence\n",
+    "  - concave surfaces"
    ]
   },
   {
@@ -264,12 +264,11 @@
     "        \"\"\"\n",
     "        return self.profile(np.asarray(r))\n",
     "\n",
-    "    def intersect(self, ray: Ray) -> tuple[float, Point, Vector] | None:\n",
+    "    def intersect(self, ray: Ray) -> tuple[float, Point, Vector]:\n",
     "        \"\"\"Find ray-surface intersection.\n",
     "\n",
-    "        Strategy: Transform to meridional plane to reduce 3D problem to 2D.\n",
-    "        The meridional plane contains the optical axis (x) and the ray origin.\n",
-    "        In this plane, we only need to solve for the axial coordinate x.\n",
+    "        Strategy: Transform the ray to the optics local coordinate system (frame).\n",
+    "        Reduce the problem to 2d by exploiting rotational symmetry of the optical element in its frame.\n",
     "\n",
     "        Args:\n",
     "            ray: Incident ray\n",
@@ -280,38 +279,31 @@
     "        \"\"\"\n",
     "        local_ray = ray.to_frame(self.frame)\n",
     "\n",
-    "        # Create meridional plane: contains optical axis and ray origin\n",
-    "        # This reduces the 3D intersection to a 2D problem (x, y_meridional)\n",
-    "        phi = np.arctan2(local_ray.origin.z, local_ray.origin.y)\n",
-    "        meridional = self.frame.make_child(\"meridional\").rotate_euler(x=phi)\n",
-    "        meri_ray = ray.to_frame(meridional)\n",
-    "\n",
-    "        # In meridional plane: only solve for x coordinate\n",
     "        def objective_function(t: float) -> float:\n",
-    "            y = meri_ray.origin.y + t * meri_ray.direction.y\n",
-    "            r = abs(y)\n",
+    "            point = local_ray.origin + t * local_ray.direction\n",
+    "            r = np.sqrt(point.y**2 + point.z**2)\n",
     "\n",
     "            if r > self.aperture_radius:\n",
     "                return 1e10 if t > 0 else -1e10\n",
     "\n",
-    "            x_ray = meri_ray.origin.x + t * meri_ray.direction.x\n",
     "            x_surface = self.profile(np.array([r]))[0]\n",
     "\n",
     "            if np.isnan(x_surface):\n",
     "                return 1e10 if t > 0 else -1e10\n",
     "\n",
-    "            return x_ray - x_surface\n",
+    "            return point.x - x_surface\n",
     "\n",
-    "        # Use bracketing method (robust for 1D problems)\n",
     "        # Search for sign change between t=0 and t=100\n",
     "        try:\n",
     "            t_hit, res = brentq(\n",
-    "                objective_function, 0.1, 100.0, xtol=1e-10, full_output=True\n",
+    "                objective_function, 1e-3, 100.0, xtol=1e-10, full_output=True\n",
     "            )\n",
-    "        except ValueError:\n",
-    "            return None\n",
+    "        except ValueError as err:\n",
+    "            raise ValueError(\n",
+    "                \"No valid intersection between ray and surface found\"\n",
+    "            ) from err\n",
     "\n",
-    "        propagated_ray = meri_ray.propagate(t_hit)\n",
+    "        propagated_ray = local_ray.propagate(t_hit)\n",
     "        hit_local = propagated_ray.origin.to_frame(self.frame)\n",
     "\n",
     "        # Compute surface normal in local frame\n",
@@ -366,7 +358,7 @@
     "        y = np.linspace(-self.aperture_radius, self.aperture_radius, n_points)\n",
     "        x = self.profile(np.abs(y))\n",
     "\n",
-    "        defaults = {\"linewidth\": 2}\n",
+    "        defaults: dict = {\"linewidth\": 2}\n",
     "        defaults.update(kwargs)\n",
     "        ax.plot(x, y, **defaults)\n",
     "\n",
@@ -582,10 +574,10 @@
     "    surface: AxisymmetricSurface,\n",
     "    n_rays: int = 10,\n",
     "    inner_radius: float = 10.0,\n",
-    "    outer_radius: float = 22.0,\n",
+    "    outer_radius: float = 20.0,\n",
     "    n1: float = 1.0,\n",
     "    n2: float = 1.5,\n",
-    "    final_x: float = 33.0,\n",
+    "    final_x: float = 38.0,\n",
     ") -> tuple[dict[int, RayHistory], tuple]:\n",
     "    \"\"\"Trace parallel rays through an optical surface.\n",
     "\n",
@@ -615,9 +607,9 @@
     "    # Generate rays in circular pattern\n",
     "    phis = np.linspace(0, 2 * np.pi, n_rays, endpoint=False)\n",
     "\n",
-    "    for radius, cmap in zip([inner_radius, outer_radius], [plt.cm.Blues, plt.cm.Reds]):\n",
+    "    for radius, cmap in zip([inner_radius, outer_radius], [plt.cm.Reds, plt.cm.Blues]):\n",
     "        rays = defaultdict(RayHistory)\n",
-    "        colors = cmap(np.linspace(0.3, 0.9, n_rays))\n",
+    "        colors = cmap(np.linspace(0.2, 0.8, n_rays))\n",
     "        for i, phi in enumerate(phis):\n",
     "            z = radius * np.cos(phi)\n",
     "            y = radius * np.sin(phi)\n",
@@ -654,14 +646,14 @@
     "                *ray_history.history[0].point[[1, 2]],\n",
     "                color=color,\n",
     "                s=80,\n",
-    "                alpha=0.3,\n",
+    "                alpha=0.5,\n",
     "                edgecolors=\"none\",\n",
     "            )\n",
     "            ax_2d.scatter(\n",
     "                *ray_history.history[-1].point[[1, 2]],\n",
     "                color=color,\n",
     "                s=80,\n",
-    "                edgecolors=\"white\",\n",
+    "                edgecolors=\"black\",\n",
     "                linewidth=1.5,\n",
     "            )\n",
     "\n",
@@ -681,7 +673,7 @@
    "id": "10",
    "metadata": {},
    "source": [
-    "# Spherical Aberration Example\n",
+    "## Example: Raytracing a lens\n",
     "\n",
     "This example demonstrates how coordinate frame transformations simplify optical raytracing.\n",
     "\n",
@@ -695,7 +687,7 @@
    "id": "11",
    "metadata": {},
    "source": [
-    "## Create root frame"
+    "### Create root frame"
    ]
   },
   {
@@ -713,7 +705,7 @@
    "id": "13",
    "metadata": {},
    "source": [
-    "## Aligned spherical lens\n",
+    "### Aligned spherical lens\n",
     "\n",
     "Start with a perfectly aligned spherical lens. Notice the **spherical aberration**: rays at different radial distances don't focus to the same point - the footprint spreads out."
    ]
@@ -749,7 +741,7 @@
    "id": "15",
    "metadata": {},
    "source": [
-    "## Misalign the lens - only one line changes!\n",
+    "### Misalign the lens - only one line changes!\n",
     "\n",
     "Now rotate the lens 5 degrees around the y-axis. **Look closely**: Only the `lens_frame` definition changes - the entire raytracing algorithm (intersection, refraction, propagation) stays identical.\n",
     "\n",
@@ -765,7 +757,7 @@
    "source": [
     "# Misaligned case: lens rotated 5 degrees\n",
     "# ↓↓↓ ONLY THIS LINE CHANGES - everything else is identical! ↓↓↓\n",
-    "lens_frame = root.make_child(\"lens\").translate(x=-20).rotate_euler(y=5, degrees=True)\n",
+    "lens_frame = root.make_child(\"lens\").translate(x=-20).rotate_euler(y=10, degrees=True)\n",
     "\n",
     "profile, derivative = spherical_profile(R=25)\n",
     "\n",
@@ -787,7 +779,7 @@
    "id": "17",
    "metadata": {},
    "source": [
-    "## Swap the surface profile - only the profile changes!\n",
+    "### Swap the surface profile - only the profile changes!\n",
     "\n",
     "Replace the spherical lens with a hyperbolic one. Again, **only the profile function changes** - frame setup and raytracing code stay identical.\n",
     "\n",
@@ -805,7 +797,7 @@
     "# ↓↓↓ ONLY THE PROFILE CHANGES - frame, raytracing, everything else identical! ↓↓↓\n",
     "lens_frame = root.make_child(\"hyperbolic_lens\").translate(x=-25)\n",
     "\n",
-    "profile, derivative = hyperbolic_profile(R=25, conic_constant=-1.0)\n",
+    "profile, derivative = hyperbolic_profile(R=25, conic_constant=-2)\n",
     "\n",
     "surface = AxisymmetricSurface(\n",
     "    frame=lens_frame,\n",
@@ -829,18 +821,11 @@
     "\n",
     "**Frame transformations for optics**: Optical elements (lenses, mirrors) live in their own frames. Repositioning or rotating them requires changing only the frame definition - all geometric calculations (intersection finding, normal computation, refraction) work in local coordinates and remain unchanged.\n",
     "\n",
-    "**Meridional plane reduction**: For rotationally symmetric surfaces, we reduce the 3D ray-surface intersection to a 1D root-finding problem by:\n",
-    "1. Transforming the ray to a meridional plane (contains optical axis + ray origin)\n",
-    "2. Solving for the single parameter `t` (propagation distance) where ray hits surface\n",
-    "This exploits symmetry to simplify the problem from 3 unknowns to 1.\n",
-    "\n",
     "**Interchangeable profiles**: Each surface profile (`spherical`, `hyperbolic`, etc.) is just a function `x = f(r)`. Swapping profiles means changing the function - the `AxisymmetricSurface` class and raytracing algorithm stay identical. This modularity allows easy exploration of different optical designs.\n",
     "\n",
     "**Immutable rays**: Like `Point` and `Vector`, `Ray` is immutable - all operations (`propagate`, `refract`, `to_frame`) return new instances. This matches hazy's philosophy: geometric primitives are values, not mutable state.\n",
     "\n",
-    "**Real raytracing**: This isn't a toy example - the meridional plane technique and conic surface formulas are used in real optical design software. Frame-based coordinates make the code match how optical engineers think about lens systems.\n",
-    "\n",
-    "**Without hazy**: Manually compose transformation matrices for each lens position, track multiplication order for nested frames (root → lens → meridional), explicitly transform every point and vector, rewrite intersection code when repositioning optics.\n",
+    "**Without hazy**: Manually compose transformation matrices for each lens position, track multiplication order for nested frames, explicitly transform every point and vector, rewrite intersection code when repositioning optics.\n",
     "\n",
     "**With hazy**: Define surfaces in natural local coordinates (`origin=(0,0,0)`, x-axis = optical axis), modify frames to position/orient elements, call `.to_frame()` for automatic transformation through arbitrarily deep hierarchies."
    ]