Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added optika.apertures.CircularSectorAperture class. #76

Merged
merged 11 commits into from
Sep 22, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 34 additions & 0 deletions optika/_tests/test_apertures.py
Original file line number Diff line number Diff line change
Expand Up @@ -142,6 +142,40 @@ def test_radius(self, a: optika.apertures.CircularAperture):
assert np.all(a.radius >= 0)


@pytest.mark.parametrize(
argnames="a",
argvalues=[
optika.apertures.CircularSectorAperture(
radius=radius,
samples_wire=21,
active=active,
inverted=inverted,
transformation=transformation,
kwargs_plot=kwargs_plot,
)
for radius in radius_parameterization
for active in active_parameterization
for inverted in inverted_parameterization
for transformation in transform_parameterization
for kwargs_plot in test_mixins.kwargs_plot_parameterization
],
)
class TestCircularSectorAperture(
AbstractTestAbstractAperture,
):
def test_radius(self, a: optika.apertures.CircularSectorAperture):
assert isinstance(a.radius, (float, u.Quantity, na.AbstractScalar))
assert np.all(a.radius >= 0)

def test_angle_start(self, a: optika.apertures.CircularSectorAperture):
assert isinstance(a.radius, (float, u.Quantity, na.AbstractScalar))
assert np.all(a.radius >= 0)

def test_angle_stop(self, a: optika.apertures.CircularSectorAperture):
assert isinstance(a.radius, (float, u.Quantity, na.AbstractScalar))
assert np.all(a.radius >= 0)


class AbstractTestAbstractPolygonalAperture(
AbstractTestAbstractAperture,
):
Expand Down
174 changes: 174 additions & 0 deletions optika/apertures.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
__all__ = [
"AbstractAperture",
"CircularAperture",
"CircularSectorAperture",
"AbstractPolygonalAperture",
"PolygonalAperture",
"AbstractRegularPolygonalAperture",
Expand Down Expand Up @@ -304,6 +305,179 @@ def wire(self, num: None | int = None) -> na.Cartesian3dVectorArray:
return result


@dataclasses.dataclass(eq=False, repr=False)
class CircularSectorAperture(
AbstractAperture,
):
"""
A `circular sector <https://en.wikipedia.org/wiki/Circular_sector>`_
aperture.

Examples
--------

Plot a single circular aperture sector

.. jupyter-execute::

import matplotlib.pyplot as plt
import astropy.units as u
import astropy.visualization
import named_arrays as na
import optika

# Define a circular aperture sector
aperture = optika.apertures.CircularSectorAperture(
radius=50 * u.mm,
angle_start=-11 * u.deg,
angle_stop=40 * u.deg,
)

# Define points to sample the aperture with
points = na.Cartesian3dVectorLinearSpace(
start=aperture.bound_lower,
stop=aperture.bound_upper,
axis=na.Cartesian3dVectorArray("x", "y", "z"),
num=na.Cartesian3dVectorArray(11, 11, 1),
)

# Compute which points are inside the aperture
where = aperture(points)

# Plot the circular aperture sector
with astropy.visualization.quantity_support():
plt.figure()
plt.gca().set_aspect("equal")
aperture.plot(components=("x", "y"), color="black")
na.plt.scatter(
points.x,
points.y,
c=where.astype(float)
)
"""

radius: u.Quantity | na.AbstractScalar = 0 * u.mm
"""
The radius of the cirucular sector.
"""

angle_start: u.Quantity | na.AbstractScalar = 0 * u.deg
r"""
The starting angle of the circular sector.
Must be between :math:`-2 \pi` and :math:`+2 \pi` radians.
"""

angle_stop: u.Quantity | na.AbstractScalar = 180 * u.deg
"""
The ending angle of the circular sector.
Must be between :math:`-2 \pi` and :math:`+2 \pi` radians and
counterclockwise from `angle_start`.
"""

samples_wire: int = 101
active: bool | na.AbstractScalar = True
inverted: bool | na.AbstractScalar = False
transformation: None | na.transformations.AbstractTransformation = None
kwargs_plot: None | dict = None

@property
def shape(self) -> dict[str, int]:
return na.broadcast_shapes(
optika.shape(self.radius),
optika.shape(self.angle_start),
optika.shape(self.angle_stop),
optika.shape(self.active),
optika.shape(self.inverted),
optika.shape(self.transformation),
)

def __call__(
self,
position: na.AbstractCartesian3dVectorArray,
) -> na.AbstractScalar:
radius = self.radius
angle_start = self.angle_start
angle_stop = self.angle_stop
active = self.active
inverted = self.inverted
if self.transformation is not None:
position = self.transformation.inverse(position)

shape = na.shape_broadcasted(
radius,
angle_start,
angle_stop,
active,
inverted,
position,
)

radius = na.broadcast_to(radius, shape)
angle_start = na.broadcast_to(angle_start, shape)
angle_stop = na.broadcast_to(angle_stop, shape)
active = na.broadcast_to(active, shape)
inverted = na.broadcast_to(inverted, shape)
position = na.broadcast_to(position, shape)

mask_radius = position.xy.length <= radius

angle = np.arctan2(position.y, position.x)
angle_positive = angle % (+2 * np.pi * u.rad)
angle_negative = angle % (-2 * np.pi * u.rad)
mask_positive = (angle_start < angle_positive) & (angle_positive < angle_stop)
mask_negative = (angle_start < angle_negative) & (angle_negative < angle_stop)
mask_angle = mask_positive | mask_negative

mask = mask_radius & mask_angle

mask[inverted] = ~mask[inverted]
mask[~active] = True

return mask

@property
def bound_lower(self) -> na.Cartesian3dVectorArray:
return self.wire().min()

@property
def bound_upper(self) -> na.Cartesian3dVectorArray:
return self.wire().max()

@property
def vertices(self) -> None:
return None

def wire(self, num: None | int = None) -> na.Cartesian3dVectorArray:
if num is None:
num = self.samples_wire
az = na.linspace(
start=self.angle_start,
stop=self.angle_stop,
axis="wire",
num=num - 2,
)
unit_radius = na.unit(self.radius)
result = na.Cartesian3dVectorArray(
x=self.radius * np.cos(az),
y=self.radius * np.sin(az),
z=0 * unit_radius if unit_radius is not None else 0,
)

vertex = na.Cartesian3dVectorArray().add_axes("wire")
vertex = vertex * unit_radius if unit_radius is not None else vertex
result = np.concatenate(
[
vertex,
result,
vertex,
],
axis="wire",
)
if self.transformation is not None:
result = self.transformation(result)
return result


@dataclasses.dataclass(eq=False, repr=False)
class AbstractPolygonalAperture(
AbstractAperture,
Expand Down
Loading