Module traceon.field

Radial series expansion in cylindrical symmetry

Let \phi_0(z) be the potential along the optical axis. We can express the potential around the optical axis as:

\phi = \phi_0(z_0) - \frac{r^2}{4} \frac{\partial \phi_0^2}{\partial z^2} + \frac{r^4}{64} \frac{\partial^4 \phi_0}{\partial z^4} - \frac{r^6}{2304} \frac{\partial \phi_0^6}{\partial z^6} + \cdots

Therefore, if we can efficiently compute the axial potential derivatives \frac{\partial \phi_0^n}{\partial z^n} we can compute the potential and therefore the fields around the optical axis. For the derivatives of \phi_0(z) closed form formulas exist in the case of radially symmetric geometries, see for example formula 13.16a in [1]. Traceon uses a recursive version of these formulas to very efficiently compute the axial derivatives of the potential.

[1] P. Hawkes, E. Kasper. Principles of Electron Optics. Volume one: Basic Geometrical Optics. 2018.

Classes

class Field
Expand source code
class Field(GeometricObject, ABC):
    def __init__(self) -> None:
        self._origin = np.array([0,0,0], dtype=np.float64)
        self._basis = np.eye(3, dtype=np.float64)
        self._update_inverse_transformation_matrix()

        self.field_bounds: Bounds3D | None = None

    def get_origin(self) -> Point3D:
        """
        Get the origin of the field in the global coordinate system. This is the position
        that the origin (0, 0, 0) was transformed to by using methods from `traceon.mesher.GeometricObject`.

        Returns
        -----------------------------
        numpy.ndarray
            Float array of shape (3,)
        """
        return self._origin.copy()
    
    def get_basis(self) -> ArrayFloat2D:
        return self._basis.copy()

    def _update_inverse_transformation_matrix(self) -> None:
        transformation_matrix = np.eye(4)
        transformation_matrix[:3, :3] = self._basis
        transformation_matrix[:3, 3] = self._origin

        assert np.linalg.det(transformation_matrix) != 0, ("Transformations of field have resulted in a two-dimensional coordinate system. "
                                                           "Please only use affine transformations.")
        self._inverse_transformation_matrix = np.linalg.inv(transformation_matrix)

    def copy(self) -> Self:
        return copy.copy(self)
    
    def map_points(self, fun: Callable[[PointLike3D], Point3D]) -> Self:
        field_copy = self.copy()
        
        field_copy._origin = fun(self._origin).astype(np.float64)
        assert field_copy._origin.shape == (3,), "Transformation of field did not map origin to a 3D point"
        
        field_copy._basis = np.array([fun(b + self._origin) - field_copy._origin for b in self._basis])
        assert field_copy._basis.shape == (3,3), "Transformation of field did not map unit vectors to a 3D vector"
        
        field_copy._update_inverse_transformation_matrix()
        return field_copy

    def map_points_to_local(self, point: PointLike3D) -> Point3D:
        """Converts a point from the global coordinate system to the local coordinate system of the field. 
    
        Parameters
        ---------------------
        point: (3,) np.ndarray of float64
            The coordinates of the point in the global coordinate system.

        Returns
        ---------------------
        (3,) np.ndarray of float64
            The coordinates of the point in the local coordinate system."""
        # represent the point in homogenous coordinates so we can do the inverse 
        # affine transformation with a single matrix multiplication.
        global_point_homogeneous = np.array([*point, 1.], dtype=np.float64)
        local_point_homogeneous = self._inverse_transformation_matrix @ global_point_homogeneous
        assert np.isclose(local_point_homogeneous[3], 1.)
        return local_point_homogeneous[:3]

    def set_bounds(self, bounds: BoundsLike3D, global_coordinates: bool = False) -> None:
        """Set the field bounds. Outside the field bounds the field always returns zero (i.e. no field). Note
        that even in 2D the field bounds needs to be specified for x,y and z axis. The trajectories in the presence
        of magnetostatic field are in general 3D even in radial symmetric geometries.
        
        Parameters
        -------------------
        bounds: (3, 2) np.ndarray of float64
            The min, max value of x, y, z respectively within the field is still computed.
        global_coordinates: bool
            If `True` the given bounds are in global coordinates and transformed to the fields local system internally.
        """
        bounds = np.array(bounds, dtype=np.float64)
        assert bounds.shape == (3,2)

        if global_coordinates:
            transformed_corners = np.array([self.map_points_to_local(corner) for corner in product(*bounds)])
            bounds = np.column_stack((transformed_corners.min(axis=0), transformed_corners.max(axis=0)))

        self.field_bounds = bounds
    
    def _within_field_bounds(self, point: PointLike3D) -> bool:
        return bool(self.field_bounds is None or np.all((self.field_bounds[:, 0] <= point) & (point <= self.field_bounds[:, 1])))

    def _matches_geometry(self, other: Field) -> bool:
        return False
    
    def field_at_point(self, point: PointLike3D) -> Vector3D:
        """Convenience function for getting the field in the case that the field is purely electrostatic
        or magneotstatic. Automatically picks one of `electrostatic_field_at_point` or `magnetostatic_field_at_point`.
        Throws an exception when the field is both electrostatic and magnetostatic.

        Parameters
        ---------------------
        point: (3,) np.ndarray of float64

        Returns
        --------------------
        (3,) np.ndarray of float64. The electrostatic field \\(\\vec{E}\\) or the magnetostatic field \\(\\vec{H}\\).
        """
        elec, mag = self.is_electrostatic(), self.is_magnetostatic()
        
        if elec and not mag:
            return self.electrostatic_field_at_point(point)
        elif not elec and mag:
            return self.magnetostatic_field_at_point(point)
         
        raise RuntimeError("Cannot use field_at_point when both electric and magnetic fields are present, " \
            "use electrostatic_field_at_point or magnetostatic_field_at_point")
     
    def potential_at_point(self, point: Point3D) -> float:
        """Convenience function for getting the potential in the case that the field is purely electrostatic
        or magneotstatic. Automatically picks one of `electrostatic_potential_at_point` or `magnetostatic_potential_at_point`.
        Throws an exception when the field is both electrostatic and magnetostatic.
         
        Parameters
        ---------------------
        point: (3,) np.ndarray of float64

        Returns
        --------------------
        float. The electrostatic potential (unit Volt) or magnetostaic scalar potential (unit Ampere)
        """
        elec, mag = self.is_electrostatic(), self.is_magnetostatic()
         
        if elec and not mag:
            return self.electrostatic_potential_at_point(point)
        elif not elec and mag:
            return self.magnetostatic_potential_at_point(point)
         
        raise RuntimeError("Cannot use potential_at_point when both electric and magnetic fields are present, " \
            "use electrostatic_potential_at_point or magnetostatic_potential_at_point")

    def electrostatic_field_at_point(self, point: PointLike3D) -> Vector3D:
        """
        Compute the electric field, \\( \\vec{E} = -\\nabla \\phi \\)
        
        Parameters
        ----------
        point: (3,) array of float64
            Position in global coordinate system at which to compute the field.
        
        Returns
        -------
        (3,) array of float64, containing the field strengths (units of V/m)
        """
        local_point = self.map_points_to_local(point)

        if self._within_field_bounds(local_point):
            return self._basis @ self.electrostatic_field_at_local_point(local_point)
        else:
             return np.array([0.,0.,0.])
        
    def magnetostatic_field_at_point(self, point: PointLike3D) -> Vector3D:
        """
        Compute the magnetic field \\( \\vec{H} \\)
        
        Parameters
        ----------
        point: (3,) array of float64
            Position in global coordinate system at which to compute the field.
             
        Returns
        -------
        (3,) np.ndarray of float64 containing the field strength (in units of A/m) in the x, y and z directions.
        """
        local_point = self.map_points_to_local(point)
        if self._within_field_bounds(local_point):
            return self._basis @ self.magnetostatic_field_at_local_point(local_point)
        else:
             return np.array([0.,0.,0.])
    
    def current_field_at_point(self, point: PointLike3D) -> Vector3D:
        """
        Compute the magnetic field produced by currents \\( \\vec{H} \\)
        
        Parameters
        ----------
        point: (3,) array of float64
            Position in global coordinate system at which to compute the field.
             
        Returns
        -------
        (3,) np.ndarray of float64 containing the field strength (in units of A/m) in the x, y and z directions.
        """
        local_point = self.map_points_to_local(point)
        if self._within_field_bounds(local_point):
            return self._basis @ self.current_field_at_local_point(local_point)
        else:
             return np.array([0.,0.,0.])

    
    def electrostatic_potential_at_point(self, point: PointLike3D) -> float:
        """
        Compute the electrostatic potential.
        
        Parameters
        ----------
        point: (3,) array of float64
            Position in global coordinate system at which to compute the field.
        
        Returns
        -------
        Potential as a float value (in units of V).
        """
        local_point = self.map_points_to_local(point)
        
        if self._within_field_bounds(local_point):
            return self.electrostatic_potential_at_local_point(local_point)
        else:
             return 0.

    def magnetostatic_potential_at_point(self, point: PointLike3D) -> float:
        """
        Compute the magnetostatic scalar potential (satisfying \\(\\vec{H} = -\\nabla \\phi \\))
        
        Parameters
        ----------
        point: (3,) array of float64
            Position in global coordinate system in local coordinate system at which to compute the field.
        
        Returns
        -------
        Potential as a float value (in units of A).
        """
        local_point = self.map_points_to_local(point)
        if self._within_field_bounds(local_point):
            return self.magnetostatic_potential_at_local_point(local_point)
        else:
             return 0.

    def is_electrostatic(self) -> bool:
        return False
     
    def is_magnetostatic(self) -> bool:
        return False
     
    def electrostatic_field_at_local_point(self, point) -> Vector3D:
        return np.zeros(3)
    
    def magnetostatic_field_at_local_point(self, point) -> Vector3D:
        return np.zeros(3)
     
    def current_field_at_local_point(self, point) -> Vector3D:
        return np.zeros(3)
    
    def electrostatic_potential_at_local_point(self, point) -> float:
        return 0.0
    
    def magnetostatic_potential_at_local_point(self, point) -> float:
        return 0.0
     
    def __add__(self, other: Field) -> Field:
        if isinstance(other, Field) and not isinstance(other, FieldSuperposition):
            return FieldSuperposition([self, other])

        return NotImplemented
    
    def __mul__(self, other: float) -> Field:
        if _is_numeric(other):
            return FieldSuperposition([self], [other])
        
        return NotImplemented
    
    def __rmul__(self, other: float) -> Field:
        return self.__mul__(other)
    
    def __neg__(self) -> Field:
        return -1*self
    
    def __sub__(self, other: Field) -> Field:
        if isinstance(other, Field):
            return self + (-other)
        
        return NotImplemented
    
    # Following function can be implemented to get a speedup while tracing. 
    # Return a field function implemented in C and a ctypes argument needed. 
    # See the field_fun variable in backend/__init__.py.
    # Note that by default it gives back a Python function, which gives no speedup.
    def get_low_level_trace_function(self) -> tuple[Callable, Any] | tuple[Callable, Any, list[Any]]:
        fun = lambda pos, vel: (self.electrostatic_field_at_point(pos), self.magnetostatic_field_at_point(pos))
        return backend.wrap_field_fun(fun), None

