Merge pull request #41 from HarrisonKramer/feat/aspheres

HarrisonKramer · web-flow · commit 11c5bcc13c71 · 2025-01-26T17:13:12.000+01:00
Feat/aspheres
diff --git a/docs/api/api_geometries.rst b/docs/api/api_geometries.rst
@@ -14,6 +14,7 @@ that require a numerical solution to find the intersection point with a ray.
    geometries.base
    geometries.chebyshev
    geometries.even_asphere
+   geometries.odd_asphere
    geometries.newton_raphson
    geometries.plane
    geometries.polynomial
diff --git a/docs/developers_guide/geometry_overview.rst b/docs/developers_guide/geometry_overview.rst
@@ -50,6 +50,7 @@ Optiland includes a wide range of built-in geometries:
 - **Planes**: Flat surfaces with infinite or finite extent.
 - **Spheres and Conics**: Defined by radius of curvature and conic constant.
 - **Even Aspheres**: Described by polynomial terms for deviations from a sphere.
+- **Odd Aspheres**: Similar to even aspheres but with additional terms for odd powers.
 - **Polynomial XY and Chebyshev Surfaces**: Useful for advanced freeform optical systems.
 - **Custom Geometries**: Users can easily extend the framework by subclassing the `BaseGeometry` (analytical geometries) or `NewtonRaphsonGeometry` (iterative geometries) classes.
 
diff --git a/docs/gallery/specialized_lenses.rst b/docs/gallery/specialized_lenses.rst
@@ -10,6 +10,7 @@ While these lenses are perhaps more complex, they are built using the same basic
     specialized_lenses/objective
     specialized_lenses/f_theta_lens
     specialized_lenses/f_theta_with_fold_mirror
+    specialized_lenses/odd_asphere
     specialized_lenses/microscope
     specialized_lenses/catadioptric
     specialized_lenses/endoscope
diff --git a/docs/gallery/specialized_lenses/odd_asphere.ipynb b/docs/gallery/specialized_lenses/odd_asphere.ipynb
diff --git a/optiland/fileio/converters.py b/optiland/fileio/converters.py
@@ -73,7 +73,7 @@ def _configure_surface_coefficients(self, data):
         surf_type = data['type']
         if surf_type == 'standard':
             return None
-        elif surf_type == 'even_asphere':
+        elif surf_type in ['even_asphere', 'odd_asphere']:
             coefficients = []
             for k in range(8):
                 coefficients.append(data[f'param_{k}'])
diff --git a/optiland/fileio/zemax_handler.py b/optiland/fileio/zemax_handler.py
@@ -418,7 +418,8 @@ def _read_surf_type(self, data):
             data (list): List of data values extracted from the Zemax file.
         """
         type_map = {'STANDARD': 'standard',
-                    'EVENASPH': 'even_asphere'}
+                    'EVENASPH': 'even_asphere',
+                    'ODDASPHE': 'odd_asphere'}
         try:
             self._current_surf_data['type'] = type_map[data[1]]
         except KeyError:
diff --git a/optiland/geometries/__init__.py b/optiland/geometries/__init__.py
@@ -5,5 +5,6 @@
 from .standard import StandardGeometry
 from .newton_raphson import NewtonRaphsonGeometry
 from .even_asphere import EvenAsphere
+from .odd_asphere import OddAsphere
 from .polynomial import PolynomialGeometry
 from .chebyshev import ChebyshevPolynomialGeometry
diff --git a/optiland/geometries/base.py b/optiland/geometries/base.py
@@ -25,6 +25,9 @@ def __init_subclass__(cls, **kwargs):
         super().__init_subclass__(**kwargs)
         BaseGeometry._registry[cls.__name__] = cls
 
+    def __str__(self):
+        return f'{self.__class__.__name__}'  # pragma: no cover
+
     @abstractmethod
     def sag(self, x=0, y=0):
         """Calculate the surface sag of the geometry.
