diff --git a/docs/index.rst b/docs/index.rst index 37ae8ad..c8165e9 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -82,6 +82,7 @@ Jupyter notebook examples on how to use :mod:`optika`. :maxdepth: 1 tutorials/prime_focus + tutorials/grazing_spectrograph API Reference diff --git a/docs/tutorials/grazing_spectrograph.ipynb b/docs/tutorials/grazing_spectrograph.ipynb new file mode 100644 index 0000000..211ef21 --- /dev/null +++ b/docs/tutorials/grazing_spectrograph.ipynb @@ -0,0 +1,1036 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Grazing-Incidence Soft X-ray Spectrograph Tutorial\n", + "\n", + "Soft X-rays (here, $\\lambda = 15$–$20$ Å, $\\sim 0.6$–$0.8$ keV) reflect efficiently\n", + "only at *grazing* incidence, so an imaging soft X-ray telescope cannot use a\n", + "normal-incidence parabola. Instead it uses a **Wolter Type-I** optic: a\n", + "**paraboloid** followed by a confocal **hyperboloid**, both struck at a shallow\n", + "grazing angle. The two reflections satisfy the Abbe sine condition far better\n", + "than a single mirror, giving good off-axis imaging.\n", + "\n", + "This tutorial designs a slitless grazing-incidence spectrograph end-to-end in\n", + "[`optika`](https://optika.readthedocs.io):\n", + "\n", + "* a **Wolter-I telescope** (paraboloid + hyperboloid),\n", + "* a **5000 line/mm transmission grating** placed in the converging beam, which\n", + " disperses the first-order spectrum across the detector while the un-diffracted\n", + " zero order forms a direct image, and\n", + "* a **2k $\\times$ 4k CMOS detector** (10 µm pixels).\n", + "\n", + "The design is driven by two numbers: a **spatial plate scale of 1 arcsec/pixel**\n", + "and a **first-order spectral range of 15–20 Å**. Neither a grazing-incidence\n", + "telescope nor a transmission grating has been built in `optika` before, so this\n", + "is also a test of the package." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import astropy.units as u\n", + "import astropy.visualization\n", + "import named_arrays as na\n", + "import optika" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Focal length from the plate scale\n", + "\n", + "The plate scale (angle subtended by one pixel) of a focusing system is\n", + "$p = w_\\mathrm{pix} / f$. Solving for the focal length with a\n", + "$1''$/pixel scale and 10 µm pixels:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "width_pixel = 10 * u.um\n", + "plate_scale = 1 * u.arcsec\n", + "\n", + "focal_length = (width_pixel / plate_scale).to(\n", + " u.mm, equivalencies=u.dimensionless_angles()\n", + ")\n", + "focal_length.to(u.m)" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "## Aperture and grazing angle\n", + "\n", + "The collecting aperture is a design choice (here a 144 mm-diameter annulus).\n", + "Once the aperture radius $r_0$ and focal length $f$ are fixed, the grazing angle\n", + "follows from the Wolter-I geometry: a ray entering at radius $r_0$ must be bent\n", + "by a total angle $4\\alpha$ (twice $2\\alpha$, once at each mirror) to reach the\n", + "focus, so\n", + "\n", + "$$r_0 = f \\tan(4\\alpha) \\quad\\Longrightarrow\\quad\n", + " \\alpha = \\tfrac14 \\arctan(r_0 / f).$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "radius_aperture = 72 * u.mm\n", + "\n", + "grazing_angle = 0.25 * np.arctan(\n", + " (radius_aperture / focal_length).to_value(u.dimensionless_unscaled)\n", + ") * u.rad\n", + "grazing_angle.to(u.deg)" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "A half-degree grazing angle is comfortably below the critical angle for\n", + "total external reflection of gold at these wavelengths (computed in the\n", + "*Reflectivity* section below), and it keeps the beam slow ($\\approx f/14$),\n", + "which — as we will see — matters for the spectral resolution." + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## Paraboloid (primary)\n", + "\n", + "`optika` represents a parabola of revolution with\n", + "`optika.sags.ParabolicSag`, whose vertex sits at the local origin and whose\n", + "focus lies a distance equal to the `focal_length` argument from that vertex.\n", + "\n", + "For grazing incidence we illuminate the parabola far from its vertex, on its\n", + "steep, nearly-cylindrical flank. A ray parallel to the axis at radius $r_0$ is\n", + "deflected by $2\\alpha$ toward the parabola focus a distance\n", + "$F_\\mathrm{p} = r_0 / \\tan(2\\alpha)$ downstream. The corresponding parabola\n", + "parameter is $f_\\mathrm{p} = r_0^2 / (4 F_\\mathrm{p})$ — here only a fraction of\n", + "a millimetre, since the illuminated flank is so steep.\n", + "\n", + "We use a negative `focal_length` so the parabola opens toward $-z$ and the\n", + "$+z$-travelling rays reflect *inward*, and we translate it so the illuminated\n", + "annulus sits near $z = 0$ and the parabola focus lands at $z = F_\\mathrm{p}$.\n", + "\n", + "The aperture is an `optika.apertures.AnnularAperture` — a grazing mirror is\n", + "physically a thin **ring** with an empty centre, not a full disk. A mirror of\n", + "axial length $L$ spans $\\Delta r = L\\tan\\alpha$ in radius, so the ring half-width\n", + "is $L\\tan\\alpha/2$. Because the grazing flank is so steep, even a short mirror is\n", + "*axially* long ($\\sim 0.2$ m here), so we keep $L$ small and place the secondary\n", + "well downstream so the two mirrors do not overlap in $z$ — otherwise some rays\n", + "could not reach the secondary while travelling forward." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "focus_parabola = (radius_aperture / np.tan(2 * grazing_angle)).to(u.mm)\n", + "f_parabola = (radius_aperture**2 / (4 * focus_parabola)).to(u.mm)\n", + "print(f\"parabola focus F_p = {focus_parabola.to(u.mm):.1f}\")\n", + "print(f\"parabola param f_p = {f_parabola.to(u.um):.1f}\")\n", + "\n", + "# A grazing mirror is a thin *annulus* (a ring with a central hole), not a full\n", + "# disk. A mirror of axial length L spans Δr = L·tan(α) in radius, so the ring\n", + "# half-width is L·tan(α)/2. Because the grazing flank is so steep, even a short\n", + "# mirror is axially long; we keep it short so the primary and secondary do not\n", + "# overlap in z (see the hyperboloid placement below).\n", + "mirror_length = 200 * u.mm\n", + "radial_halfwidth = (mirror_length * np.tan(grazing_angle) / 2).to(u.mm)\n", + "clearance = radial_halfwidth + 10 * u.mm\n", + "\n", + "paraboloid = optika.surfaces.Surface(\n", + " name=\"paraboloid\",\n", + " sag=optika.sags.ParabolicSag(focal_length=-f_parabola),\n", + " material=optika.materials.Mirror(),\n", + " aperture=optika.apertures.AnnularAperture(\n", + " radius_inner=radius_aperture - clearance,\n", + " radius_outer=radius_aperture + clearance,\n", + " ),\n", + " transformation=na.transformations.Cartesian3dTranslation(z=focus_parabola),\n", + " is_pupil_stop=True,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "## Hyperboloid (secondary)\n", + "\n", + "The hyperboloid is **confocal** with the paraboloid: one of its two foci\n", + "coincides with the parabola focus $F_\\mathrm{p}$ (its *far* focus), and the\n", + "other is the real system focus $F_\\mathrm{sys}$, where the detector sits (its\n", + "*near* focus). A conic mirror images one focus perfectly onto the other, so the\n", + "converging cone aimed at $F_\\mathrm{p}$ is re-imaged to a perfect point at\n", + "$F_\\mathrm{sys}$ — on-axis the Wolter-I has *zero* geometric aberration.\n", + "\n", + "A confocal conic is fixed by its two foci plus one point it passes through. We\n", + "require the hyperboloid to intercept the converging beam at an axial position\n", + "`z_intercept`; the beam radius there sets the third constraint. From the foci\n", + "$z_N, z_F$ and the through-point we recover the standard conic parameters:\n", + "semi-major axis $a$ (from $2a = |d_F - d_N|$), eccentricity $e = c/a$ with\n", + "$c = (z_F - z_N)/2$, vertex radius of curvature $R = a(e^2 - 1)$, and conic\n", + "constant $k = -e^2 < -1$.\n", + "\n", + "> **Note.** The parabola's true focus is one $f_\\mathrm{p}$ *short* of its\n", + "> vertex translation, i.e. at $z = F_\\mathrm{p} - f_\\mathrm{p}$. Using that\n", + "> exact value as the hyperboloid's far focus is what makes the on-axis image\n", + "> collapse to a point." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "z_system_focus = focal_length # detector / near focus\n", + "z_parabola_focus = focus_parabola - f_parabola # exact parabola focus / far focus\n", + "\n", + "z_center = 0.5 * (z_system_focus + z_parabola_focus)\n", + "c_hyp = 0.5 * (z_parabola_focus - z_system_focus) # half focus separation = a*e\n", + "\n", + "# Intercept the converging beam here. The hyperboloid is also axially long, so\n", + "# we place it well downstream of the primary's z-extent: the two mirrors must\n", + "# not overlap in z, or part of the beam could not reach the secondary travelling\n", + "# forward.\n", + "z_intercept = 500 * u.mm\n", + "radius_intercept = radius_aperture * (z_parabola_focus - z_intercept) / z_parabola_focus\n", + "\n", + "d_far = np.sqrt(radius_intercept**2 + (z_intercept - z_parabola_focus)**2)\n", + "d_near = np.sqrt(radius_intercept**2 + (z_intercept - z_system_focus)**2)\n", + "\n", + "a_hyp = 0.5 * np.abs(d_far - d_near)\n", + "e_hyp = (c_hyp / a_hyp).to_value(u.dimensionless_unscaled)\n", + "radius_hyp = a_hyp * (e_hyp**2 - 1)\n", + "conic_hyp = -(e_hyp**2)\n", + "z_vertex_hyp = z_center - a_hyp\n", + "\n", + "print(f\"eccentricity e = {e_hyp:.6f}\")\n", + "print(f\"conic k = {conic_hyp:.6f}\")\n", + "print(f\"vertex radius R = {radius_hyp.to(u.mm):.4f}\")\n", + "print(f\"vertex z = {z_vertex_hyp.to(u.mm):.2f}\")\n", + "\n", + "hyperboloid = optika.surfaces.Surface(\n", + " name=\"hyperboloid\",\n", + " sag=optika.sags.ConicSag(radius=-radius_hyp, conic=conic_hyp),\n", + " material=optika.materials.Mirror(),\n", + " aperture=optika.apertures.AnnularAperture(\n", + " radius_inner=radius_intercept - clearance,\n", + " radius_outer=radius_intercept + clearance,\n", + " ),\n", + " transformation=na.transformations.Cartesian3dTranslation(z=z_vertex_hyp),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "## Transmission grating\n", + "\n", + "A 5000 line/mm grating has a groove period $d = 200$ nm. In first order the\n", + "diffraction angle is $\\sin\\beta = m\\lambda/d$. Placed a distance $\\Lambda$\n", + "upstream of the focus, the grating shifts the first-order image laterally by\n", + "$\\approx \\Lambda\\,\\lambda/d$, giving a linear dispersion\n", + "$\\mathrm{d}x/\\mathrm{d}\\lambda = \\Lambda/d$.\n", + "\n", + "We model it as a flat, **transmissive** surface (the default `Vacuum` material\n", + "passes light straight through) carrying a `ConstantRulingSpacing` whose `normal`\n", + "points along $x$ — so the spectrum disperses along the detector's long axis.\n", + "`optika` applies the grating equation automatically during propagation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "density = 5000 / u.mm\n", + "spacing_grating = (1 / density).to(u.nm)\n", + "z_grating = 800 * u.mm # just downstream of the secondary, in the converging beam\n", + "\n", + "# the grating must be large enough to pass the converging beam at z_grating\n", + "radius_grating = (\n", + " radius_aperture * (z_parabola_focus - z_grating) / z_parabola_focus + 10 * u.mm\n", + ")\n", + "\n", + "grating = optika.surfaces.Surface(\n", + " name=\"grating\",\n", + " aperture=optika.apertures.CircularAperture(radius_grating),\n", + " rulings=optika.rulings.Rulings(\n", + " spacing=optika.rulings.ConstantRulingSpacing(\n", + " constant=spacing_grating,\n", + " normal=na.Cartesian3dVectorArray(1, 0, 0),\n", + " ),\n", + " diffraction_order=1,\n", + " ),\n", + " transformation=na.transformations.Cartesian3dTranslation(z=z_grating),\n", + ")\n", + "\n", + "dispersion = ((z_system_focus - z_grating) / spacing_grating).to(u.mm / u.AA)\n", + "print(f\"grating clear radius = {radius_grating.to(u.mm):.1f}\")\n", + "print(f\"linear dispersion = {dispersion:.4f} ({(width_pixel / dispersion).to(u.mAA):.1f}/pixel)\")" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "## Detector\n", + "\n", + "A 2k $\\times$ 4k CMOS array with 10 µm pixels. Real detectors of this kind are\n", + "*buttable*, so we make the dispersion axis the long (4k $\\approx$ 41 mm) one to\n", + "give the spectrum room. We offset the array along the dispersion direction so\n", + "the zero-order image lands near one end and the first-order spectrum (which\n", + "disperses toward $-x$) fills the rest of the chip." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "num_pixel = na.Cartesian2dVectorArray(4096, 2048) # (dispersion, cross-dispersion)\n", + "detector_offset = -12 * u.mm\n", + "\n", + "sensor = optika.sensors.ImagingSensor(\n", + " name=\"sensor\",\n", + " width_pixel=width_pixel,\n", + " axis_pixel=na.Cartesian2dVectorArray(\"detector_x\", \"detector_y\"),\n", + " num_pixel=num_pixel,\n", + " transformation=na.transformations.Cartesian3dTranslation(\n", + " x=detector_offset, z=z_system_focus,\n", + " ),\n", + " is_field_stop=True,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "## Source\n", + "\n", + "The Sun is effectively at infinity, so the source is specified as a *direction*\n", + "rather than a physical aperture: a circular patch of sky $0.5^\\circ$ across.\n", + "Following the convention used elsewhere in `optika`, giving the aperture radius\n", + "as a **dimensionless cosine** marks the surface as being at infinity." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16", + "metadata": {}, + "outputs": [], + "source": [ + "radius_sun = 0.25 * u.deg # 0.5 deg full angular diameter\n", + "\n", + "# Place the source well upstream of the telescope. The paraboloid's outer (most\n", + "# upstream) illuminated radius sits near z = -0.3 m on the steep flank, so a\n", + "# source at, say, z = -0.2 m would leave those rays travelling *backwards* to\n", + "# reach the mirror. z = -0.7 m keeps every ray moving forward (+z).\n", + "source = optika.surfaces.Surface(\n", + " name=\"solar disk\",\n", + " aperture=optika.apertures.CircularAperture(radius=np.cos(radius_sun)),\n", + " transformation=na.transformations.Cartesian3dTranslation(z=-700 * u.mm),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "## Input rays\n", + "\n", + "The entrance pupil of a Wolter-I telescope is a thin **annulus** — only the\n", + "grazing ring is illuminated, the centre is empty. We therefore build the pupil\n", + "grid directly in physical coordinates as a ring at radius $r_0$, with a small\n", + "radial half-width set by the mirror's axial length\n", + "($\\Delta r \\approx L_\\mathrm{mirror}\\tan\\alpha / 2$)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "# the pupil grid is the illuminated ring itself: radius r0 +/- radial_halfwidth\n", + "# (the mirror extent defined with the paraboloid above).\n", + "azimuth = na.linspace(0, 2 * np.pi, axis=\"pupil_az\", num=24, endpoint=False) * u.rad\n", + "radius_pupil = na.linspace(\n", + " radius_aperture - radial_halfwidth,\n", + " radius_aperture + radial_halfwidth,\n", + " axis=\"pupil_r\",\n", + " num=5,\n", + ")\n", + "pupil = na.Cartesian2dVectorArray(\n", + " x=radius_pupil * np.cos(azimuth),\n", + " y=radius_pupil * np.sin(azimuth),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "The wavelength axis samples the 15–20 Å band, and a small angular field\n", + "grid (a few arcminutes) will be used for the spot diagrams." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [ + "wavelength = na.linspace(15, 20, axis=\"wavelength\", num=6) * u.AA\n", + "\n", + "# the detector half-height corresponds to this angular half-field ...\n", + "half_field_detector = np.arctan(\n", + " (num_pixel.y * width_pixel / (2 * focal_length)).to_value(u.dimensionless_unscaled)\n", + ") * u.rad\n", + "print(f\"detector half-field = {half_field_detector.to(u.arcmin):.1f}\")\n", + "\n", + "# ... but the well-corrected field is only a few arcmin, so sample that for the\n", + "# spot diagrams.\n", + "field_halfwidth = 5 * u.arcmin\n", + "field = na.Cartesian2dVectorLinearSpace(\n", + " start=-field_halfwidth,\n", + " stop=field_halfwidth,\n", + " axis=na.Cartesian2dVectorArray(\"fx\", \"fy\"),\n", + " num=3,\n", + " centers=True,\n", + ")\n", + "\n", + "grid_input = optika.vectors.ObjectVectorArray(\n", + " wavelength=wavelength, field=field, pupil=pupil,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "21", + "metadata": {}, + "source": [ + "## Building the system\n", + "\n", + "`optika.systems.SequentialSystem` assembles the surfaces, the source, the\n", + "sensor, and the input grid.\n", + "\n", + "Because the entrance pupil is an annulus specified in **physical** coordinates\n", + "(not the normalized $[-1, 1]^2$ square that the automatic backward-raytrace\n", + "assumes), we compute the default ray function directly with\n", + "`normalized_*=False` and cache it — the same workaround used in the Fresnel\n", + "zone-plate tutorial for a non-standard pupil." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22", + "metadata": {}, + "outputs": [], + "source": [ + "system = optika.systems.SequentialSystem(\n", + " object=source,\n", + " surfaces=[paraboloid, hyperboloid, grating],\n", + " sensor=sensor,\n", + " grid_input=grid_input,\n", + ")\n", + "\n", + "rays_default = system.rayfunction(\n", + " wavelength=wavelength,\n", + " field=field,\n", + " pupil=pupil,\n", + " normalized_field=False,\n", + " normalized_pupil=False,\n", + ")\n", + "system.__dict__[\"rayfunction_default\"] = rays_default\n", + "print(\"object at infinity:\", system.object_is_at_infinity)" + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "## Layout\n", + "\n", + "We plot the on-axis ray bundle in the cross-dispersion ($z$–$y$) plane. The beam\n", + "enters from the left, grazes the paraboloid and hyperboloid (the shallow double\n", + "fold near the top and bottom of the aperture), passes through the grating, and\n", + "converges to the focus.\n", + "\n", + "The mirrors are drawn (black) as their **meridional sag profiles** — the actual\n", + "surface cross-section over the illuminated radial band, found from the on-axis\n", + "trace. (`system.plot` instead draws each surface's aperture *rim*; on these\n", + "near-cylindrical grazing sags the inner and outer rims fall hundreds of mm apart\n", + "in $z$, so the rim outline is not a faithful picture of the optic.)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24", + "metadata": {}, + "outputs": [], + "source": [ + "raytrace_layout = system.raytrace(\n", + " wavelength=17.5 * u.AA,\n", + " field=na.Cartesian2dVectorArray(0, 0) * u.deg,\n", + " pupil=pupil,\n", + " normalized_field=False,\n", + " normalized_pupil=False,\n", + ")\n", + "\n", + "\n", + "def hit_radius_range(index_surface):\n", + " \"\"\"Radial extent illuminated on a surface, from the on-axis raytrace.\"\"\"\n", + " position = raytrace_layout.outputs.position[{system.axis_surface: index_surface}]\n", + " radius = position.xy.length\n", + " return radius.min(), radius.max()\n", + "\n", + "\n", + "def plot_mirror(ax, surface, radius_range, component, num=51, **kwargs):\n", + " \"\"\"Draw a mirror's meridional sag profile in the (z, ``component``) plane.\"\"\"\n", + " radius = na.linspace(radius_range[0], radius_range[1], axis=\"profile\", num=num)\n", + " for sign in [+1, -1]:\n", + " coordinate = sign * radius\n", + " zero = 0 * u.mm\n", + " if component == \"y\":\n", + " point = na.Cartesian3dVectorArray(x=zero, y=coordinate, z=zero)\n", + " else:\n", + " point = na.Cartesian3dVectorArray(x=coordinate, y=zero, z=zero)\n", + " point = na.Cartesian3dVectorArray(x=point.x, y=point.y, z=surface.sag(point))\n", + " if surface.transformation is not None:\n", + " point = surface.transformation(point)\n", + " na.plt.plot(point, ax=ax, axis=\"profile\", components=(\"z\", component), **kwargs)\n", + "\n", + "\n", + "# surface order in the trace: object=0, paraboloid=1, hyperboloid=2, grating=3, sensor=4\n", + "with astropy.visualization.quantity_support():\n", + " fig, ax = plt.subplots(figsize=(9, 4), constrained_layout=True)\n", + " na.plt.plot(\n", + " raytrace_layout.outputs.position,\n", + " ax=ax,\n", + " axis=system.axis_surface,\n", + " components=(\"z\", \"y\"),\n", + " color=\"tab:blue\",\n", + " alpha=0.4,\n", + " )\n", + " plot_mirror(ax, paraboloid, hit_radius_range(1), \"y\", color=\"black\", linewidth=2)\n", + " plot_mirror(ax, hyperboloid, hit_radius_range(2), \"y\", color=\"black\", linewidth=2)\n", + " ax.set_xlabel(f\"$z$ ({u.mm:latex_inline})\")\n", + " ax.set_ylabel(f\"$y$ ({u.mm:latex_inline})\")\n", + " ax.set_title(\"Cross-dispersion plane: mirror profiles (black) and rays\")" + ] + }, + { + "cell_type": "markdown", + "id": "25", + "metadata": {}, + "source": [ + "### Dispersion layout\n", + "\n", + "The layout above is the cross-dispersion ($z$–$y$) plane, where all diffraction\n", + "orders overlap. To *see* the dispersion we switch to the $z$–$x$ plane and trace\n", + "**three orders** through the grating: the un-diffracted **zero order** ($m = 0$,\n", + "grey), which forms the direct image at $x \\approx 0$; the **first order**\n", + "($m = 1$, solid); and the **second order** ($m = 2$, dashed), each traced at 15,\n", + "17.5, and 20 Å. The grating fans the orders out to different positions on the\n", + "detector — second order at twice the deflection of first. The annular pupil\n", + "appears as two pencils, one from each side of the ring, and the mirror profiles\n", + "are again drawn in black. The right panel zooms in on the focal plane, where the\n", + "three orders land in distinct strips along the black detector bar." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "26", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib as mpl\n", + "\n", + "# meridional (dispersion-plane) slice of the annular pupil: two pencils at +/- r0\n", + "azimuth_layout = na.linspace(0, np.pi, axis=\"pupil_az\", num=2) * u.rad\n", + "pupil_layout = na.Cartesian2dVectorArray(\n", + " x=radius_pupil * np.cos(azimuth_layout),\n", + " y=radius_pupil * np.sin(azimuth_layout),\n", + ")\n", + "\n", + "\n", + "def raytrace_order(wavelength, order):\n", + " grating_order = optika.surfaces.Surface(\n", + " name=\"grating\",\n", + " aperture=optika.apertures.CircularAperture(radius_grating),\n", + " rulings=optika.rulings.Rulings(\n", + " spacing=optika.rulings.ConstantRulingSpacing(\n", + " constant=spacing_grating, normal=na.Cartesian3dVectorArray(1, 0, 0),\n", + " ),\n", + " diffraction_order=order,\n", + " ),\n", + " transformation=na.transformations.Cartesian3dTranslation(z=z_grating),\n", + " )\n", + " sys_order = optika.systems.SequentialSystem(\n", + " object=source, surfaces=[paraboloid, hyperboloid, grating_order],\n", + " sensor=sensor, grid_input=grid_input,\n", + " )\n", + " return sys_order.raytrace(\n", + " wavelength=wavelength,\n", + " field=na.Cartesian2dVectorArray(0, 0) * u.deg,\n", + " pupil=pupil_layout, normalized_field=False, normalized_pupil=False,\n", + " )\n", + "\n", + "\n", + "cmap = mpl.colormaps[\"viridis\"]\n", + "norm = mpl.colors.Normalize(15, 20)\n", + "\n", + "\n", + "def draw(ax):\n", + " # zero order is achromatic: straight through to the direct image near x = 0\n", + " na.plt.plot(\n", + " raytrace_order(17.5 * u.AA, order=0).outputs.position, ax=ax,\n", + " axis=system.axis_surface, components=(\"z\", \"x\"),\n", + " color=\"0.6\", linewidth=0.8,\n", + " )\n", + " # first order (solid) and second order (dashed), fanned out by wavelength.\n", + " # (pass colours as strings: na.plt.plot would mistake an RGBA tuple for a\n", + " # length-4 array.)\n", + " for order, linestyle in [(1, \"-\"), (2, \"--\")]:\n", + " for lam_value in [15.0, 17.5, 20.0]:\n", + " na.plt.plot(\n", + " raytrace_order(lam_value * u.AA, order=order).outputs.position, ax=ax,\n", + " axis=system.axis_surface, components=(\"z\", \"x\"),\n", + " color=mpl.colors.to_hex(cmap(norm(lam_value))),\n", + " linewidth=0.8, linestyle=linestyle,\n", + " )\n", + " # detector extent along the dispersion axis\n", + " x_half = (num_pixel.x * width_pixel / 2).to(u.mm)\n", + " ax.plot(\n", + " u.Quantity([z_system_focus, z_system_focus]).to(u.mm),\n", + " u.Quantity([detector_offset - x_half, detector_offset + x_half]).to(u.mm),\n", + " color=\"black\", linewidth=4, solid_capstyle=\"butt\", zorder=5,\n", + " )\n", + "\n", + "\n", + "with astropy.visualization.quantity_support():\n", + " fig, (ax, axz) = plt.subplots(\n", + " 1, 2, figsize=(12, 4), gridspec_kw=dict(width_ratios=[2, 1]),\n", + " constrained_layout=True,\n", + " )\n", + " draw(ax)\n", + " draw(axz)\n", + "\n", + " # add the mirror profiles to the full panel (they sit outside the zoom window)\n", + " plot_mirror(ax, paraboloid, hit_radius_range(1), \"x\", color=\"black\", linewidth=1.5)\n", + " plot_mirror(ax, hyperboloid, hit_radius_range(2), \"x\", color=\"black\", linewidth=1.5)\n", + "\n", + " # zoom the right panel on the focal region to separate the orders\n", + " axz.set_xlim(1850 * u.mm, 2100 * u.mm)\n", + " axz.set_ylim(-28 * u.mm, 6 * u.mm)\n", + "\n", + " sm = mpl.cm.ScalarMappable(norm=norm, cmap=cmap)\n", + " sm.set_array([])\n", + " fig.colorbar(sm, ax=axz, label=f\"wavelength ({u.AA:latex_inline})\")\n", + "\n", + " for a in (ax, axz):\n", + " a.set_xlabel(f\"$z$ ({u.mm:latex_inline})\")\n", + " ax.set_ylabel(f\"$x$ ({u.mm:latex_inline}) — dispersion\")\n", + " ax.set_title(\"Full system: $m=0$ (grey), $m=1$ (solid), $m=2$ (dashed)\")\n", + " axz.set_title(\"Detector region (zoom)\")" + ] + }, + { + "cell_type": "markdown", + "id": "27", + "metadata": {}, + "source": [ + "## Aperture footprints\n", + "\n", + "A useful check is *where* the rays actually land on each optic relative to its\n", + "clear aperture. Tracing the full field and pupil (at 17.5 Å), the blue points\n", + "below are the ray intercepts and the black outline is each aperture. The two\n", + "grazing mirrors are illuminated in a thin **annular footprint** that sits\n", + "comfortably inside their annular stops (the central hole is the obscuration),\n", + "while the still-converging beam under-fills the circular transmission grating." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28", + "metadata": {}, + "outputs": [], + "source": [ + "# trace the full field + pupil (one wavelength) to map the footprint on each optic\n", + "raytrace_footprint = system.raytrace(\n", + " wavelength=17.5 * u.AA,\n", + " field=field,\n", + " pupil=pupil,\n", + " normalized_field=False,\n", + " normalized_pupil=False,\n", + ")\n", + "\n", + "# (index in the surface axis, surface, label)\n", + "footprint_surfaces = [\n", + " (1, paraboloid, \"paraboloid\"),\n", + " (2, hyperboloid, \"hyperboloid\"),\n", + " (3, grating, \"grating\"),\n", + "]\n", + "\n", + "with astropy.visualization.quantity_support():\n", + " fig, axs = plt.subplots(1, 3, figsize=(12, 4.5), constrained_layout=True)\n", + " for (index, surface, name), ax in zip(footprint_surfaces, axs):\n", + " position = raytrace_footprint.outputs.position[{system.axis_surface: index}]\n", + " unvignetted = raytrace_footprint.outputs.unvignetted[{system.axis_surface: index}]\n", + " na.plt.scatter(\n", + " position.x.to(u.mm), position.y.to(u.mm),\n", + " ax=ax, s=2, where=unvignetted, color=\"tab:blue\",\n", + " )\n", + " na.plt.plot(\n", + " surface.aperture.wire(), ax=ax, axis=\"wire\",\n", + " components=(\"x\", \"y\"), color=\"black\",\n", + " )\n", + " ax.set_aspect(\"equal\")\n", + " ax.set_title(name)\n", + " ax.set_xlabel(f\"$x$ ({u.mm:latex_inline})\")\n", + " axs[0].set_ylabel(f\"$y$ ({u.mm:latex_inline})\")" + ] + }, + { + "cell_type": "markdown", + "id": "29", + "metadata": {}, + "source": [ + "## Dispersion\n", + "\n", + "Tracing on-axis at 15, 17.5, and 20 Å, the first-order images fall in a line\n", + "along the dispersion ($x$) axis, evenly spaced at the linear dispersion derived\n", + "above ($\\approx 0.63$ mm/Å — coarser than a compact layout would give, because\n", + "the longer, non-overlapping telescope shortens the grating-to-focus lever arm).\n", + "\n", + "The positions below are in **sensor-local** coordinates. Because we offset the\n", + "detector by $-12$ mm, the un-diffracted **zero order** lands at $x = +12$ mm —\n", + "near one end of the 41 mm-long axis — while the 15–20 Å first-order band sits\n", + "near $x \\approx 0$. Longer wavelengths disperse further toward $-x$, so the\n", + "buttable 4k axis records the zero-order direct image and a broad first-order\n", + "spectrum at the same time.\n", + "\n", + "The same grating also diffracts into higher orders. The **second order**\n", + "($m = 2$) of the 15–20 Å band lands further out (near $x \\approx -7$ to $-13$ mm)\n", + "— still on the chip, and with *twice* the dispersion, giving roughly twice the\n", + "spectral resolution. (Third order would start to run off the detector past\n", + "$\\sim 16$ Å.) In a slitless spectrograph the orders overlap *spectrally* —\n", + "second order of $\\lambda$ falls where first order of $2\\lambda$ would — so for an\n", + "isolated band like this one the higher order is a free, higher-resolution copy of\n", + "the same spectrum in a separate strip of the detector." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30", + "metadata": {}, + "outputs": [], + "source": [ + "def centroid_x(wavelength, order=1):\n", + " grating_order = optika.surfaces.Surface(\n", + " name=\"grating\",\n", + " aperture=optika.apertures.CircularAperture(radius_grating),\n", + " rulings=optika.rulings.Rulings(\n", + " spacing=optika.rulings.ConstantRulingSpacing(\n", + " constant=spacing_grating, normal=na.Cartesian3dVectorArray(1, 0, 0),\n", + " ),\n", + " diffraction_order=order,\n", + " ),\n", + " transformation=na.transformations.Cartesian3dTranslation(z=z_grating),\n", + " )\n", + " sys_order = optika.systems.SequentialSystem(\n", + " object=source, surfaces=[paraboloid, hyperboloid, grating_order],\n", + " sensor=sensor, grid_input=grid_input,\n", + " )\n", + " out = sys_order.rayfunction(\n", + " wavelength=wavelength,\n", + " field=na.Cartesian2dVectorArray(0, 0) * u.deg,\n", + " pupil=pupil, normalized_field=False, normalized_pupil=False,\n", + " ).outputs\n", + " x = out.position.x.to(u.mm).ndarray.value\n", + " where = out.unvignetted.ndarray\n", + " if not where.any():\n", + " return np.nan # diffracted off the detector\n", + " return np.mean(x[where]) # mm, unvignetted rays only\n", + "\n", + "print(f\" zero order: x = {centroid_x(17.5 * u.AA, order=0):8.3f} mm\")\n", + "for order in [1, 2]:\n", + " print(f\"order m = {order}:\")\n", + " for lam in [15, 17.5, 20] * u.AA:\n", + " print(f\" {lam:>5}: x = {centroid_x(lam, order=order):8.3f} mm\")" + ] + }, + { + "cell_type": "markdown", + "id": "31", + "metadata": {}, + "source": [ + "## Spot diagrams\n", + "\n", + "`spot_diagram` shows the geometric PSF at each of the nine field points; the\n", + "grating dispersion is removed (each spot is re-centred on its own centroid) and\n", + "the points are coloured by wavelength. The **on-axis** panel (centre) collapses\n", + "to a near point — the Wolter-I is aberration-free on-axis — and the six\n", + "wavelengths sit almost on top of one another, since the grating adds only a\n", + "small blur in the converging beam.\n", + "\n", + "Moving off-axis, the spots quickly grow into the comatic spirals characteristic\n", + "of a two-mirror grazing telescope. Even at the few-arcminute field sampled here\n", + "the blur reaches hundreds of microns (tens of pixels), so the off-axis\n", + "aberrations — not the grating — set the usable field of view, which for this\n", + "simple design is only a couple of arcminutes wide." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = system.spot_diagram()\n", + "fig.set_size_inches(8, 8)" + ] + }, + { + "cell_type": "markdown", + "id": "33", + "metadata": {}, + "source": [ + "## Spatial plate scale\n", + "\n", + "A clean check of the plate scale uses the **cross-dispersion** ($y$) direction,\n", + "which the grating leaves untouched. Displacing the source by an angle\n", + "$\\theta_y$ shifts the image by $f_\\mathrm{eff}\\,\\theta_y$, so both the effective\n", + "focal length and the plate scale $w_\\mathrm{pix}/f_\\mathrm{eff}$ follow from the\n", + "measured shift.\n", + "\n", + "For a *two-mirror* telescope the effective focal length is measured from the rear\n", + "principal plane, which sits down among the mirrors — not from the coordinate\n", + "origin we used to place the focus. Now that the primary and secondary are spread\n", + "out over $\\sim 0.7$ m so they don't overlap, that principal plane sits well\n", + "inside the telescope, and the measured $f_\\mathrm{eff} \\approx 1.78$ m comes out\n", + "$\\sim 14\\%$ short of the 2.06 m back-focal distance. The plate scale is therefore\n", + "$\\approx 1.16''$/pixel rather than the $1''$/pixel target — a direct consequence\n", + "of the long grazing telescope. Moving the detector out $\\sim 0.3$ m (or\n", + "rescaling the whole design) would recover it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34", + "metadata": {}, + "outputs": [], + "source": [ + "theta_y = 30 * u.arcsec\n", + "\n", + "out_on = system.rayfunction(\n", + " wavelength=17.5 * u.AA, field=na.Cartesian2dVectorArray(0, 0) * u.deg,\n", + " pupil=pupil, normalized_field=False, normalized_pupil=False,\n", + ").outputs\n", + "out_off = system.rayfunction(\n", + " wavelength=17.5 * u.AA,\n", + " field=na.Cartesian2dVectorArray(0 * u.deg, theta_y.to(u.deg)),\n", + " pupil=pupil, normalized_field=False, normalized_pupil=False,\n", + ").outputs\n", + "\n", + "dy = np.nanmean(out_off.position.y.ndarray.value) - np.nanmean(out_on.position.y.ndarray.value)\n", + "dy = np.abs(dy) * out_on.position.y.unit\n", + "\n", + "f_eff = (dy / np.tan(theta_y)).to(u.m)\n", + "measured_scale = (theta_y / (dy / width_pixel)).to(u.arcsec)\n", + "\n", + "print(f\"image shift = {dy.to(u.um):.2f} ({(dy / width_pixel).to_value(u.dimensionless_unscaled):.2f} pixels)\")\n", + "print(f\"effective focal len = {f_eff:.3f} (back-focal distance was {focal_length.to(u.m):.3f})\")\n", + "print(f\"measured plate scale = {measured_scale:.3f} / pixel (target 1.000 arcsec / pixel)\")" + ] + }, + { + "cell_type": "markdown", + "id": "35", + "metadata": {}, + "source": [ + "## Mirror reflectivity at grazing incidence\n", + "\n", + "The geometric model above assumed perfect mirrors. In reality the grazing-angle\n", + "reflectivity sets the throughput. Using `optika`'s tabulated optical constants\n", + "$n = 1 - \\delta + i\\beta$ (CXRO data), the critical angle for total external\n", + "reflection is $\\theta_c = \\sqrt{2\\delta}$, and the single-surface reflectivity\n", + "follows from the Fresnel equation\n", + "\n", + "$$R(\\theta) = \\left|\\frac{\\sin\\theta - \\sqrt{n^2 - \\cos^2\\theta}}\n", + " {\\sin\\theta + \\sqrt{n^2 - \\cos^2\\theta}}\\right|^2 ,$$\n", + "\n", + "with $\\theta$ measured from the *surface*." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36", + "metadata": {}, + "outputs": [], + "source": [ + "au = optika.chemicals.Chemical(\"Au\")\n", + "ir = optika.chemicals.Chemical(\"Ir\")\n", + "\n", + "def reflectivity(chem, graze, wavelength):\n", + " n = chem.n(wavelength)\n", + " s, c = np.sin(graze), np.cos(graze)\n", + " root = np.sqrt(n**2 - c**2)\n", + " return np.abs((s - root) / (s + root)) ** 2\n", + "\n", + "# critical angle of gold across the band\n", + "delta_au = 1 - np.real(au.n(wavelength).ndarray)\n", + "theta_c = (np.sqrt(2 * delta_au) * u.rad).to(u.deg)\n", + "print(f\"Au critical angle over 15-20 A: {theta_c.min():.2f} - {theta_c.max():.2f}\")\n", + "print(f\"chosen grazing angle: {grazing_angle.to(u.deg):.2f}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37", + "metadata": {}, + "outputs": [], + "source": [ + "graze = na.linspace(0, 4, axis=\"graze\", num=400) * u.deg\n", + "\n", + "with astropy.visualization.quantity_support():\n", + " fig, (axL, axR) = plt.subplots(1, 2, figsize=(11, 4), constrained_layout=True)\n", + "\n", + " # reflectivity vs grazing angle at mid-band\n", + " for chem, name, color in [(au, \"Au\", \"tab:orange\"), (ir, \"Ir\", \"tab:purple\")]:\n", + " R = reflectivity(chem, graze, 17.5 * u.AA)\n", + " na.plt.plot(graze, R, ax=axL, axis=\"graze\", label=name, color=color)\n", + " axL.axvline(grazing_angle.to_value(u.deg), color=\"k\", ls=\"--\", lw=1,\n", + " label=f\"design ({grazing_angle.to(u.deg):.2f})\")\n", + " axL.set_xlabel(\"grazing angle (deg)\")\n", + " axL.set_ylabel(\"reflectivity (single bounce)\")\n", + " axL.set_title(r\"Reflectivity vs angle ($\\lambda = 17.5$ Å)\")\n", + " axL.legend()\n", + "\n", + " # two-bounce throughput vs wavelength at the design angle\n", + " wl = na.linspace(15, 20, axis=\"w\", num=101) * u.AA\n", + " for chem, name, color in [(au, \"Au\", \"tab:orange\"), (ir, \"Ir\", \"tab:purple\")]:\n", + " R = reflectivity(chem, grazing_angle, wl)\n", + " na.plt.plot(wl, R**2, ax=axR, axis=\"w\", label=name, color=color)\n", + " axR.set_xlabel(f\"wavelength ({u.AA:latex_inline})\")\n", + " axR.set_ylabel(\"throughput (two bounces)\")\n", + " axR.set_title(f\"Telescope throughput at {grazing_angle.to(u.deg):.2f}\")\n", + " axR.legend()" + ] + }, + { + "cell_type": "markdown", + "id": "38", + "metadata": {}, + "source": [ + "## Summary\n", + "\n", + "`optika` handled the full grazing-incidence spectrograph without modification:\n", + "\n", + "* **Wolter-I geometry.** A paraboloid (`ParabolicSag`) and a confocal hyperboloid\n", + " (`ConicSag`, $k < -1$) illuminated at large radius give exact on-axis imaging;\n", + " the only subtlety is matching the hyperboloid's far focus to the parabola's\n", + " *true* focus at $F_\\mathrm{p} - f_\\mathrm{p}$. The sag solvers were stable even\n", + " at sub-millimetre vertex radii.\n", + "* **The mirrors must not overlap in z.** Because the grazing flank is so steep,\n", + " each mirror is axially long, so the secondary has to sit well downstream of the\n", + " primary — otherwise part of the beam cannot reach it travelling forward. That\n", + " makes the telescope $\\sim 0.7$ m long, which is why the effective focal length\n", + " (rearward principal plane) and the grating lever arm — and hence the plate\n", + " scale ($\\approx 1.16''$/pix) and dispersion ($\\approx 0.63$ mm/Å) — come out a\n", + " bit coarser than a naively compact layout would suggest.\n", + "* **Dispersion.** A constant-pitch transmission grating in the converging beam\n", + " produces a linear, well-behaved first-order spectrum with the zero order on the\n", + " same detector.\n", + "* **Performance.** The telescope alone is diffraction-limited on-axis\n", + " ($\\ll 1$ pixel); the spectral resolution ($R \\sim 10^3$) is set mostly by the\n", + " aberration of the plane grating working in the converging beam — a real design\n", + " trade-off that a faster beam would worsen. A blazed/curved grating or a tilted\n", + " detector would recover more resolution.\n", + "\n", + "Natural next steps: fold the measured reflectivity and grating efficiency into\n", + "the ray intensities for an end-to-end throughput estimate, and study the\n", + "off-axis field aberrations that ultimately limit the usable field of view." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/optika/apertures/__init__.py b/optika/apertures/__init__.py index 084c24d..95f7e60 100644 --- a/optika/apertures/__init__.py +++ b/optika/apertures/__init__.py @@ -6,6 +6,7 @@ AbstractAperture, CircularAperture, CircularSectorAperture, + AnnularAperture, EllipticalAperture, AbstractPolygonalAperture, PolygonalAperture, @@ -22,6 +23,7 @@ "AbstractAperture", "CircularAperture", "CircularSectorAperture", + "AnnularAperture", "EllipticalAperture", "AbstractPolygonalAperture", "PolygonalAperture", diff --git a/optika/apertures/_apertures.py b/optika/apertures/_apertures.py index f30acd1..0a8b693 100644 --- a/optika/apertures/_apertures.py +++ b/optika/apertures/_apertures.py @@ -532,6 +532,191 @@ def wire(self, num: None | int = None) -> na.Cartesian3dVectorArray: return result +@dataclasses.dataclass(eq=False, repr=False) +class AnnularAperture( + AbstractAperture, +): + r""" + An annular (ring-shaped) aperture, optionally restricted to a sector. + + Defined by an inner and outer radius. By default it spans the full ring, + but the optional ``angle_start`` / ``angle_stop`` arguments restrict it to an + annular sector, analogous to :class:`CircularSectorAperture`. + Setting ``radius_inner=0`` recovers a circular sector. + + Examples + -------- + + Plot an annular aperture and an annular 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 full annulus and a 90-degree annular sector + aperture = optika.apertures.AnnularAperture( + radius_inner=30 * u.mm, + radius_outer=50 * u.mm, + angle_start=na.ScalarArray([0, 0] * u.deg, axes="aperture"), + angle_stop=na.ScalarArray([360, 90] * u.deg, axes="aperture"), + ) + + # 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(21, 21, 1), + ) + + # Compute which points are inside the aperture + where = aperture(points) + + # Plot the apertures + with astropy.visualization.quantity_support(): + fig, ax = plt.subplots(ncols=2, figsize=(8, 4), constrained_layout=True) + ax = na.ScalarArray(ax, axes="aperture") + for a in [ax[dict(aperture=0)], ax[dict(aperture=1)]]: + a.ndarray.set_aspect("equal") + aperture.plot(ax=ax, components=("x", "y"), color="black") + na.plt.scatter(points.x, points.y, ax=ax, c=where.astype(float)) + """ + + radius_inner: u.Quantity | na.AbstractScalar = 0 * u.mm + """The inner radius of the annulus.""" + + radius_outer: u.Quantity | na.AbstractScalar = 0 * u.mm + """The outer radius of the annulus.""" + + angle_start: u.Quantity | na.AbstractScalar = 0 * u.deg + r""" + The starting angle of the annular sector. + Must be between :math:`-2 \pi` and :math:`+2 \pi` radians. + """ + + angle_stop: u.Quantity | na.AbstractScalar = 360 * u.deg + r""" + The ending angle of the annular sector, + counterclockwise from :attr:`angle_start`. + The default spans the full ring. + """ + + @property + def shape(self) -> dict[str, int]: + return na.broadcast_shapes( + optika.shape(self.radius_inner), + optika.shape(self.radius_outer), + 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_inner = self.radius_inner + radius_outer = self.radius_outer + 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_inner, + radius_outer, + angle_start, + angle_stop, + active, + inverted, + position, + ) + + radius_inner = na.broadcast_to(radius_inner, shape) + radius_outer = na.broadcast_to(radius_outer, 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) + + radius = position.xy.length + mask_radius = (radius_inner <= radius) & (radius <= radius_outer) + + angle = np.arctan2(position.y, position.x) + angle_relative = (angle - angle_start) % (360 * u.deg) + mask_angle = angle_relative <= (angle_stop - angle_start) + + 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 + + # split the samples between the outer and inner arcs, reserving one + # point to close the loop + num_arc = num - 1 + num_outer = num_arc // 2 + num_inner = num_arc - num_outer + + unit_radius = na.unit(self.radius_outer) + z = 0 * unit_radius if unit_radius is not None else 0 + + azimuth_outer = na.linspace( + self.angle_start, self.angle_stop, axis="wire", num=num_outer + ) + azimuth_inner = na.linspace( + self.angle_stop, self.angle_start, axis="wire", num=num_inner + ) + + outer = na.Cartesian3dVectorArray( + x=self.radius_outer * np.cos(azimuth_outer), + y=self.radius_outer * np.sin(azimuth_outer), + z=z, + ) + inner = na.Cartesian3dVectorArray( + x=self.radius_inner * np.cos(azimuth_inner), + y=self.radius_inner * np.sin(azimuth_inner), + z=z, + ) + + result = np.concatenate([outer, inner], axis="wire") + # close the loop back to the first point + result = np.concatenate( + [result, result[dict(wire=slice(0, 1))]], + axis="wire", + ) + + if self.transformation is not None: + result = self.transformation(result) + return result + + @dataclasses.dataclass(eq=False, repr=False) class EllipticalAperture( AbstractAperture, diff --git a/optika/apertures/_apertures_test.py b/optika/apertures/_apertures_test.py index c6ae0a8..9838991 100644 --- a/optika/apertures/_apertures_test.py +++ b/optika/apertures/_apertures_test.py @@ -178,6 +178,53 @@ def test_angle_stop(self, a: optika.apertures.CircularSectorAperture): assert np.all(a.radius >= 0) +@pytest.mark.parametrize( + argnames="a", + argvalues=[ + optika.apertures.AnnularAperture( + radius_inner=radius_inner, + radius_outer=radius_outer, + angle_start=angle_start, + angle_stop=angle_stop, + samples_wire=21, + active=active, + inverted=inverted, + transformation=transformation, + kwargs_plot=kwargs_plot, + ) + for radius_inner, radius_outer in [ + (50 * u.mm, 100 * u.mm), + (na.linspace(20, 40, axis="radius", num=4) * u.mm, 100 * u.mm), + ] + for angle_start, angle_stop in [ + (0 * u.deg, 360 * u.deg), + (30 * u.deg, 120 * u.deg), + ] + 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 TestAnnularAperture( + AbstractTestAbstractAperture, +): + def test_radius_inner(self, a: optika.apertures.AnnularAperture): + assert isinstance(a.radius_inner, (float, u.Quantity, na.AbstractScalar)) + assert np.all(a.radius_inner >= 0) + + def test_radius_outer(self, a: optika.apertures.AnnularAperture): + assert isinstance(a.radius_outer, (float, u.Quantity, na.AbstractScalar)) + assert np.all(a.radius_outer >= a.radius_inner) + + def test_angle_start(self, a: optika.apertures.AnnularAperture): + assert na.unit_normalized(a.angle_start).is_equivalent(u.deg) + + def test_angle_stop(self, a: optika.apertures.AnnularAperture): + assert na.unit_normalized(a.angle_stop).is_equivalent(u.deg) + assert np.all(a.angle_stop >= a.angle_start) + + class AbstractTestAbstractPolygonalAperture( AbstractTestAbstractAperture, ):