The Mesh class (and the classes defined in traceon.geometry) are subclasses of GeometricObject. This means that they all can be moved, rotated, mirrored.

Ancestors

Subclasses

Methods

def current_field_at_local_point(self, point)
Expand source code
def current_field_at_local_point(self, point) -> Vector3D:
    return np.zeros(3)
def current_field_at_point(self, point)
Expand source code
def current_field_at_point(self, point: PointLike3D) -> Vector3D:
    """
    Compute the magnetic field produced by currents \\( \\vec{H} \\)
    
    Parameters
    ----------
    point: (3,) array of float64
        Position in global coordinate system at which to compute the field.
         
    Returns
    -------
    (3,) np.ndarray of float64 containing the field strength (in units of A/m) in the x, y and z directions.
    """
    local_point = self.map_points_to_local(point)
    if self._within_field_bounds(local_point):
        return self._basis @ self.current_field_at_local_point(local_point)
    else:
         return np.array([0.,0.,0.])

Compute the magnetic field produced by currents \vec{H}

Parameters

point : (3,) array of float64
Position in global coordinate system at which to compute the field.

Returns

(3,) np.ndarray of float64 containing the field strength (in units of A/m) in the x, y and z directions.

def electrostatic_field_at_local_point(self, point)
Expand source code
def electrostatic_field_at_local_point(self, point) -> Vector3D:
    return np.zeros(3)
def electrostatic_field_at_point(self, point)
Expand source code
def electrostatic_field_at_point(self, point: PointLike3D) -> Vector3D:
    """
    Compute the electric field, \\( \\vec{E} = -\\nabla \\phi \\)
    
    Parameters
    ----------
    point: (3,) array of float64
        Position in global coordinate system at which to compute the field.
    
    Returns
    -------
    (3,) array of float64, containing the field strengths (units of V/m)
    """
    local_point = self.map_points_to_local(point)

    if self._within_field_bounds(local_point):
        return self._basis @ self.electrostatic_field_at_local_point(local_point)
    else:
         return np.array([0.,0.,0.])

Compute the electric field, \vec{E} = -\nabla \phi

Parameters

point : (3,) array of float64
Position in global coordinate system at which to compute the field.

Returns

(3,) array of float64, containing the field strengths (units of V/m)

def electrostatic_potential_at_local_point(self, point)
Expand source code
def electrostatic_potential_at_local_point(self, point) -> float:
    return 0.0
def electrostatic_potential_at_point(self, point)
Expand source code
def electrostatic_potential_at_point(self, point: PointLike3D) -> float:
    """
    Compute the electrostatic potential.
    
    Parameters
    ----------
    point: (3,) array of float64
        Position in global coordinate system at which to compute the field.
    
    Returns
    -------
    Potential as a float value (in units of V).
    """
    local_point = self.map_points_to_local(point)
    
    if self._within_field_bounds(local_point):
        return self.electrostatic_potential_at_local_point(local_point)
    else:
         return 0.

Compute the electrostatic potential.

Parameters

point : (3,) array of float64
Position in global coordinate system at which to compute the field.

Returns

Potential as a float value (in units of V).

def field_at_point(self, point)
Expand source code
def field_at_point(self, point: PointLike3D) -> Vector3D:
    """Convenience function for getting the field in the case that the field is purely electrostatic
    or magneotstatic. Automatically picks one of `electrostatic_field_at_point` or `magnetostatic_field_at_point`.
    Throws an exception when the field is both electrostatic and magnetostatic.

    Parameters
    ---------------------
    point: (3,) np.ndarray of float64

    Returns
    --------------------
    (3,) np.ndarray of float64. The electrostatic field \\(\\vec{E}\\) or the magnetostatic field \\(\\vec{H}\\).
    """
    elec, mag = self.is_electrostatic(), self.is_magnetostatic()
    
    if elec and not mag:
        return self.electrostatic_field_at_point(point)
    elif not elec and mag:
        return self.magnetostatic_field_at_point(point)
     
    raise RuntimeError("Cannot use field_at_point when both electric and magnetic fields are present, " \
        "use electrostatic_field_at_point or magnetostatic_field_at_point")

Convenience function for getting the field in the case that the field is purely electrostatic or magneotstatic. Automatically picks one of electrostatic_field_at_point or magnetostatic_field_at_point. Throws an exception when the field is both electrostatic and magnetostatic.

Parameters

point : (3,) np.ndarray of float64
 

Returns

(3,) np.ndarray of float64. The electrostatic field \vec{E} or the magnetostatic field \vec{H}.

def get_basis(self)
Expand source code
def get_basis(self) -> ArrayFloat2D:
    return self._basis.copy()
def get_origin(self)
Expand source code
def get_origin(self) -> Point3D:
    """
    Get the origin of the field in the global coordinate system. This is the position
    that the origin (0, 0, 0) was transformed to by using methods from `traceon.mesher.GeometricObject`.

    Returns
    -----------------------------
    numpy.ndarray
        Float array of shape (3,)
    """
    return self._origin.copy()

Get the origin of the field in the global coordinate system. This is the position that the origin (0, 0, 0) was transformed to by using methods from GeometricObject.

Returns

numpy.ndarray
Float array of shape (3,)
def is_electrostatic(self)
Expand source code
def is_electrostatic(self) -> bool:
    return False
def is_magnetostatic(self)
Expand source code
def is_magnetostatic(self) -> bool:
    return False