diff --git a/optiland/geometries/chebyshev.py b/optiland/geometries/chebyshev.py
@@ -66,6 +66,9 @@ def __init__(self, coordinate_system, radius, conic=0.0,
         self.norm_y = norm_y
         self.is_symmetric = False
 
+    def __str__(self):
+        return 'Chebyshev Polynomial'
+
     def sag(self, x=0, y=0):
         """
         Calculates the sag of the Chebyshev polynomial surface at the given
diff --git a/optiland/geometries/even_asphere.py b/optiland/geometries/even_asphere.py
@@ -49,6 +49,10 @@ def __init__(self, coordinate_system, radius, conic=0.0,
         super().__init__(coordinate_system, radius, conic, tol, max_iter)
         self.c = coefficients
         self.is_symmetric = True
+        self.order = 2  # used for optimization scaling
+
+    def __str__(self):
+        return 'Even Asphere'
 
     def sag(self, x=0, y=0):
         """
diff --git a/optiland/geometries/newton_raphson.py b/optiland/geometries/newton_raphson.py
@@ -35,6 +35,9 @@ def __init__(self, coordinate_system, radius, conic=0.0, tol=1e-10,
         self.tol = tol
         self.max_iter = max_iter
 
+    def __str__(self):
+        return 'Newton Raphson'  # pragma: no cover
+
     @abstractmethod
     def sag(self, x=0, y=0):
         """Calculate the surface sag of the geometry.
diff --git a/optiland/geometries/odd_asphere.py b/optiland/geometries/odd_asphere.py
@@ -0,0 +1,116 @@
+"""Odd Asphere Geometry
+
+The Even Asphere geometry represents a surface defined by an even asphere in
+two dimensions. The surface is defined as:
+
+z = r^2 / (R * (1 + sqrt(1 - (1 + k) * r^2 / R^2))) + sum(Ci * r^i)
+
+where
+- r^2 = x^2 + y^2
+- R is the radius of curvature
+- k is the conic constant
+- Ci are the aspheric coefficients
+
+Kramer Harrison, 2025
+"""
+import warnings
+import numpy as np
+from optiland.geometries.even_asphere import EvenAsphere
+
+
+class OddAsphere(EvenAsphere):
+    """
+    Represents an odd asphere geometry defined as:
+
+    z = r^2 / (R * (1 + sqrt(1 - (1 + k) * r^2 / R^2))) + sum(Ci * r^i)
+
+    where
+    - r^2 = x^2 + y^2
+    - R is the radius of curvature
+    - k is the conic constant
+    - Ci are the aspheric coefficients
+    - i is an integer from 1 to n
+
+    Args:
+        coordinate_system (str): The coordinate system used for the geometry.
+        radius (float): The radius of curvature of the geometry.
+        conic (float, optional): The conic constant of the geometry.
+            Defaults to 0.0.
+        tol (float, optional): The tolerance value used in calculations.
+            Defaults to 1e-10.
+        max_iter (int, optional): The maximum number of iterations used in
+            calculations. Defaults to 100.
+        coefficients (list, optional): The coefficients of the asphere.
+            Defaults to an empty list, indicating no aspheric coefficients are
+            used.
+    """
+
+    def __init__(self, coordinate_system, radius, conic=0.0,
+                 tol=1e-10, max_iter=100, coefficients=[]):
+        super().__init__(coordinate_system, radius, conic,
+                         tol, max_iter, coefficients)
+        self.order = 1  # used for optimization scaling
+
+    def __str__(self):
+        return 'Odd Asphere'
+
+    def sag(self, x=0, y=0):
+        """
+        Calculates the sag of the asphere at the given coordinates.
+
+        Args:
+            x (float, np.ndarray, optional): The x-coordinate(s).
+                Defaults to 0.
+            y (float, np.ndarray, optional): The y-coordinate(s).
+                Defaults to 0.
+
+        Returns:
+            float: The sag value at the given coordinates.
+        """
+        r2 = x**2 + y**2
+        r = np.sqrt(r2)
+        z = r2 / (self.radius *
+                  (1 + np.sqrt(1 - (1 + self.k) * r2 / self.radius**2)))
+        for i, Ci in enumerate(self.c):
+            z += Ci * r ** (i + 1)
+
+        return z
+
+    def _surface_normal(self, x, y):
+        """
+        Calculates the surface normal of the asphere at the given x and y
+        position.
+
+        Args:
+            x (np.ndarray): The x values to use for calculation.
+            y (np.ndarray): The y values to use for calculation.
+
+        Returns:
+            tuple: The surface normal components (nx, ny, nz).
+        """
+        r2 = x**2 + y**2
+        r = np.sqrt(r2)
+
+        denom = self.radius * np.sqrt(1 - (1 + self.k)*r2 / self.radius**2)
+        dfdx = x / denom
+        dfdy = y / denom
+
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            for i, Ci in enumerate(self.c):
+                x_term = (i + 1) * x * Ci * r**(i - 1)
+                y_term = (i + 1) * y * Ci * r**(i - 1)
+
+                x_term[~np.isfinite(x_term)] = 0
+                y_term[~np.isfinite(y_term)] = 0
+
+                dfdx += x_term
+                dfdy += y_term
+
+        mag = np.sqrt(dfdx**2 + dfdy**2 + 1)
+
+        nx = dfdx / mag
+        ny = dfdy / mag
+        nz = -1 / mag
+
+        return nx, ny, nz
diff --git a/optiland/geometries/plane.py b/optiland/geometries/plane.py
@@ -23,6 +23,9 @@ def __init__(self, coordinate_system):
         self.radius = np.inf
         self.is_symmetric = True
 
+    def __str__(self):
+        return 'Planar'
+
     def sag(self, x=0, y=0):
         """Calculate the surface sag of the plane geometry.
 
diff --git a/optiland/geometries/polynomial.py b/optiland/geometries/polynomial.py
@@ -59,6 +59,9 @@ def __init__(self, coordinate_system, radius, conic=0.0,
         if len(self.c) == 0:
             self.c = np.zeros((1, 1))
 
+    def __str__(self):
+        return 'Polynomial XY'
+
     def sag(self, x=0, y=0):
         """
         Calculates the sag of the polynomial surface at the given coordinates.
diff --git a/optiland/geometries/standard.py b/optiland/geometries/standard.py
@@ -44,6 +44,9 @@ def __init__(self, coordinate_system, radius, conic=0.0):
         self.k = conic
         self.is_symmetric = True
 
+    def __str__(self):
+        return 'Standard'
+
     def sag(self, x=0, y=0):
         """Calculate the surface sag of the geometry at the given coordinates.
 
diff --git a/optiland/optimization/variable/asphere_coeff.py b/optiland/optimization/variable/asphere_coeff.py
@@ -32,6 +32,10 @@ def __init__(self, optic, surface_number, coeff_number, apply_scaling=True,
         super().__init__(optic, surface_number, apply_scaling, **kwargs)
         self.coeff_number = coeff_number
 
+        # Scaling changes with the order of the asphere per coefficient
+        surf = self._surfaces.surfaces[self.surface_number]
+        self.order = surf.geometry.order
+
     def get_value(self):
         """
         Get the current value of the aspheric coefficient.
@@ -64,7 +68,7 @@ def scale(self, value):
         Args:
             value: The value to scale
         """
-        return value * 10 ** (4 + 2 * self.coeff_number)
+        return value * 10 ** (4 + self.order * self.coeff_number)
 
     def inverse_scale(self, scaled_value):
         """
@@ -73,7 +77,7 @@ def inverse_scale(self, scaled_value):
         Args:
             scaled_value: The scaled value to inverse scale
         """
-        return scaled_value / 10 ** (4 + 2 * self.coeff_number)
+        return scaled_value / 10 ** (4 + self.order * self.coeff_number)
 
     def __str__(self):
         """
diff --git a/optiland/surfaces/surface_factory.py b/optiland/surfaces/surface_factory.py
@@ -17,6 +17,7 @@
     Plane,
     StandardGeometry,
     EvenAsphere,
+    OddAsphere,
     PolynomialGeometry,
     ChebyshevPolynomialGeometry)
 from optiland.surfaces.object_surface import ObjectSurface
@@ -82,6 +83,10 @@ def create_surface(self, surface_type, index, is_stop, material,
                 'geometry': self._configure_even_asphere_geometry,
                 'expected_params': ['radius', 'conic', 'coefficients']
             },
+            'odd_asphere': {
+                'geometry': self._configure_odd_asphere_geometry,
+                'expected_params': ['radius', 'conic', 'coefficients']
+            },
             'polynomial': {
                 'geometry': self._configure_polynomial_geometry,
                 'expected_params': ['radius', 'conic', 'coefficients',
@@ -207,6 +212,29 @@ def _configure_even_asphere_geometry(cs, **kwargs):
 
         return geometry
 
+    @staticmethod
+    def _configure_odd_asphere_geometry(cs, **kwargs):
+        """
+        Configures an odd asphere geometry based on the given parameters.
+
+        Parameters:
+            cs: The coordinate system for the geometry.
+            **kwargs: Additional keyword arguments for the geometry. Options
+                include radius, conic, and coefficients.
+
+        Returns:
+            geometry: The configured geometry object.
+        """
+        radius = kwargs.get('radius', np.inf)
+        conic = kwargs.get('conic', 0)
+        tol = kwargs.get('tol', 1e-6)
+        max_iter = kwargs.get('max_iter', 100)
+        coefficients = kwargs.get('coefficients', [])
+
+        geometry = OddAsphere(cs, radius, conic, tol, max_iter, coefficients)
+
+        return geometry
+
     @staticmethod
     def _configure_polynomial_geometry(cs, **kwargs):
         """
diff --git a/optiland/visualization/visualization.py b/optiland/visualization/visualization.py
@@ -13,7 +13,7 @@
 import pandas as pd
 import matplotlib.pyplot as plt
 import vtk
-from optiland import geometries, materials
+from optiland import materials
 from optiland.visualization.rays import Rays2D, Rays3D
 from optiland.visualization.system import OpticalSystem
 
@@ -193,12 +193,11 @@ def view(self):
 
         surf_type = []
         for surf in self.optic.surface_group.surfaces:
-            if isinstance(surf.geometry, geometries.EvenAsphere):
-                surf_type.append('Even Asphere')
-            elif isinstance(surf.geometry, geometries.Plane):
-                surf_type.append('Planar')
-            elif isinstance(surf.geometry, geometries.StandardGeometry):
-                surf_type.append('Standard')
+            g = surf.geometry
+
+            # check if __str__ method exists
+            if type(g).__dict__.get('__str__'):
+                surf_type.append(str(surf.geometry))
             else:
                 raise ValueError('Unknown surface type')
 
diff --git a/tests/test_geometries.py b/tests/test_geometries.py
diff --git a/tests/test_surface_factory.py b/tests/test_surface_factory.py
