Introduction — Modeling a Spectroscopic Instrument¶
The Central Question¶
Given:
- a source at sky position $(\alpha, \delta)$
- emitting light at wavelength $\lambda$
Predict: the pixel position $(x, y)$ where that photon ends on the detector.
$$ (x,\, y) = f\!\left(\alpha,\, \delta,\, \lambda\right) $$
This function is called the dispersion model.
The Instrument Pipeline¶
Each element performs one physical transformation. The output of one is the input of the next.
sky angle (α, δ) [rad]
↓
Telescope angle → position [mm]
↓
Collimator position → angle [rad]
↓
Grism angle → dispersed angle (λ-dep.) [rad]
↓
Camera angle → position [mm]
↓
Detector position → pixel
↓
pixel (x, y)
This notebook builds that pipeline step by step, starting from a single lens equation and ending with a full 2D optical chain.
1. Setup and Imports¶
We will use:
- NumPy for vectorized array math (applying one formula to thousands of rays or wavelengths at once)
- matplotlib for plotting curves (refractive index, dispersion, spectral traces)
import numpy as np
import matplotlib.pyplot as plt
2. Ray Representation and the Thin Lens Equation¶
Describing a Light Ray¶
At any cross-section along the optical axis, a ray is described by two numbers:
- $x$ — height above the optical axis [m]
- $\theta$ — angle to the optical axis [rad]
As the ray passes through each element, $(x, \theta)$ changes.
Small angle approximation:
$$ \sin\theta \approx \theta \qquad \tan\theta \approx \theta $$
This makes every transformation linear. Geometric distortions beyond paraxial will be handled separately.
The Thin Lens Equation¶
A thin lens of focal length $f$ deflects the angle of a ray:
$$ \theta' = \theta - \frac{x}{f} \qquad x' = x $$
The position $x$ is unchanged at the lens. Only the angle changes.
Physical meaning:
- Ray through the center ($x = 0$): not deflected.
- Ray at height $x$: bent toward the focal point by $x/f$.
- All parallel rays ($\theta = 0$): converge to the focal point.
This one equation drives most of the instrument model. Collimator, Camera, and Telescope are all thin lenses — each used in a specific configuration of $(x, \theta)$.
The Paraxial Hypothesis¶
Under the paraxial approximation, every element applies a linear transformation: $\text{output} = a \cdot \text{input}$ for some real coefficient $a$. This has a crucial consequence: the two transverse axes ($x$ and $y$) are decoupled.
- Dispersion direction (x): $\alpha \to \text{telescope} \to \text{collimator} \to \text{grism} \to \text{camera} \to x_\text{pix}$. The grism adds $\delta(\lambda)$ along this axis only: $x_\text{pix} = f(\alpha, \lambda)$.
- Cross-dispersion direction (y): $\delta \to \text{telescope} \to \text{collimator} \to \text{camera} \to y_\text{pix}$. No grism effect — pure magnification: $y_\text{pix} = g(\delta)$.
Sections 3–10 focus on the dispersion direction (x axis) only. We trace how one sky angle $\alpha$ maps to one pixel $x$ as a function of wavelength. This is the minimal model to understand dispersion — the y axis adds nothing new yet (it comes back in Section 11).
# One ray
x = 0.5 # height above optical axis [m]
theta = 0.0 # angle [rad] — this ray is parallel to the axis
# One lens
f = 2.0 # focal length [m]
# Apply the thin lens equation
theta_out = theta - x / f # angle is deflected
x_out = x # position is unchanged at the lens
print("theta_out =", theta_out) # -> -0.25
print("x_out =", x_out) # -> 0.5
theta_out = -0.25 x_out = 0.5
Problem: if we have 3 lenses, we copy these 2 lines 3 times. A mistake in the formula must be corrected in 3 places.
Step 2 — Write a Function¶
Python concept: function (write once, call anywhere)
A function has a name, receives inputs (parameters), and produces outputs. Write it once, fix it once.
def apply_lens(x, theta, f):
"""
Apply a thin lens of focal length f to a ray.
Parameters
----------
x : float — height above the optical axis [m]
theta : float — angle to the optical axis [rad]
f : float — focal length [m]
Returns
-------
x_out, theta_out : float, float
"""
return x, theta - x / f
# Same result, now reusable:
x_out, theta_out = apply_lens(x=0.5, theta=0.0, f=2.0)
print(theta_out) # -> -0.25 rad
-0.25
NumPy — Arrays and Vectorization¶
Python concept: numerical arrays
NumPy is the standard Python library for numerical computing. Its core object is the ndarray — a typed, multi-dimensional array on which arithmetic operates without any loop. This is called vectorization: one call processes thousands of rays or wavelengths at once, as fast as C.
Documentation: numpy.org/doc/stable — in particular NumPy quickstart and Broadcasting.
x = np.array([0.0, 0.5, 1.0]) # 1D array of 3 floats
print(x * 2) # -> [0. 1. 2. ] element-wise
print(x - x / 2) # -> [0. 0.25 0.5] element-wise
print(np.sqrt(x)) # -> [0. 0.707 1. ]
[0. 1. 2.] [0. 0.25 0.5 ] [0. 0.70710678 1. ]
Step 3 — Test the Function¶
Exercise: run this in your notebook
Apply apply_lens to three parallel rays at once, using arrays, and verify they converge to the focal point.
# Three parallel rays at different heights
x_rays = np.array([0.0, 0.5, 1.0])
theta_rays = np.array([0.0, 0.0, 0.0]) # all parallel to axis
f = 2.0
x_out, theta_out = apply_lens(x_rays, theta_rays, f)
print("Output angles:", theta_out)
# Expected: [ 0. -0.25 -0.5 ]
# A ray farther from the axis is bent more.
# Verify: propagate to the focal plane and check convergence
x_focal = x_out + f * theta_out
print("Positions at focal plane:", x_focal)
# Expected: all close to 0.0 — they converge to the focal point
Output angles: [ 0. -0.25 -0.5 ] Positions at focal plane: [0. 0. 0.]
4. Collimator and Camera Functions¶
Collimator — Position → Angle¶
Role: placed at distance $f_\text{coll}$ from the telescope focal plane, it converts each focal-plane point into a parallel beam.
A source at position $x$ in the focal plane produces a beam at angle:
$$ \theta = -\frac{x}{f_\text{coll}} $$
This is the thin lens formula applied to a point source at the focal plane: $\theta_\text{in} = 0$, so $\theta_\text{out} = 0 - x/f = -x/f$.
Camera — Angle → Position¶
Role: re-focuses a parallel beam onto the detector. A collimated beam at angle $\theta$ converges to a point at:
$$ x = -f_\text{cam}\,\theta $$
Collimator + Camera together:
$$ x_\text{det} = \frac{f_\text{cam}}{f_\text{coll}}\,x_\text{focal} $$
The image is magnified by $\gamma = f_\text{cam}/f_\text{coll}$.
The Grism acts between them. With the grism: $\theta \to \theta + \delta(\lambda) \to \text{camera} \to$ dispersed image. The extra $\delta(\lambda)$ shifts the image by $-f_\text{cam}\cdot\delta(\lambda)$, a different amount for each wavelength.
Named Functions¶
Python concept: one function, one clear responsibility
Each function takes one physical quantity as input and returns one. Named, single-purpose functions make the chain explicit.
def collimator(x_focal, f_coll):
"""
Collimator: focal-plane position → beam direction.
Parameters
----------
x_focal : float — position [m]
f_coll : float — focal length [m]
Returns
-------
theta : float — beam direction [rad]
"""
return -x_focal / f_coll
def camera(theta, f_cam):
"""
Camera: beam direction → detector position.
Parameters
----------
theta : float — beam direction [rad]
f_cam : float — focal length [m]
Returns
-------
x : float — detector position [m]
"""
return -f_cam * theta
f_coll = 0.5 # [m]
f_cam = 0.3 # [m]
# A source at x = 1 mm in the focal plane
x_focal = 1e-3 # [m]
# Step 1: collimator
theta = collimator(x_focal, f_coll)
print(f"beam angle : {theta*1e3:.2f} mrad") # -> -2.00 mrad
# Step 2: no grism (delta = 0)
theta_disp = theta + 0.0
# Step 3: camera
x_det = camera(theta_disp, f_cam)
print(f"detector x : {x_det*1e3:.2f} mm") # -> 0.60 mm
# Check: magnification gamma = f_cam / f_coll = 0.3/0.5 = 0.6
gamma = f_cam / f_coll
print(f"gamma : {gamma:.2f}")
print(f"expected x : {x_focal * gamma * 1e3:.2f} mm")
beam angle : -2.00 mrad detector x : 0.60 mm gamma : 0.60 expected x : 0.60 mm
5. Telescope and Class-Based Design¶
Telescope¶
Role: convert a sky angle into a position in the focal plane. The telescope is pointed at a centre $\alpha_0$ (the optical axis). A source at sky angle $\alpha$ produces an angular offset $\alpha - \alpha_0$:
$$ x_\text{focal} = -f_\text{tel}\,(\alpha - \alpha_0) $$
This is the same formula as the Camera: $x = -f \cdot \theta$ with $\theta = \alpha - \alpha_0$.
Euclid: $f_\text{tel} = 24.5$ m → 1 arcmin offset ≈ 7.1 mm in the focal plane.
The Camera Class¶
Python concept: a class bundles a parameter and a method
__init__ runs when the object is created — stores parameters. self is the object itself — it carries its own data. forward() is a method: a function that belongs to the class.
class Camera:
"""
Camera: beam angle [rad] → focal-plane position [m].
x = -f * theta
Parameters
----------
f : float — focal length [m]
"""
def __init__(self, f):
self.f = f # store focal length as an attribute
def forward(self, theta):
"""Beam angle [rad] → position [m]."""
return -self.f * theta
# Create a Camera object (stores f once):
cam = Camera(f=0.3)
x = cam.forward(theta=0.01) # 10 mrad -> -3 mm
print(f"x = {x*1e3:.1f} mm") # -> -3.0 mm
x = -3.0 mm
Telescope — Inheriting from Camera¶
Python concept: inheritance + extending with a new parameter
super().__init__(f) calls Camera.__init__. Inheritance lets us extend a class without rewriting it. A Telescope is a Camera (same formula), so it inherits from it and just adds the pointing offset.
class Telescope(Camera):
"""
Telescope: sky angle [rad] → focal-plane position [m].
Inherits Camera (x = -f * angle).
Adds a pointing centre alpha0: only the offset (alpha - alpha0) matters.
Parameters
----------
f : float — focal length [m]
alpha0 : float — pointing centre [rad] (default: 0)
"""
def __init__(self, f, alpha0=0.0):
super().__init__(f) # call Camera.__init__(f) to set self.f
self.alpha0 = alpha0 # store pointing centre
def forward(self, alpha):
"""Sky angle [rad] → focal-plane position [m]."""
return super().forward(alpha - self.alpha0) # subtract offset, then focus
# On-axis pointing (default):
tel = Telescope(f=1.2)
print(tel.forward(0.001)) # -> -1.2e-3 m
# Off-axis pointing (centre at 5 mrad):
tel = Telescope(f=1.2, alpha0=5e-3)
print(tel.forward(5e-3)) # -> 0.0 m (pointing centre lands on axis)
print(tel.forward(6e-3)) # -> -1.2e-3 m (1 mrad offset from pointing)
-0.0012 -0.0 -0.0012
Collimator as a Class¶
Every element has the same interface: .forward(input)
Reading the instrument chain becomes as natural as reading the pipeline diagram.
class Collimator:
"""
Convert a focal-plane position into a beam direction.
theta = -x / f
Parameters
----------
f : float — focal length [m]
"""
def __init__(self, f):
self.f = f
def forward(self, x):
"""Position [m] → beam direction [rad]."""
return -x / self.f
# Test the tel -> coll chain:
tel = Telescope(f=1.2)
coll = Collimator(f=0.5)
alpha = 1e-3 # sky angle [rad]
x_focal = tel.forward(alpha) # sky angle → focal position
theta = coll.forward(x_focal) # focal position → beam direction
print(f"focal-plane : {x_focal*1e3:.2f} mm")
print(f"direction : {theta*1e6:.1f} µrad")
# expected: -1.20 mm, -2400.0 µrad
focal-plane : -1.20 mm direction : 2400.0 µrad
6. Refractive Index (Material Class)¶
Light in a Medium¶
Light slows down in a transparent medium: $v = c/n$. At each interface, the ray bends (Snell's law).
Dispersion: $n = n(\lambda)$. $n$ depends on wavelength. Blue (short $\lambda$) slows more than red (long $\lambda$) → different wavelengths exit at different angles. This is the physical origin of spectral dispersion.
Cauchy model (1836):
$$ n(\lambda) = n_0 + \frac{k}{\lambda^2} $$
$\lambda$ in µm; $n_0$ and $k$ depend on the glass type.
class Material:
"""
Optical glass: Cauchy dispersion model.
Parameters
----------
n0 : float — base refractive index
k : float — dispersion coefficient [µm²]
"""
def __init__(self, n0, k):
self.n0 = n0
self.k = k
def refractive_index(self, wavelength):
"""Compute n(λ). wavelength in µm."""
return self.n0 + self.k / wavelength**2
# Test:
glass = Material(n0=1.44, k=0.004)
print(glass.refractive_index(1.0)) # -> 1.444
print(glass.refractive_index(2.0)) # -> 1.441 (n decreases with λ)
1.444 1.4409999999999998
Visualization¶
matplotlib is the standard Python plotting library. The pyplot sub-module gives a MATLAB-like interface: create a figure, add plots, label axes, then show or save.
n decreases with $\lambda$ (normal dispersion). The variation from 0.9 to 2.0 µm is small (~0.003) but sufficient to spread the spectrum across many pixels on the detector.
glass = Material(n0=1.44, k=0.004)
wavelength = np.linspace(0.9, 2.0, 200) # NISP NIR range [µm]
n = glass.refractive_index(wavelength)
plt.figure(figsize=(6, 3))
plt.plot(wavelength, n, color='steelblue')
plt.xlabel("Wavelength [µm]")
plt.ylabel("n(λ)")
plt.title("Refractive index — Cauchy model")
plt.tight_layout()
plt.show()
7. The Prism (Composition)¶
Snell's law (paraxial form): $n_1\,\theta_1 = n_2\,\theta_2$
Net deviation of a thin prism (apex angle $A$):
$$ \delta_\text{prism}(\lambda) = A\,\bigl(n(\lambda) - 1\bigr) $$
Because $n$ depends on $\lambda$, different wavelengths are deflected by different amounts — the dispersive action of the prism.
Python concept: Prism has-a Material (composition)
self.material = material — the Prism stores a reference to a Material object. This is composition: one object contains another.
class Prism:
"""
Thin prism: angular deviation δ(λ) = A * (n(λ) - 1).
Parameters
----------
material : Material — the glass (provides n(λ))
A : float — apex angle [rad]
"""
def __init__(self, material, A):
self.material = material # Prism stores a Material object
self.A = A
def deviation(self, wavelength):
"""Angular deviation [rad]. wavelength in µm."""
n = self.material.refractive_index(wavelength)
return self.A * (n - 1)
# Test:
glass = Material(n0=1.44, k=0.004)
prism = Prism(material=glass, A=0.05) # apex angle 0.05 rad ~ 2.9 deg
print(np.degrees(prism.deviation(1.0))) # deviation at 1 µm [deg]
1.2719663051904275
8. The Diffraction Grating¶
Grating equation:
$$ d\,(\sin\theta_m - \sin\theta_i) = m\,\lambda $$
- $d$ = groove spacing [µm]
- $m$ = diffraction order
- $\theta_i$ = incident angle, $\theta_m$ = diffracted angle
Paraxial, normal incidence ($\theta_i \approx 0$):
$$ \delta_\text{grating}(\lambda) \approx \frac{m\,\lambda}{d} $$
Unlike the prism, grating dispersion is linear in $\lambda$.
Datasheets usually quote the groove density $\rho$ (grooves/mm) rather than the spacing: $d\,[\mu\text{m}] = 1000/\rho\,[\text{mm}^{-1}]$.
class Grating:
"""
Diffraction grating: angular deviation δ(λ) = m * λ / d.
Parameters
----------
m : int — diffraction order (can be negative)
rho : float — groove density [grooves/mm]
"""
def __init__(self, m, rho):
self.m = m
self.rho = rho # groove density [grooves/mm]
self.d = 1e3 / rho # groove spacing [µm]
def deviation(self, wavelength):
"""Angular deviation [rad]. wavelength in µm."""
return self.m * wavelength / self.d
# Test (Euclid/NISP grating: ~13.75 grooves/mm):
grating = Grating(m=1, rho=13.75)
d1 = np.degrees(grating.deviation(1.0))
d2 = np.degrees(grating.deviation(2.0))
print(f"groove spacing d = {grating.d:.2f} µm")
print(f"δ(1.0 µm) = {d1:.3f}°, δ(2.0 µm) = {d2:.3f}°")
# Grating dispersion is linear: δ(2µm) = 2 × δ(1µm)
groove spacing d = 72.73 µm δ(1.0 µm) = 0.788°, δ(2.0 µm) = 1.576°
9. The Grism (Composing Prism + Grating)¶
A grism = a grating bonded to a prism.
The prism is chosen so that at the central wavelength $\lambda_0$, prism and grating deviations cancel:
$$ \delta_\text{prism}(\lambda_0) - \delta_\text{grating}(\lambda_0) = 0 $$
→ the beam passes straight through at $\lambda_0$. Other $\lambda$ are deflected proportionally to $(\lambda - \lambda_0)$.
Total deviation:
$$ \delta_\text{grism}(\lambda) = \delta_\text{prism}(\lambda) - \delta_\text{grating}(\lambda) $$
The minus sign is physical, not a free choice: the grating is ruled on the prism's exit face, and both Prism.deviation and Grating.deviation are defined relative to the optical axis. Refraction (prism) and diffraction (grating) bend the beam in opposite rotational senses relative to that shared axis, so the grating term must be subtracted. With this combination, a grism cancels at its design wavelength using its ordinary, positive diffraction order (e.g. $m=+1$, matching real grism hardware) — there is no need to flip the sign of $m$ to make the two contributions oppose each other.
Python concept: Grism has-a Prism and has-a Grating
Grism does not inherit from Prism or Grating. A Grism is not a special kind of prism — it contains both. This is composition.
Inheritance vs Composition:
- Inheritance — "is-a":
class Telescope(Camera). A Telescope IS a Camera (same formula, different name) → inherit. - Composition — "has-a":
Grism(prism, grating). A Grism HAS a Prism and HAS a Grating; its behavior combines two components → store both, combine their deviations.
Rule: use inheritance when the new class specializes an existing one. Use composition when the new class is built from existing components.
class Grism:
"""
Grism: prism + grating.
The grating is ruled on the prism's exit face, so — with both angles
measured relative to the optical axis — refraction (prism) and
diffraction (grating) bend the beam in opposite rotational senses.
The grating term therefore subtracts from the prism term.
Total deviation = prism deviation - grating deviation.
Parameters
----------
prism : Prism
grating : Grating
"""
def __init__(self, prism, grating):
self.prism = prism # has-a Prism
self.grating = grating # has-a Grating
def deviation(self, wavelength):
"""Total angular deviation [rad]. wavelength in µm."""
return (self.prism.deviation(wavelength)
- self.grating.deviation(wavelength))
def forward(self, theta, wavelength):
"""
Apply the grism: beam angle → dispersed angle.
theta : float or array — beam angle [rad]
wavelength : float or array — wavelength [µm]
"""
return theta + self.deviation(wavelength)
Grism — Dispersion Curve¶
The wavelength where $\delta = 0$ is the blaze wavelength $\lambda_0$.
glass = Material(n0=1.44, k=0.004)
prism = Prism(material=glass, A=np.radians(2.145)) # Euclid/NISP prism apex angle
grating = Grating(m=1, rho=13.75) # ordinary (+1) order; Grism.deviation subtracts it from the prism term
grism = Grism(prism=prism, grating=grating)
wavelength = np.linspace(0.9, 2.0, 200)
delta = np.degrees(grism.deviation(wavelength))
plt.figure(figsize=(6, 3))
plt.plot(wavelength, delta, color='darkorange')
plt.axhline(0, color='gray', lw=0.8, linestyle='--')
plt.xlabel("Wavelength [µm]")
plt.ylabel("Deviation [deg]")
plt.title("Grism angular deviation vs wavelength")
plt.tight_layout()
plt.show()
class Detector:
"""
Detector: focal-plane position [m] → pixel coordinate.
pixel = (x + x0) / pixel_size
Parameters
----------
pixel_size : float — pixel size [m/pix]
x0 : float — detector center offset [m]
"""
def __init__(self, pixel_size, x0=0.0):
self.pixel_size = pixel_size
self.x0 = x0
def forward(self, x):
"""Physical position [m] → pixel coordinate."""
return (x + self.x0) / self.pixel_size
# Test (NISP: 18 µm pixels):
det = Detector(pixel_size=18e-6, x0=0.0)
pix = det.forward(x=1.8e-3) # 1.8 mm -> 100.0 pixels
print(f"pixel: {pix:.1f}") # -> 100.0
pixel: 100.0
The Complete Chain¶
Every element now has the same interface: .forward(input). The chain reads almost like the pipeline diagram.
def trace_instrument(alpha, wavelength, tel, coll, grism, cam, det):
"""
Forward model: sky field angle + wavelength → detector pixel.
Parameters
----------
alpha : float or array — sky field angle [rad]
wavelength : float or array — wavelength [µm]
tel : Telescope — sky angle → focal position [m]
coll : Collimator — position → beam angle [rad]
grism : Grism — angle → dispersed angle (λ-dep.)
cam : Camera — angle → detector position [m]
det : Detector — position → pixel
Returns
-------
pixel : float or array
"""
x = tel.forward(alpha) # sky angle → focal plane [m]
theta = coll.forward(x) # position → beam angle [rad]
theta = grism.forward(theta, wavelength) # dispersed angle (λ-dep.)
x = cam.forward(theta) # angle → detector position [m]
return det.forward(x) # position → pixel
Full Instrument — Build and Plot the Dispersion Curve¶
We assemble the instrument using approximate Euclid/NISP parameters:
| Element | Parameter | Value |
|---|---|---|
| Telescope | focal length | 24.5 m |
| Collimator | focal length | 2.0 m |
| Camera | focal length | 1.0 m |
| Prism | apex angle $A$ | 2.145° |
| Grating | groove density $\rho$ | 13.75 grooves/mm |
| Detector | pixel size | 18 µm |
| — | operating wavelength range | 1.2 – 1.9 µm |
The grating is used in its ordinary order $m=+1$ (see Section 9) — Grism.deviation subtracts the grating term from the prism term, so the two naturally cancel near the design wavelength (checked below). The glass is the same Cauchy-fit material introduced in Section 6, used here as an approximation of the real fused-silica prism/grating.
# Build the instrument (Euclid/NISP-like parameters)
tel = Telescope(f=24.5)
coll = Collimator(f=2.0)
glass = Material(n0=1.44, k=0.004)
prism = Prism(material=glass, A=np.radians(2.145))
grating = Grating(m=1, rho=13.75)
grism = Grism(prism, grating)
cam = Camera(f=1.0)
det = Detector(pixel_size=18e-6, x0=0.0)
# Source at field center (alpha = 0), vary wavelength
alpha = 0.0
wavelength = np.linspace(1.2, 1.9, 200) # [µm]
pixel = trace_instrument(alpha, wavelength, tel, coll, grism, cam, det)
plt.figure(figsize=(6, 3))
plt.plot(wavelength, pixel)
plt.xlabel("Wavelength [µm]")
plt.ylabel("Detector position [pix]")
plt.title("Euclid/NISP dispersion curve — α = 0 (field center)")
plt.tight_layout()
plt.show()
Verification — Undeflected Wavelength¶
The wavelength where the grism's total deviation crosses zero (its ordinary order $m=+1$ subtracted from the prism deviation via Grism.deviation) should be close to 1.2 µm.
lam_fine = np.linspace(1.0, 1.4, 4000)
dev_fine = grism.deviation(lam_fine)
lambda0 = lam_fine[np.argmin(np.abs(dev_fine))]
print(f"undeflected wavelength ~ {lambda0:.3f} µm (target: ~1.2 µm)")
undeflected wavelength ~ 1.205 µm (target: ~1.2 µm)
Verification — Spectral Resolution¶
The local spectral dispersion $d\lambda/d\text{pixel}$ should be close to 1.372 nm/pixel.
dlambda_dpixel = np.gradient(wavelength, pixel) * 1000 # µm -> nm per pixel
resolution = np.mean(np.abs(dlambda_dpixel))
print(f"spectral resolution ~ {resolution:.3f} nm/pix (target: ~1.372 nm/pix)")
spectral resolution ~ 1.301 nm/pix (target: ~1.372 nm/pix)
What the Model Gives Us¶
We have a 1D forward model: given sky angle $\alpha$ and $\lambda$ → predict detector pixel $x$.
$$ x_\text{pix} = f\!\left(\alpha,\, \lambda\right) $$
The same chain applied to $\delta$ gives $y$ independently. Section 11 shows how to combine them into a single 2D model.
This simplified model neglects geometric distortions and lateral color (present in real instruments as higher-order terms). Such residuals are typically corrected later by machine learning. The verifications above already show a small gap between our simplified Cauchy-glass model and the real instrument's quoted spectral resolution — exactly the kind of residual a hybrid physical + ML model is designed to correct.
11. Extension to 2D (Vectors, Rotation, Pointing)¶
Where the 1D Model Breaks¶
In Section 2 we argued that $x$ and $y$ are decoupled under the paraxial approximation — but this relied on the grism being perfectly aligned with the x-axis. In reality, it is not. The detector position depends on both $\alpha$ and $\delta$:
$$ x_\text{pix} = f(\alpha, \delta, \lambda) \qquad y_\text{pix} = g(\alpha, \delta, \lambda) $$
The 1D model of Sections 3–10 is a pedagogical simplification. It is a good first step — but a real calibration model must track both axes together.
A Cleaner Design — 2D Vectors¶
Each quantity becomes a NumPy array of shape (2,) or (2, N)
Rather than writing two methods per element (one for x, one for y — error-prone and hard to keep in sync), we redefine each element to operate on a 2D vector directly. One method per element, one call per step. NumPy broadcasts over the wavelength axis automatically.
def rotation_matrix(angle):
"""2D rotation matrix for a given angle [rad]."""
c, s = np.cos(angle), np.sin(angle)
return np.array([[c, -s],
[s, c]])
class Camera2D:
"""
Camera: beam direction [rad] → focal-plane position [m], in 2D.
x = -f * theta
Parameters
----------
f : float — focal length [m]
"""
def __init__(self, f):
self.f = f
def forward(self, theta):
"""theta : (2,) or (2, N) array → position, same shape."""
return -self.f * theta
class Collimator2D:
"""
Collimator: focal-plane position [m] → beam direction [rad], in 2D.
theta = -x / f
Parameters
----------
f : float — focal length [m]
"""
def __init__(self, f):
self.f = f
def forward(self, x):
"""x : (2,) or (2, N) array → direction, same shape."""
return -x / self.f
class Detector2D:
"""
Detector: focal-plane position [m] → pixel coordinates, in 2D.
pixel = (x + offset) / pixel_size
Parameters
----------
pixel_size : float — pixel size [m/pix]
offset : (2,) array — detector center offset [m]
"""
def __init__(self, pixel_size, offset=None):
self.pixel_size = pixel_size
self.offset = np.zeros(2) if offset is None else np.asarray(offset, dtype=float)
def forward(self, x):
"""x : (2,) or (2, N) array → pixel coordinates, same shape."""
x = np.asarray(x)
offset = self.offset if x.ndim == 1 else self.offset[:, None]
return (x + offset) / self.pixel_size
Telescope2D — Pointing and Field Rotation¶
Sky position $[\alpha, \delta]$ → focal-plane position $[x, y]$.
The pointing is the sky position $[\alpha_0, \delta_0]$ of the field centre. Additionally the focal plane may be rotated relative to the equatorial axes by a roll angle $\psi$.
Telescope2D inherits from Camera2D (same formula $x = -f\,\theta$), and adds the pointing offset and roll rotation before calling the parent's forward.
class Telescope2D(Camera2D):
"""
2D telescope: sky position → focal-plane position.
Parameters
----------
f : float — focal length [m]
sky0 : (2,) array — pointing centre [alpha0, delta0] [rad]
angle : float — field rotation / roll angle [rad]
"""
def __init__(self, f, sky0=None, angle=0.0):
super().__init__(f)
self.sky0 = np.zeros(2) if sky0 is None else np.asarray(sky0, dtype=float)
self.angle = angle
self.R = rotation_matrix(angle)
def forward(self, sky):
"""
sky : (2,) or (2, N) array [alpha, delta] [rad]
Returns : same shape — focal-plane position [m]
"""
sky = np.asarray(sky)
sky0 = self.sky0 if sky.ndim == 1 else self.sky0[:, None]
offset = sky - sky0 # angular offset from pointing
return super().forward(self.R @ offset) # rotate, then focus
# No rotation - 1D result recovered on each axis:
tel2d = Telescope2D(f=24.5, sky0=np.array([0., 0.]))
print(tel2d.forward(np.array([1e-4, 2e-4]))) # -> [-2.45e-3, -4.90e-3] m
# 45 deg roll - the alpha-offset now contributes to both x and y:
tel2d_45 = Telescope2D(f=24.5, sky0=np.zeros(2), angle=np.radians(45))
print(tel2d_45.forward(np.array([1e-4, 0.])))
[-0.00245 -0.0049 ] [-0.00173241 -0.00173241]
Grism2D — Rotation + Dispersion¶
The prism tilt is a 2D rotation; dispersion adds to one axis (the grating grooves run along the other). The beam is rotated to align with the prism apex, dispersion is added along one axis, then rotated back.
class Grism2D:
"""
Grism (2D): prism + grating, with a tilt angle coupling x and y.
As in the 1D Grism (Section 9), the grating term subtracts from the
prism term: refraction and diffraction bend the beam in opposite
rotational senses when both are measured relative to the optical axis.
Parameters
----------
prism : Prism
grating : Grating
tilt : float — prism/grating tilt relative to the x-axis [rad]
"""
def __init__(self, prism, grating, tilt=0.0):
self.prism = prism
self.grating = grating
self.R = rotation_matrix(tilt)
def deviation(self, wavelength):
"""Total angular deviation [rad] (scalar or array). wavelength in µm."""
return (self.prism.deviation(wavelength)
- self.grating.deviation(wavelength))
def forward(self, angle, wavelength):
"""
angle : (2,) or (2, N) array — beam direction [rad]
wavelength : float or (N,) array — wavelength [µm]
Returns : (2,) or (2, N) array — dispersed beam direction [rad]
"""
angle = np.asarray(angle, dtype=float)
rotated = self.R @ angle # align with prism apex
dev = np.asarray(self.deviation(wavelength), dtype=float)
if rotated.ndim == 1 and dev.ndim == 1 and dev.shape[0] > 1:
# single ray, many wavelengths -> broadcast to (2, N)
rotated = np.broadcast_to(rotated[:, None], (2, dev.shape[0])).copy()
else:
rotated = rotated.copy()
rotated[1] = rotated[1] + dev # dispersion on the y' axis
return self.R.T @ rotated # rotate back
Full 2D Chain¶
One function chains all elements, mirroring the 1D trace_instrument but operating on 2D vectors throughout.
def trace_instrument_2d(sky, wavelength, tel, coll, grism, cam, det):
"""
2D forward model: sky position + wavelength → detector pixel (x, y).
Parameters
----------
sky : (2,) array — [alpha, delta] [rad]
wavelength : float or (N,) array — [µm]
Returns
-------
pixel : (2,) or (2, N) array — [pixel_x, pixel_y]
"""
pos_foc = tel.forward(sky) # sky → focal plane [m]
angle_col = coll.forward(pos_foc) # position → beam direction [rad]
angle_gr = grism.forward(angle_col, wavelength) # dispersed angle (λ-dep.)
pos_cam = cam.forward(angle_gr) # angle → detector position [m]
return det.forward(pos_cam) # position → pixel
Full 2D Chain — Test¶
Build the 2D instrument and plot spectral traces for several sources at different $(\alpha, \delta)$ offsets. The $\delta$-offset source traces a tilted line — not a horizontal one. This is the signature of the prism tilt coupling the two axes.
# Instrument parameters (same physical values as Section 10, in 2D)
tel2d = Telescope2D(f=24.5)
coll2d = Collimator2D(f=2.0)
grism2d = Grism2D(prism, grating, tilt=np.radians(2))
cam2d = Camera2D(f=1.0)
det2d = Detector2D(pixel_size=18e-6)
wl = np.linspace(1.2, 1.9, 200) # [µm]
# A few sources: (alpha, delta) in rad
sources = {
"on-axis" : np.array([0., 0. ]),
"alpha offset" : np.array([1e-4, 0. ]),
"delta offset" : np.array([0., 1e-4 ]),
"both offsets" : np.array([1e-4, 1e-4 ]),
}
fig, ax = plt.subplots(figsize=(6, 5))
for label, sky in sources.items():
pix = trace_instrument_2d(sky, wl, tel2d, coll2d, grism2d, cam2d, det2d)
ax.plot(pix[0], pix[1], label=label) # pix[0] = x, pix[1] = y
ax.set_xlabel("pixel x")
ax.set_ylabel("pixel y")
ax.legend()
ax.set_title("Spectral traces — 2D model")
plt.tight_layout()
plt.show()
12. Exercises¶
The exercises below reinforce the concepts covered above. Add code cells as needed and run them to check your results.
Exercise 1 — Implement the Grating class from scratch¶
Without looking back at Section 8, write your own Grating class with:
__init__(self, m, rho)storing the diffraction order and groove density [grooves/mm] (compute the spacingd = 1000 / rhointernally)deviation(self, wavelength)returning $\delta(\lambda) = m\lambda/d$
Verify that deviation(2.0) is exactly twice deviation(1.0) (linearity).
# Exercise 1: implement your own Grating class here
# class Grating:
# ...
# then test it:
# g = Grating(m=1, rho=13.75)
# assert np.isclose(g.deviation(2.0), 2 * g.deviation(1.0))
Exercise 2 — Build the full 2D chain from the 1D version¶
Starting from the 1D classes of Section 10 (Telescope, Collimator, Camera, Detector, Grism), extend each one so that forward accepts a (2,) (or (2, N)) array and returns the same shape — this is exactly what Telescope2D, Collimator2D, Camera2D, Detector2D, and Grism2D do above. Confirm you understand why each method changes (or doesn't) compared to its 1D counterpart.
Exercise 3 — A 3×3 grid of sources¶
Define a 3×3 grid of sources: $\alpha \in \{-1, 0, +1\}\cdot10^{-4}$ rad $\times$ $\delta \in \{-1, 0, +1\}\cdot10^{-4}$ rad (9 sources total). Plot all 9 spectral traces on the same detector-pixel plot, as was done for 4 sources above.
# Exercise 3: build and plot a 3x3 grid of sources
offsets = np.array([-1e-4, 0., 1e-4])
fig, ax = plt.subplots(figsize=(6, 5))
for a in offsets:
for d in offsets:
sky = np.array([a, d])
pix = trace_instrument_2d(sky, wl, tel2d, coll2d, grism2d, cam2d, det2d)
ax.plot(pix[0], pix[1], label=f"α={a:.0e}, δ={d:.0e}")
ax.set_xlabel("pixel x")
ax.set_ylabel("pixel y")
ax.set_title("Spectral traces — 3x3 source grid")
ax.legend(fontsize=6, ncol=2)
plt.tight_layout()
plt.show()
Exercise 4 — Explore the effect of the tilt angle¶
Re-run the 3×3 grid plot for a few different values of tilt passed to Grism2D (e.g. 0°, 2°, 10°, 45°). How does the tilt angle affect:
- the slope of each spectral trace?
- the coupling between $\alpha$ and $\delta$ offsets (i.e. how much a $\delta$-only offset shifts the trace in $x$)?
# Exercise 4: explore the effect of the grism tilt angle
tilt_values_deg = [0, 2, 10, 45]
fig, axes = plt.subplots(1, len(tilt_values_deg), figsize=(4 * len(tilt_values_deg), 4), sharey=True)
for ax, tilt_deg in zip(axes, tilt_values_deg):
grism_test = Grism2D(prism, grating, tilt=np.radians(tilt_deg))
for a in offsets:
for d in offsets:
sky = np.array([a, d])
pix = trace_instrument_2d(sky, wl, tel2d, coll2d, grism_test, cam2d, det2d)
ax.plot(pix[0], pix[1])
ax.set_title(f"tilt = {tilt_deg}°")
ax.set_xlabel("pixel x")
axes[0].set_ylabel("pixel y")
plt.tight_layout()
plt.show()
13. From Notebook to Library¶
The vector-capable (2D) classes above — Telescope2D, Collimator2D, Grism2D,
Camera2D, Detector2D, trace_instrument_2d — are the real, general model;
the scalar 1D classes of Sections 3–10 were only a stepping stone to derive the
physics. They have been promoted into the dispcraft library under
dispcraft/optics/ (plain names, no 2D suffix — there's no 1D sibling left
to disambiguate from), with their physical parameters moved out of the
notebook into models/stage1_instrument.toml.
This section rebuilds the same instrument from that config through
dispcraft.optics, and re-runs the Section 11 spectral-trace plot and the
Section 10 verification checks to confirm the promoted code reproduces the
same physics. No new logic is introduced here — this is a thin driver over
the library.
import tomllib
from pathlib import Path
from dispcraft.optics import build_instrument_from_config
config_path = Path("..") / "models" / "stage1_instrument.toml"
with open(config_path, "rb") as f:
instrument_config = tomllib.load(f)
tel_lib, coll_lib, grism_lib, cam_lib, det_lib = build_instrument_from_config(instrument_config)
from dispcraft.optics import trace_instrument
# Same sources as Section 11's "Full 2D Chain — Test"
sources = {
"on-axis" : np.array([0., 0. ]),
"alpha offset" : np.array([1e-4, 0. ]),
"delta offset" : np.array([0., 1e-4 ]),
"both offsets" : np.array([1e-4, 1e-4 ]),
}
fig, ax = plt.subplots(figsize=(6, 5))
for label, sky in sources.items():
pix = trace_instrument(sky, wl, tel_lib, coll_lib, grism_lib, cam_lib, det_lib)
ax.plot(pix[0], pix[1], label=label)
ax.set_xlabel("pixel x")
ax.set_ylabel("pixel y")
ax.legend()
ax.set_title("Spectral traces — dispcraft.optics (library-built instrument)")
plt.tight_layout()
plt.show()
Verification — Same Numbers as Section 10¶
Same two checks as Section 10, now computed from the library-built instrument.
Dispersion runs along pixel_y in this 2D model (the grism disperses along the
"y' axis" of its own, possibly tilted, frame — see Grism.forward — and with
the small tilt_deg = 2.0 from the config that axis stays close to pixel_y),
so the resolution check below reads off pixel[1] instead of pixel[0].
# Undeflected wavelength
lam_fine = np.linspace(1.0, 1.4, 4000)
dev_fine = grism_lib.deviation(lam_fine)
lambda0_lib = lam_fine[np.argmin(np.abs(dev_fine))]
print(f"undeflected wavelength ~ {lambda0_lib:.3f} µm (target: ~1.2 µm)")
# Spectral resolution (on-axis source)
pixel_onaxis = trace_instrument(np.array([0., 0.]), wl, tel_lib, coll_lib, grism_lib, cam_lib, det_lib)
dlambda_dpixel_lib = np.gradient(wl, pixel_onaxis[1]) * 1000 # µm -> nm per pixel
resolution_lib = np.mean(np.abs(dlambda_dpixel_lib))
print(f"spectral resolution ~ {resolution_lib:.3f} nm/pix (target: ~1.372 nm/pix)")
undeflected wavelength ~ 1.205 µm (target: ~1.2 µm) spectral resolution ~ 1.301 nm/pix (target: ~1.372 nm/pix)
This is the tested, promoted version of the 2D chain built above in Section
11 — same numbers, now backed by dispcraft/ (tested in tests/test_optics.py)
and models/stage1_instrument.toml, instead of inline notebook classes and
hardcoded parameters.