def magnetostatic_field_at_local_point(self, point)
Expand source code
def magnetostatic_field_at_local_point(self, point) -> Vector3D:
    return np.zeros(3)
def magnetostatic_field_at_point(self, point)
Expand source code
def magnetostatic_field_at_point(self, point: PointLike3D) -> Vector3D:
    """
    Compute the magnetic field \\( \\vec{H} \\)
    
    Parameters
    ----------
    point: (3,) array of float64
        Position in global coordinate system at which to compute the field.
         
    Returns
    -------
    (3,) np.ndarray of float64 containing the field strength (in units of A/m) in the x, y and z directions.
    """
    local_point = self.map_points_to_local(point)
    if self._within_field_bounds(local_point):
        return self._basis @ self.magnetostatic_field_at_local_point(local_point)
    else:
         return np.array([0.,0.,0.])

Compute the magnetic field \vec{H}

Parameters

point : (3,) array of float64
Position in global coordinate system at which to compute the field.

Returns

(3,) np.ndarray of float64 containing the field strength (in units of A/m) in the x, y and z directions.

def magnetostatic_potential_at_local_point(self, point)
Expand source code
def magnetostatic_potential_at_local_point(self, point) -> float:
    return 0.0
def magnetostatic_potential_at_point(self, point)
Expand source code
def magnetostatic_potential_at_point(self, point: PointLike3D) -> float:
    """
    Compute the magnetostatic scalar potential (satisfying \\(\\vec{H} = -\\nabla \\phi \\))
    
    Parameters
    ----------
    point: (3,) array of float64
        Position in global coordinate system in local coordinate system at which to compute the field.
    
    Returns
    -------
    Potential as a float value (in units of A).
    """
    local_point = self.map_points_to_local(point)
    if self._within_field_bounds(local_point):
        return self.magnetostatic_potential_at_local_point(local_point)
    else:
         return 0.

Compute the magnetostatic scalar potential (satisfying \vec{H} = -\nabla \phi )

Parameters

point : (3,) array of float64
Position in global coordinate system in local coordinate system at which to compute the field.

Returns

Potential as a float value (in units of A).

def map_points_to_local(self, point)
Expand source code
def map_points_to_local(self, point: PointLike3D) -> Point3D:
    """Converts a point from the global coordinate system to the local coordinate system of the field. 

    Parameters
    ---------------------
    point: (3,) np.ndarray of float64
        The coordinates of the point in the global coordinate system.

    Returns
    ---------------------
    (3,) np.ndarray of float64
        The coordinates of the point in the local coordinate system."""
    # represent the point in homogenous coordinates so we can do the inverse 
    # affine transformation with a single matrix multiplication.
    global_point_homogeneous = np.array([*point, 1.], dtype=np.float64)
    local_point_homogeneous = self._inverse_transformation_matrix @ global_point_homogeneous
    assert np.isclose(local_point_homogeneous[3], 1.)
    return local_point_homogeneous[:3]

Converts a point from the global coordinate system to the local coordinate system of the field.

Parameters

point : (3,) np.ndarray of float64
The coordinates of the point in the global coordinate system.

Returns

(3,) np.ndarray of float64 The coordinates of the point in the local coordinate system.

def potential_at_point(self, point)
Expand source code
def potential_at_point(self, point: Point3D) -> float:
    """Convenience function for getting the potential in the case that the field is purely electrostatic
    or magneotstatic. Automatically picks one of `electrostatic_potential_at_point` or `magnetostatic_potential_at_point`.
    Throws an exception when the field is both electrostatic and magnetostatic.
     
    Parameters
    ---------------------
    point: (3,) np.ndarray of float64

    Returns
    --------------------
    float. The electrostatic potential (unit Volt) or magnetostaic scalar potential (unit Ampere)
    """
    elec, mag = self.is_electrostatic(), self.is_magnetostatic()
     
    if elec and not mag:
        return self.electrostatic_potential_at_point(point)
    elif not elec and mag:
        return self.magnetostatic_potential_at_point(point)
     
    raise RuntimeError("Cannot use potential_at_point when both electric and magnetic fields are present, " \
        "use electrostatic_potential_at_point or magnetostatic_potential_at_point")

Convenience function for getting the potential in the case that the field is purely electrostatic or magneotstatic. Automatically picks one of electrostatic_potential_at_point or magnetostatic_potential_at_point. Throws an exception when the field is both electrostatic and magnetostatic.

Parameters

point : (3,) np.ndarray of float64
 

Returns

float. The electrostatic potential (unit Volt) or magnetostaic scalar potential (unit Ampere)
 
def set_bounds(self, bounds, global_coordinates=False)
Expand source code
def set_bounds(self, bounds: BoundsLike3D, global_coordinates: bool = False) -> None:
    """Set the field bounds. Outside the field bounds the field always returns zero (i.e. no field). Note
    that even in 2D the field bounds needs to be specified for x,y and z axis. The trajectories in the presence
    of magnetostatic field are in general 3D even in radial symmetric geometries.
    
    Parameters
    -------------------
    bounds: (3, 2) np.ndarray of float64
        The min, max value of x, y, z respectively within the field is still computed.
    global_coordinates: bool
        If `True` the given bounds are in global coordinates and transformed to the fields local system internally.
    """
    bounds = np.array(bounds, dtype=np.float64)
    assert bounds.shape == (3,2)

    if global_coordinates:
        transformed_corners = np.array([self.map_points_to_local(corner) for corner in product(*bounds)])
        bounds = np.column_stack((transformed_corners.min(axis=0), transformed_corners.max(axis=0)))

    self.field_bounds = bounds

Set the field bounds. Outside the field bounds the field always returns zero (i.e. no field). Note that even in 2D the field bounds needs to be specified for x,y and z axis. The trajectories in the presence of magnetostatic field are in general 3D even in radial symmetric geometries.

Parameters

bounds : (3, 2) np.ndarray of float64
The min, max value of x, y, z respectively within the field is still computed.
global_coordinates : bool
If True the given bounds are in global coordinates and transformed to the fields local system internally.

Inherited members

class FieldAxial (field, z, electrostatic_coeffs=None, magnetostatic_coeffs=None)
Expand source code
class FieldAxial(Field, ABC):
    """An electrostatic field resulting from a radial series expansion around the optical axis. You should
    not initialize this class yourself, but it is used as a base class for the fields returned by the `axial_derivative_interpolation` methods. 
    This base class overloads the +,*,- operators so it is very easy to take a superposition of different fields."""
    
    def __init__(self, 
                 field: FieldBEM, 
                 z: ArrayFloat1D, 
                 electrostatic_coeffs: ArrayFloat3D | None = None, 
                 magnetostatic_coeffs: ArrayFloat3D | None = None):
        
        super().__init__()
        self.field = field
        self._origin = field._origin
        self._basis = field._basis
        self._update_inverse_transformation_matrix()

        N = len(z)
        assert z.shape == (N,)
        assert electrostatic_coeffs is None or len(electrostatic_coeffs)== N-1
        assert magnetostatic_coeffs is None or len(magnetostatic_coeffs) == N-1
        assert electrostatic_coeffs is not None or magnetostatic_coeffs is not None
        
        assert z[0] < z[-1], "z values in axial interpolation should be ascending"

        self.z = z
        self.electrostatic_coeffs = electrostatic_coeffs if electrostatic_coeffs is not None else np.zeros_like(magnetostatic_coeffs)
        self.magnetostatic_coeffs = magnetostatic_coeffs if magnetostatic_coeffs is not None else np.zeros_like(electrostatic_coeffs)
        
        self.has_electrostatic = bool(np.any(self.electrostatic_coeffs != 0.))
        self.has_magnetostatic = bool(np.any(self.magnetostatic_coeffs != 0.))
     
    def is_electrostatic(self) -> bool:
        return self.has_electrostatic

    def is_magnetostatic(self) -> bool:
        return self.has_magnetostatic
    
    def _matches_geometry(self, other: Field) -> bool:
        if (self.__class__ != other.__class__):
            return False
        else:
            other = cast(FieldAxial, other)
            return(np.allclose(self._origin, other._origin) 
                and np.allclose(self._basis, other._basis)
                and self.z.shape == other.z.shape and np.allclose(self.z, other.z))
    
    def __str__(self) -> str:
        name = self.__class__.__name__
        return f'<Traceon {name}, zmin={self.z[0]} mm, zmax={self.z[-1]} mm,\n\tNumber of samples on optical axis: {len(self.z)}>'
     
    def __add__(self, other: Field) -> Field:
        if self._matches_geometry(other):
            other = cast(FieldAxial, other)
            field_copy = self.copy()
            field_copy.electrostatic_coeffs = self.electrostatic_coeffs + other.electrostatic_coeffs
            field_copy.magnetostatic_coeffs = self.magnetostatic_coeffs + other.magnetostatic_coeffs
            return field_copy
        else:
            return super().__add__(other)
    
    def __sub__(self, other: Field) -> Field:
        if isinstance(other, Field):
            return self.__add__(-other)

        return NotImplemented

    def __radd__(self, other: Field) -> Field:
        return self.__add__(other)
     
    def __mul__(self, other: float) -> Field:
        if _is_numeric(other):
            field_copy = self.copy()
            field_copy.electrostatic_coeffs = other * self.electrostatic_coeffs
            field_copy.magnetostatic_coeffs = other * self.electrostatic_coeffs
            return field_copy
        else:
            return super().__mul__(other)
     
    def __neg__(self) -> Field:
        return -1*self
    
    def __rmul__(self, other: float) -> Field:
        return self.__mul__(other)

An electrostatic field resulting from a radial series expansion around the optical axis. You should not initialize this class yourself, but it is used as a base class for the fields returned by the axial_derivative_interpolation methods. This base class overloads the +,*,- operators so it is very easy to take a superposition of different fields.

Ancestors

Subclasses

Methods

def is_electrostatic(self)
Expand source code
def is_electrostatic(self) -> bool:
    return self.has_electrostatic
def is_magnetostatic(self)
Expand source code
def is_magnetostatic(self) -> bool:
    return self.has_magnetostatic

Inherited members

class FieldBEM (electrostatic_point_charges, magnetostatic_point_charges, current_point_charges)
Expand source code
class FieldBEM(Field, ABC):
    """An electrostatic field (resulting from surface charges) as computed from the Boundary Element Method. You should
    not initialize this class yourself, but it is used as a base class for the fields returned by the `solve_direct` function. 
    This base class overloads the +,*,- operators so it is very easy to take a superposition of different fields."""
    
    def __init__(self, 
                 electrostatic_point_charges: EffectivePointCharges, 
                 magnetostatic_point_charges: EffectivePointCharges, 
                 current_point_charges: EffectivePointCharges):
        
        super().__init__()
        
        self.electrostatic_point_charges = electrostatic_point_charges
        self.magnetostatic_point_charges = magnetostatic_point_charges
        self.current_point_charges = current_point_charges
        
    def is_electrostatic(self) -> bool:
        return len(self.electrostatic_point_charges) > 0

    def is_magnetostatic(self) -> bool:
        return len(self.magnetostatic_point_charges) > 0 or len(self.current_point_charges) > 0 
    
    def _matches_geometry(self, other: Field) -> bool:
        return (self.__class__ == other.__class__
            and np.allclose(self._origin, other._origin) 
            and np.allclose(self._basis, other._basis))


    def __add__(self, other: Field) -> Field:
        if self._matches_geometry(other):
            other = cast(FieldBEM, other)
            field_copy = self.copy()
            field_copy.electrostatic_point_charges = self.electrostatic_point_charges + other.electrostatic_point_charges
            field_copy.magnetostatic_point_charges = self.magnetostatic_point_charges + other.magnetostatic_point_charges
            field_copy.current_point_charges = self.current_point_charges + other.current_point_charges
            return field_copy
        else:
            return super().__add__(other)
    
    def __sub__(self, other: Field) -> Field:
        if isinstance(other, Field):
            return self.__add__(-other)
        
        return NotImplemented

    def __radd__(self, other: Field) -> Field:
        return self.__add__(other)
        
    def __mul__(self, other: float) -> Field:
        if _is_numeric(other):
           field_copy = self.copy()
           field_copy.electrostatic_point_charges = self.electrostatic_point_charges * other
           field_copy.magnetostatic_point_charges = self.magnetostatic_point_charges * other
           field_copy.current_point_charges = self.current_point_charges * other
           return field_copy
        else:
            return super().__mul__(other)
    
    def __neg__(self) -> FieldBEM:
        return self.__class__(
            self.electrostatic_point_charges.__neg__(),
            self.magnetostatic_point_charges.__neg__(),
            self.current_point_charges.__neg__())
     
    def __rmul__(self, other: float) -> Field:
        return self.__mul__(other)
      
    def area_of_elements(self, indices: ArrayLikeInt1D):
        """Compute the total area of the elements at the given indices.
        
        Parameters
        ------------
        indices: int iterable
            Indices giving which elements to include in the area calculation.

        Returns
        ---------------
        The sum of the area of all elements with the given indices.
        """
        return sum(self.area_of_element(i) for i in indices) 

    @abstractmethod
    def area_of_element(self, i: int) -> float:
        ...
    
    def charge_on_element(self, i: int) -> float:
        return self.area_of_element(i) * self.electrostatic_point_charges.charges[i]
    
    def charge_on_elements(self, indices: ArrayLikeInt1D) -> float:
        """Compute the sum of the charges present on the elements with the given indices. To
        get the total charge of a physical group use `names['name']` for indices where `names` 
        is returned by `traceon.excitation.Excitation.get_electrostatic_active_elements()`.

        Parameters
        ----------
        indices: (N,) array of int
            indices of the elements contributing to the charge sum. 
         
        Returns
        -------
        The sum of the charge. See the note about units on the front page."""
        return sum(self.charge_on_element(i) for i in indices)
    
    def __str__(self) -> str:
        name = self.__class__.__name__
        return f'<Traceon {name}\n' \
            f'\tNumber of electrostatic points: {len(self.electrostatic_point_charges)}\n' \
            f'\tNumber of magnetizable points: {len(self.magnetostatic_point_charges)}\n' \
            f'\tNumber of current rings: {len(self.current_point_charges)}>'

An electrostatic field (resulting from surface charges) as computed from the Boundary Element Method. You should not initialize this class yourself, but it is used as a base class for the fields returned by the solve_direct function. This base class overloads the +,*,- operators so it is very easy to take a superposition of different fields.

Ancestors

Subclasses

Methods

def area_of_element(self, i)
Expand source code
@abstractmethod
def area_of_element(self, i: int) -> float:
    ...
def area_of_elements(self, indices)
Expand source code
def area_of_elements(self, indices: ArrayLikeInt1D):
    """Compute the total area of the elements at the given indices.
    
    Parameters
    ------------
    indices: int iterable
        Indices giving which elements to include in the area calculation.

    Returns
    ---------------
    The sum of the area of all elements with the given indices.
    """
    return sum(self.area_of_element(i) for i in indices) 

Compute the total area of the elements at the given indices.

Parameters

indices : int iterable
Indices giving which elements to include in the area calculation.

Returns

The sum of the area of all elements with the given indices.

def charge_on_element(self, i)
Expand source code
def charge_on_element(self, i: int) -> float:
    return self.area_of_element(i) * self.electrostatic_point_charges.charges[i]
def charge_on_elements(self, indices)
Expand source code
def charge_on_elements(self, indices: ArrayLikeInt1D) -> float:
    """Compute the sum of the charges present on the elements with the given indices. To
    get the total charge of a physical group use `names['name']` for indices where `names` 
    is returned by `traceon.excitation.Excitation.get_electrostatic_active_elements()`.

    Parameters
    ----------
    indices: (N,) array of int
        indices of the elements contributing to the charge sum. 
     
    Returns
    -------
    The sum of the charge. See the note about units on the front page."""
    return sum(self.charge_on_element(i) for i in indices)

Compute the sum of the charges present on the elements with the given indices. To get the total charge of a physical group use names['name'] for indices where names is returned by Excitation.get_electrostatic_active_elements().

Parameters

indices : (N,) array of int
indices of the elements contributing to the charge sum.

Returns

The sum of the charge. See the note about units on the front page.

def is_electrostatic(self)
Expand source code
def is_electrostatic(self) -> bool:
    return len(self.electrostatic_point_charges) > 0
def is_magnetostatic(self)
Expand source code
def is_magnetostatic(self) -> bool:
    return len(self.magnetostatic_point_charges) > 0 or len(self.current_point_charges) > 0 

Inherited members

class FieldRadialAxial (field, zmin, zmax, N=None)
Expand source code
class FieldRadialAxial(FieldAxial):
    def __init__(self, field: FieldRadialBEM, zmin: float, zmax: float, N: int | None = None) -> None:
        """
        Produces a field which uses an axial interpolation to very quickly compute the field around the z-axis.
        Note that the approximation degrades as the point at which the field is computed is further from the z-axis.
        Also note that the fields produced by current and the magnetizable material are merged into one interpolation,
        so the method `current_field_at_point` will always return zero.

        Parameters
        -----------------------
        field: `traceon.field.FieldRadialBEM`
            Field for which to compute the axial interpolation
        zmin : float
            Location on the optical axis where to start sampling the radial expansion coefficients.
        zmax : float
            Location on the optical axis where to stop sampling the radial expansion coefficients. Any field
            evaluation outside [zmin, zmax] will return a zero field strength.
        N: int, optional
            Number of samples to take on the optical axis, if N=None the amount of samples
            is determined by taking into account the number of elements in the mesh.
        """
        assert isinstance(field, FieldRadialBEM)

        z, electrostatic_coeffs, magnetostatic_coeffs = FieldRadialAxial._get_interpolation_coefficients(field, zmin, zmax, N=N)
        
        super().__init__(field, z, electrostatic_coeffs, magnetostatic_coeffs)

        assert self.electrostatic_coeffs.shape == (len(z)-1, backend.DERIV_2D_MAX, 6)
        assert self.magnetostatic_coeffs.shape == (len(z)-1, backend.DERIV_2D_MAX, 6)
    
    @staticmethod
    def _get_interpolation_coefficients(field: FieldRadialBEM, 
                                        zmin: float, 
                                        zmax: float, 
                                        N: int | None = None) -> tuple[ArrayFloat1D, ArrayFloat3D, ArrayFloat3D]:
        
        assert zmax > zmin, "zmax should be bigger than zmin"

        N_charges = max(len(field.electrostatic_point_charges.charges), len(field.magnetostatic_point_charges.charges))
        N = N if N is not None else int(FACTOR_AXIAL_DERIV_SAMPLING_2D*N_charges)
        z = np.linspace(zmin, zmax, N)
        
        st = time.time()
        elec_derivs = np.concatenate(util.split_collect(field.get_electrostatic_axial_potential_derivatives, z), axis=0)
        elec_coeffs = _quintic_spline_coefficients(z, elec_derivs.T)
        
        mag_derivs = np.concatenate(util.split_collect(field.get_magnetostatic_axial_potential_derivatives, z), axis=0)
        mag_coeffs = _quintic_spline_coefficients(z, mag_derivs.T)
        
        logging.log_info(f'Computing derivative interpolation took {(time.time()-st)*1000:.2f} ms ({len(z)} items)')

        return z, elec_coeffs, mag_coeffs
     
    def electrostatic_field_at_local_point(self, point: PointLike3D) -> Vector3D:
        """
        Compute the electric field, \\( \\vec{E} = -\\nabla \\phi \\)
        
        Parameters
        ----------
        point: (3,) array of float64
            Position in local coordinate system at which to compute the field.
             
        Returns
        -------
        Numpy array containing the field strengths (in units of V/mm) in the r and z directions.
        """
        point = np.array(point, dtype=np.float64)
        assert point.shape == (3,), "Please supply a three dimensional point"
        return backend.field_radial_derivs(point, self.z, self.electrostatic_coeffs)
    
    def magnetostatic_field_at_local_point(self, point: PointLike3D) -> Vector3D:
        """
        Compute the magnetic field \\( \\vec{H} \\)
        
        Parameters
        ----------
        point: (3,) array of float64
            Position in local coordinate system at which to compute the field.
             
        Returns
        -------
        (3,) np.ndarray of float64 containing the field strength (in units of A/m) in the x, y and z directions.
        """
        point = np.array(point, dtype=np.float64)
        assert point.shape == (3,), "Please supply a three dimensional point"
        return backend.field_radial_derivs(point, self.z, self.magnetostatic_coeffs)
     
    def electrostatic_potential_at_local_point(self, point: PointLike3D) -> float:
        """
        Compute the electrostatic potential (close to the axis).

        Parameters
        ----------
        point: (3,) array of float64
            Position at which to compute the potential.
        
        Returns
        -------
        Potential as a float value (in units of V).
        """
        point = np.array(point)
        assert point.shape == (3,), "Please supply a three dimensional point"
        return backend.potential_radial_derivs(point, self.z, self.electrostatic_coeffs)
    
    def magnetostatic_potential_at_local_point(self, point: PointLike3D) -> float:
        """
        Compute the magnetostatic scalar potential (satisfying \\(\\vec{H} = -\\nabla \\phi \\)) close to the axis
        
        Parameters
        ----------
        point: (3,) array of float64
            Position in local coordinate system at which to compute the potential.
        
        Returns
        -------
        Potential as a float value (in units of A).
        """
        point = np.array(point, dtype=np.float64)
        assert point.shape == (3,), "Please supply a three dimensional point"
        return backend.potential_radial_derivs(point, self.z, self.magnetostatic_coeffs)
    
    def get_tracer(self, bounds: BoundsLike3D) -> Tracer:
        return T.Tracer(self, bounds)
    
    def get_low_level_trace_function(self) -> tuple[Callable, Any]:
        args = backend.FieldDerivsArgs(self.z, self.electrostatic_coeffs, self.magnetostatic_coeffs)
        return backend.field_fun(("field_radial_derivs_traceable", backend.backend_lib)), args

An electrostatic field resulting from a radial series expansion around the optical axis. You should not initialize this class yourself, but it is used as a base class for the fields returned by the axial_derivative_interpolation methods. This base class overloads the +,*,- operators so it is very easy to take a superposition of different fields.

Produces a field which uses an axial interpolation to very quickly compute the field around the z-axis. Note that the approximation degrades as the point at which the field is computed is further from the z-axis. Also note that the fields produced by current and the magnetizable material are merged into one interpolation, so the method current_field_at_point will always return zero.

Parameters

field : FieldRadialBEM
Field for which to compute the axial interpolation
zmin : float
Location on the optical axis where to start sampling the radial expansion coefficients.
zmax : float
Location on the optical axis where to stop sampling the radial expansion coefficients. Any field evaluation outside [zmin, zmax] will return a zero field strength.
N : int, optional
Number of samples to take on the optical axis, if N=None the amount of samples is determined by taking into account the number of elements in the mesh.

Ancestors

Methods

def electrostatic_field_at_local_point(self, point)
Expand source code
def electrostatic_field_at_local_point(self, point: PointLike3D) -> Vector3D:
    """
    Compute the electric field, \\( \\vec{E} = -\\nabla \\phi \\)
    
    Parameters
    ----------
    point: (3,) array of float64
        Position in local coordinate system at which to compute the field.
         
    Returns
    -------
    Numpy array containing the field strengths (in units of V/mm) in the r and z directions.
    """
    point = np.array(point, dtype=np.float64)
    assert point.shape == (3,), "Please supply a three dimensional point"
    return backend.field_radial_derivs(point, self.z, self.electrostatic_coeffs)

Compute the electric field, \vec{E} = -\nabla \phi

Parameters

point : (3,) array of float64
Position in local coordinate system at which to compute the field.

Returns

Numpy array containing the field strengths (in units of V/mm) in the r and z directions.

def electrostatic_potential_at_local_point(self, point)
Expand source code
def electrostatic_potential_at_local_point(self, point: PointLike3D) -> float:
    """
    Compute the electrostatic potential (close to the axis).

    Parameters
    ----------
    point: (3,) array of float64
        Position at which to compute the potential.
    
    Returns
    -------
    Potential as a float value (in units of V).
    """
    point = np.array(point)
    assert point.shape == (3,), "Please supply a three dimensional point"
    return backend.potential_radial_derivs(point, self.z, self.electrostatic_coeffs)

Compute the electrostatic potential (close to the axis).

Parameters

point : (3,) array of float64
Position at which to compute the potential.

Returns

Potential as a float value (in units of V).

def get_tracer(self, bounds)
Expand source code
def get_tracer(self, bounds: BoundsLike3D) -> Tracer:
    return T.Tracer(self, bounds)
def magnetostatic_field_at_local_point(self, point)
Expand source code
def magnetostatic_field_at_local_point(self, point: PointLike3D) -> Vector3D:
    """
    Compute the magnetic field \\( \\vec{H} \\)
    
    Parameters
    ----------
    point: (3,) array of float64
        Position in local coordinate system at which to compute the field.
         
    Returns
    -------
    (3,) np.ndarray of float64 containing the field strength (in units of A/m) in the x, y and z directions.
    """
    point = np.array(point, dtype=np.float64)
    assert point.shape == (3,), "Please supply a three dimensional point"
    return backend.field_radial_derivs(point, self.z, self.magnetostatic_coeffs)

Compute the magnetic field \vec{H}

Parameters

point : (3,) array of float64
Position in local coordinate system at which to compute the field.

Returns

(3,) np.ndarray of float64 containing the field strength (in units of A/m) in the x, y and z directions.

def magnetostatic_potential_at_local_point(self, point)
Expand source code
def magnetostatic_potential_at_local_point(self, point: PointLike3D) -> float:
    """
    Compute the magnetostatic scalar potential (satisfying \\(\\vec{H} = -\\nabla \\phi \\)) close to the axis
    
    Parameters
    ----------
    point: (3,) array of float64
        Position in local coordinate system at which to compute the potential.
    
    Returns
    -------
    Potential as a float value (in units of A).
    """
    point = np.array(point, dtype=np.float64)
    assert point.shape == (3,), "Please supply a three dimensional point"
    return backend.potential_radial_derivs(point, self.z, self.magnetostatic_coeffs)

Compute the magnetostatic scalar potential (satisfying \vec{H} = -\nabla \phi ) close to the axis

Parameters

point : (3,) array of float64
Position in local coordinate system at which to compute the potential.

Returns

Potential as a float value (in units of A).

Inherited members

class FieldRadialBEM (electrostatic_point_charges=None,
magnetostatic_point_charges=None,
current_point_charges=None)
Expand source code
class FieldRadialBEM(FieldBEM):
    """A radially symmetric electrostatic field. The field is a result of the surface charges as computed by the
    `solve_direct` function. See the comments in `FieldBEM`."""
    
    def __init__(self, 
                 electrostatic_point_charges: EffectivePointCharges | None = None, 
                 magnetostatic_point_charges: EffectivePointCharges | None = None, 
                 current_point_charges: EffectivePointCharges | None = None) -> None:
        
        if electrostatic_point_charges is None:
            electrostatic_point_charges = EffectivePointCharges.empty_2d()
        if magnetostatic_point_charges is None:
            magnetostatic_point_charges = EffectivePointCharges.empty_2d()
        if current_point_charges is None:
            current_point_charges = EffectivePointCharges.empty_3d()
        
        assert all([isinstance(eff, EffectivePointCharges) for eff in [electrostatic_point_charges,
                                                                       magnetostatic_point_charges,
                                                                       current_point_charges]])
        self.symmetry = E.Symmetry.RADIAL
        super().__init__(electrostatic_point_charges, magnetostatic_point_charges, current_point_charges)
         
    def current_field_at_local_point(self, point: PointLike3D) -> Vector3D:
        point = np.array(point, dtype=np.float64)
        assert point.shape == (3,), "Please supply a three dimensional point"
            
        currents = self.current_point_charges.charges
        jacobians = self.current_point_charges.jacobians
        positions = self.current_point_charges.positions
        return backend.current_field_radial(point, currents, jacobians, positions)
     
    def electrostatic_field_at_local_point(self, point: PointLike3D) -> Vector3D:
        """
        Compute the electric field, \\( \\vec{E} = -\\nabla \\phi \\)
        
        Parameters
        ----------
        point: (3,) array of float64
            Position in local coordinate system at which to compute the field.
        
        Returns
        -------
        (3,) array of float64, containing the field strengths (units of V/m)
        """
        point = np.array(point, dtype=np.float64)
        assert point.shape == (3,), "Please supply a three dimensional point"
          
        charges = self.electrostatic_point_charges.charges
        jacobians = self.electrostatic_point_charges.jacobians
        positions = self.electrostatic_point_charges.positions
        return backend.field_radial(point, charges, jacobians, positions)
    
    def magnetostatic_field_at_local_point(self, point: PointLike3D) -> Vector3D:
        """
        Compute the magnetic field \\( \\vec{H} \\)
        
        Parameters
        ----------
        point: (3,) array of float64
            Position in local coordinate system at which to compute the field.
             
        Returns
        -------
        (3,) np.ndarray of float64 containing the field strength (in units of A/m) in the x, y and z directions.
        """
        point = np.array(point, dtype=np.float64)
        assert point.shape == (3,), "Please supply a three dimensional point"
        current_field = self.current_field_at_point(point)
        
        charges = self.magnetostatic_point_charges.charges
        jacobians = self.magnetostatic_point_charges.jacobians
        positions = self.magnetostatic_point_charges.positions
        
        mag_field = backend.field_radial(point, charges, jacobians, positions)

        return current_field + mag_field

    def electrostatic_potential_at_local_point(self, point: PointLike3D) -> float:
        """
        Compute the electrostatic potential.
        
        Parameters
        ----------
        point: (3,) array of float64
            Position in local coordinate system at which to compute the field.
        
        Returns
        -------
        Potential as a float value (in units of V).
        """
        point = np.array(point, dtype=np.float64)
        assert point.shape == (3,), "Please supply a three dimensional point"
        charges = self.electrostatic_point_charges.charges
        jacobians = self.electrostatic_point_charges.jacobians
        positions = self.electrostatic_point_charges.positions
        return backend.potential_radial(point, charges, jacobians, positions)
     
    def magnetostatic_potential_at_local_point(self, point: PointLike3D) -> float:
        """
        Compute the magnetostatic scalar potential (satisfying \\(\\vec{H} = -\\nabla \\phi \\))
        
        Parameters
        ----------
        point: (3,) array of float64
            Position in local coordinate system in local coordinate system at which to compute the field.
        
        Returns
        -------
        Potential as a float value (in units of A).
        """
        point = np.array(point, dtype=np.float64)
        assert point.shape == (3,), "Please supply a three dimensional point"
        charges = self.magnetostatic_point_charges.charges
        jacobians = self.magnetostatic_point_charges.jacobians
        positions = self.magnetostatic_point_charges.positions
        return backend.potential_radial(point, charges, jacobians, positions)
    
    def current_potential_axial(self, z: float) -> float:
        assert isinstance(z, float)
        currents = self.current_point_charges.charges
        jacobians = self.current_point_charges.jacobians
        positions = self.current_point_charges.positions
        return backend.current_potential_axial(z, currents, jacobians, positions)
     
    def get_electrostatic_axial_potential_derivatives(self, z: ArrayLikeFloat1D) -> ArrayFloat2D:
        """
        Compute the derivatives of the electrostatic potential a points on the optical axis (z-axis). 
         
        Parameters
        ----------
        z : (N,) np.ndarray of float64
            Positions on the optical axis at which to compute the derivatives.

        Returns
        ------- 
        Numpy array of shape (N, 9) containing the derivatives. At index i one finds the i-th derivative (so
        at position 0 the potential itself is returned). The highest derivative returned is a 
        constant currently set to 9."""
        z = np.array(z, dtype=np.float64)
        charges = self.electrostatic_point_charges.charges
        jacobians = self.electrostatic_point_charges.jacobians
        positions = self.electrostatic_point_charges.positions
        return backend.axial_derivatives_radial(z, charges, jacobians, positions)
    
    def get_magnetostatic_axial_potential_derivatives(self, z):
        """
        Compute the derivatives of the magnetostatic potential at points on the optical axis (z-axis). 
         
        Parameters
        ----------
        z : (N,) np.ndarray of float64
            Positions on the optical axis at which to compute the derivatives.

        Returns
        ------- 
        Numpy array of shape (N, 9) containing the derivatives. At index i one finds the i-th derivative (so
        at position 0 the potential itself is returned). The highest derivative returned is a 
        constant currently set to 9."""
        charges = self.magnetostatic_point_charges.charges
        jacobians = self.magnetostatic_point_charges.jacobians
        positions = self.magnetostatic_point_charges.positions
         
        derivs_magnetic = backend.axial_derivatives_radial(z, charges, jacobians, positions)
        derivs_current = self.get_current_axial_potential_derivatives(z)
        return derivs_magnetic + derivs_current
     
    def get_current_axial_potential_derivatives(self, z: ArrayLikeFloat1D) -> ArrayFloat2D:
        """
        Compute the derivatives of the current magnetostatic scalar potential at points on the optical axis.
         
        Parameters
        ----------
        z : (N,) np.ndarray of float64
            Positions on the optical axis at which to compute the derivatives.
         
        Returns
        ------- 
        Numpy array of shape (N, 9) containing the derivatives. At index i one finds the i-th derivative (so
        at position 0 the potential itself is returned). The highest derivative returned is a 
        constant currently set to 9."""
        z = np.array(z, dtype=np.float64)
        currents = self.current_point_charges.charges
        jacobians = self.current_point_charges.jacobians
        positions = self.current_point_charges.positions
        return backend.current_axial_derivatives_radial(z, currents, jacobians, positions)
      
    def area_of_element(self, i: int) -> float:
        jacobians = self.electrostatic_point_charges.jacobians
        positions = self.electrostatic_point_charges.positions
        return 2*np.pi*np.sum(jacobians[i] * positions[i, :, 0])
    
    def get_tracer(self, bounds: BoundsLike3D)-> Tracer:
        return T.Tracer(self, bounds)
    
    def get_low_level_trace_function(self) -> tuple[Callable, Any]:
        args = backend.FieldEvaluationArgsRadial(self.electrostatic_point_charges, self.magnetostatic_point_charges, self.current_point_charges, self.field_bounds)
        return backend.field_fun(("field_radial_traceable", backend.backend_lib)), args

A radially symmetric electrostatic field. The field is a result of the surface charges as computed by the solve_direct function. See the comments in FieldBEM.

Ancestors

Methods

def area_of_element(self, i)
Expand source code
def area_of_element(self, i: int) -> float:
    jacobians = self.electrostatic_point_charges.jacobians
    positions = self.electrostatic_point_charges.positions
    return 2*np.pi*np.sum(jacobians[i] * positions[i, :, 0])
def current_field_at_local_point(self, point)
Expand source code
def current_field_at_local_point(self, point: PointLike3D) -> Vector3D:
    point = np.array(point, dtype=np.float64)
    assert point.shape == (3,), "Please supply a three dimensional point"
        
    currents = self.current_point_charges.charges
    jacobians = self.current_point_charges.jacobians
    positions = self.current_point_charges.positions
    return backend.current_field_radial(point, currents, jacobians, positions)
def current_potential_axial(self, z)
Expand source code
def current_potential_axial(self, z: float) -> float:
    assert isinstance(z, float)
    currents = self.current_point_charges.charges
    jacobians = self.current_point_charges.jacobians
    positions = self.current_point_charges.positions
    return backend.current_potential_axial(z, currents, jacobians, positions)
def electrostatic_field_at_local_point(self, point)
Expand source code
def electrostatic_field_at_local_point(self, point: PointLike3D) -> Vector3D:
    """
    Compute the electric field, \\( \\vec{E} = -\\nabla \\phi \\)
    
    Parameters
    ----------
    point: (3,) array of float64
        Position in local coordinate system at which to compute the field.
    
    Returns
    -------
    (3,) array of float64, containing the field strengths (units of V/m)
    """
    point = np.array(point, dtype=np.float64)
    assert point.shape == (3,), "Please supply a three dimensional point"
      
    charges = self.electrostatic_point_charges.charges
    jacobians = self.electrostatic_point_charges.jacobians
    positions = self.electrostatic_point_charges.positions
    return backend.field_radial(point, charges, jacobians, positions)

Compute the electric field, \vec{E} = -\nabla \phi

Parameters

point : (3,) array of float64
Position in local coordinate system at which to compute the field.

Returns

(3,) array of float64, containing the field strengths (units of V/m)

def electrostatic_potential_at_local_point(self, point)
Expand source code
def electrostatic_potential_at_local_point(self, point: PointLike3D) -> float:
    """
    Compute the electrostatic potential.
    
    Parameters
    ----------
    point: (3,) array of float64
        Position in local coordinate system at which to compute the field.
    
    Returns
    -------
    Potential as a float value (in units of V).
    """
    point = np.array(point, dtype=np.float64)
    assert point.shape == (3,), "Please supply a three dimensional point"
    charges = self.electrostatic_point_charges.charges
    jacobians = self.electrostatic_point_charges.jacobians
    positions = self.electrostatic_point_charges.positions
    return backend.potential_radial(point, charges, jacobians, positions)

Compute the electrostatic potential.

Parameters

point : (3,) array of float64
Position in local coordinate system at which to compute the field.

Returns

Potential as a float value (in units of V).

def get_current_axial_potential_derivatives(self, z)
Expand source code
def get_current_axial_potential_derivatives(self, z: ArrayLikeFloat1D) -> ArrayFloat2D:
    """
    Compute the derivatives of the current magnetostatic scalar potential at points on the optical axis.
     
    Parameters
    ----------
    z : (N,) np.ndarray of float64
        Positions on the optical axis at which to compute the derivatives.
     
    Returns
    ------- 
    Numpy array of shape (N, 9) containing the derivatives. At index i one finds the i-th derivative (so
    at position 0 the potential itself is returned). The highest derivative returned is a 
    constant currently set to 9."""
    z = np.array(z, dtype=np.float64)
    currents = self.current_point_charges.charges
    jacobians = self.current_point_charges.jacobians
    positions = self.current_point_charges.positions
    return backend.current_axial_derivatives_radial(z, currents, jacobians, positions)

Compute the derivatives of the current magnetostatic scalar potential at points on the optical axis.

Parameters

z : (N,) np.ndarray of float64
Positions on the optical axis at which to compute the derivatives.

Returns

Numpy array of shape (N, 9) containing the derivatives. At index i one finds the i-th derivative (so at position 0 the potential itself is returned). The highest derivative returned is a constant currently set to 9.

def get_electrostatic_axial_potential_derivatives(self, z)
Expand source code
def get_electrostatic_axial_potential_derivatives(self, z: ArrayLikeFloat1D) -> ArrayFloat2D:
    """
    Compute the derivatives of the electrostatic potential a points on the optical axis (z-axis). 
     
    Parameters
    ----------
    z : (N,) np.ndarray of float64
        Positions on the optical axis at which to compute the derivatives.

    Returns
    ------- 
    Numpy array of shape (N, 9) containing the derivatives. At index i one finds the i-th derivative (so
    at position 0 the potential itself is returned). The highest derivative returned is a 
    constant currently set to 9."""
    z = np.array(z, dtype=np.float64)
    charges = self.electrostatic_point_charges.charges
    jacobians = self.electrostatic_point_charges.jacobians
    positions = self.electrostatic_point_charges.positions
    return backend.axial_derivatives_radial(z, charges, jacobians, positions)

Compute the derivatives of the electrostatic potential a points on the optical axis (z-axis).

Parameters

z : (N,) np.ndarray of float64
Positions on the optical axis at which to compute the derivatives.

Returns

Numpy array of shape (N, 9) containing the derivatives. At index i one finds the i-th derivative (so at position 0 the potential itself is returned). The highest derivative returned is a constant currently set to 9.

def get_magnetostatic_axial_potential_derivatives(self, z)
Expand source code
def get_magnetostatic_axial_potential_derivatives(self, z):
    """
    Compute the derivatives of the magnetostatic potential at points on the optical axis (z-axis). 
     
    Parameters
    ----------
    z : (N,) np.ndarray of float64
        Positions on the optical axis at which to compute the derivatives.

    Returns
    ------- 
    Numpy array of shape (N, 9) containing the derivatives. At index i one finds the i-th derivative (so
    at position 0 the potential itself is returned). The highest derivative returned is a 
    constant currently set to 9."""
    charges = self.magnetostatic_point_charges.charges
    jacobians = self.magnetostatic_point_charges.jacobians
    positions = self.magnetostatic_point_charges.positions
     
    derivs_magnetic = backend.axial_derivatives_radial(z, charges, jacobians, positions)
    derivs_current = self.get_current_axial_potential_derivatives(z)
    return derivs_magnetic + derivs_current

Compute the derivatives of the magnetostatic potential at points on the optical axis (z-axis).

Parameters

z : (N,) np.ndarray of float64
Positions on the optical axis at which to compute the derivatives.

Returns

Numpy array of shape (N, 9) containing the derivatives. At index i one finds the i-th derivative (so at position 0 the potential itself is returned). The highest derivative returned is a constant currently set to 9.

def get_tracer(self, bounds)
Expand source code
def get_tracer(self, bounds: BoundsLike3D)-> Tracer:
    return T.Tracer(self, bounds)
def magnetostatic_field_at_local_point(self, point)
Expand source code
def magnetostatic_field_at_local_point(self, point: PointLike3D) -> Vector3D:
    """
    Compute the magnetic field \\( \\vec{H} \\)
    
    Parameters
    ----------
    point: (3,) array of float64
        Position in local coordinate system at which to compute the field.
         
    Returns
    -------
    (3,) np.ndarray of float64 containing the field strength (in units of A/m) in the x, y and z directions.
    """
    point = np.array(point, dtype=np.float64)
    assert point.shape == (3,), "Please supply a three dimensional point"
    current_field = self.current_field_at_point(point)
    
    charges = self.magnetostatic_point_charges.charges
    jacobians = self.magnetostatic_point_charges.jacobians
    positions = self.magnetostatic_point_charges.positions
    
    mag_field = backend.field_radial(point, charges, jacobians, positions)

    return current_field + mag_field

Compute the magnetic field \vec{H}

Parameters

point : (3,) array of float64
Position in local coordinate system at which to compute the field.

Returns

(3,) np.ndarray of float64 containing the field strength (in units of A/m) in the x, y and z directions.

def magnetostatic_potential_at_local_point(self, point)
Expand source code
def magnetostatic_potential_at_local_point(self, point: PointLike3D) -> float:
    """
    Compute the magnetostatic scalar potential (satisfying \\(\\vec{H} = -\\nabla \\phi \\))
    
    Parameters
    ----------
    point: (3,) array of float64
        Position in local coordinate system in local coordinate system at which to compute the field.
    
    Returns
    -------
    Potential as a float value (in units of A).
    """
    point = np.array(point, dtype=np.float64)
    assert point.shape == (3,), "Please supply a three dimensional point"
    charges = self.magnetostatic_point_charges.charges
    jacobians = self.magnetostatic_point_charges.jacobians
    positions = self.magnetostatic_point_charges.positions
    return backend.potential_radial(point, charges, jacobians, positions)

Compute the magnetostatic scalar potential (satisfying \vec{H} = -\nabla \phi )

Parameters

point : (3,) array of float64
Position in local coordinate system in local coordinate system at which to compute the field.

Returns

Potential as a float value (in units of A).

Inherited members

class FieldSuperposition (fields, factors=None)
Expand source code
class FieldSuperposition(Field):
    """Representing a linear combination of fields (superposition). Will be automatically created if fields are added
    together (field1 + field2) and the underlying field classes do not implement a specialized add method."""
    
    def __init__(self, fields: Iterable[Field], factors: Iterable[float] | Iterable[np.floating] | None = None) -> None:
        super().__init__()
        
        assert all([isinstance(f, Field) for f in fields])
        self.fields: List[Field] = list(fields)
         
        self.factors: ArrayFloat1D = np.ones(len(self.fields)) if factors is None else np.array(factors, dtype=np.float64)

    @staticmethod
    def _concat_factors(f1: ArrayFloat1D, f2: ArrayFloat1D) -> ArrayFloat1D:
        return np.concatenate( (f1, f2) )

    def map_points(self, fun: Callable[[PointLike3D], Point3D]) -> FieldSuperposition:
        return FieldSuperposition([f.map_points(fun) for f in self.fields], self.factors)
    
    def electrostatic_field_at_local_point(self, point: PointLike3D) -> Vector3D:
        return np.sum([fa*f.electrostatic_field_at_point(point) for fa, f in zip(self.factors, self.fields)], axis=0)

    def magnetostatic_field_at_local_point(self, point: PointLike3D) -> Vector3D:
        return np.sum([fa*f.magnetostatic_field_at_point(point) for fa, f in zip(self.factors, self.fields)], axis=0)
    
    def current_field_at_local_point(self, point: PointLike3D) -> Vector3D:
        return np.sum([fa*f.current_field_at_point(point) for fa, f in zip(self.factors, self.fields)], axis=0)
    
    def electrostatic_potential_at_local_point(self, point: PointLike3D) -> float:
        return sum([fa.item()*f.electrostatic_potential_at_point(point) for fa, f in zip(self.factors, self.fields)])

    def magnetostatic_potential_at_local_point(self, point: PointLike3D) -> float:
        return sum([fa.item()*f.magnetostatic_potential_at_point(point) for fa, f in zip(self.factors, self.fields)])
    
    def is_electrostatic(self) -> bool:
        return any(f.is_electrostatic() for f in self.fields)
    
    def is_magnetostatic(self) -> bool:
        return any(f.is_magnetostatic() for f in self.fields)

    def get_tracer(self, bounds: BoundsLike3D) -> Tracer:
        return T.Tracer(self, bounds)

    def __add__(self, other: Field) -> FieldSuperposition:
        if isinstance(other, FieldSuperposition):
            return FieldSuperposition(self.fields + other.fields, FieldSuperposition._concat_factors(self.factors, other.factors) )
        else:
            return FieldSuperposition(self.fields + [other], FieldSuperposition._concat_factors(self.factors, np.array([1.])))

    def __radd__(self, other: Field) -> FieldSuperposition:
        return FieldSuperposition([other]+self.fields, [1.0]+list(self.factors))
     
    def __iadd__(self, other: Field) -> FieldSuperposition:
        self.fields = (self + other).fields
        return self

    def __mul__(self, other: float) -> FieldSuperposition:
        if _is_numeric(other):
            return FieldSuperposition(self.fields, other*self.factors)
        else:
            return NotImplemented
    
    def __rmul__(self, other: float) -> FieldSuperposition :
        return self.__mul__(other)
    
    def __getitem__(self, index: int | slice) -> Field:
        if isinstance(index, slice):
            fields: List[Field] = np.array(self.fields, dtype=object).__getitem__(index).tolist() # type: ignore
            return FieldSuperposition(fields, self.factors[index])
        elif isinstance(index, int):
            return self.factors[index] * self.fields[index]

        return NotImplemented
     
    def __len__(self) -> int:
        return len(self.fields)

    def __iter__(self) -> Iterator[Field]:
        for fa, f in zip(self.factors, self.fields):
            yield f*fa.item()
     
    def __str__(self) -> str:
        field_strs = ''.join(f'\n\t{f.__class__.__name__} (times factor {fa})' for fa, f in zip(self.factors, self.fields))
        return f"<FieldSuperposition with fields: {field_strs}>"

Representing a linear combination of fields (superposition). Will be automatically created if fields are added together (field1 + field2) and the underlying field classes do not implement a specialized add method.

Ancestors

Methods

def current_field_at_local_point(self, point)
Expand source code
def current_field_at_local_point(self, point: PointLike3D) -> Vector3D:
    return np.sum([fa*f.current_field_at_point(point) for fa, f in zip(self.factors, self.fields)], axis=0)
def electrostatic_field_at_local_point(self, point)
Expand source code
def electrostatic_field_at_local_point(self, point: PointLike3D) -> Vector3D:
    return np.sum([fa*f.electrostatic_field_at_point(point) for fa, f in zip(self.factors, self.fields)], axis=0)
def electrostatic_potential_at_local_point(self, point)
Expand source code
def electrostatic_potential_at_local_point(self, point: PointLike3D) -> float:
    return sum([fa.item()*f.electrostatic_potential_at_point(point) for fa, f in zip(self.factors, self.fields)])
def get_tracer(self, bounds)
Expand source code
def get_tracer(self, bounds: BoundsLike3D) -> Tracer:
    return T.Tracer(self, bounds)
def is_electrostatic(self)
Expand source code
def is_electrostatic(self) -> bool:
    return any(f.is_electrostatic() for f in self.fields)
def is_magnetostatic(self)
Expand source code
def is_magnetostatic(self) -> bool:
    return any(f.is_magnetostatic() for f in self.fields)
def magnetostatic_field_at_local_point(self, point)
Expand source code
def magnetostatic_field_at_local_point(self, point: PointLike3D) -> Vector3D:
    return np.sum([fa*f.magnetostatic_field_at_point(point) for fa, f in zip(self.factors, self.fields)], axis=0)
def magnetostatic_potential_at_local_point(self, point)
Expand source code
def magnetostatic_potential_at_local_point(self, point: PointLike3D) -> float:
    return sum([fa.item()*f.magnetostatic_potential_at_point(point) for fa, f in zip(self.factors, self.fields)])

Inherited members