Introduction

Modeling a Spectroscopic Instrument

1. Dispersion model

The Central Question

Given:
· a source at sky position (α, δ)
· emitting light at wavelength λ

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)

2. Geometrical Optics

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]
· θ — angle to the optical axis [rad]

As the ray passes through each element, (x, θ) 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.
f x θ = 0 θ' lens (f) focal plane x'

A ray enters at height x parallel to the axis (θ = 0); the lens deflects it by −x/f toward the focal point.

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 (θ = 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, θ).

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)
α → telescope → collimator → grism → camera → xpix

The grism adds δ(λ) along this axis only.
\[ x_\text{pix} = f(\alpha, \lambda) \]
Cross-dispersion direction (y)
δ → telescope → collimator → camera → ypix

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 α 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.

3. First Implementation

Step 1 — Naive Implementation

Python concept: variable, expression

Start simple: one ray, one lens. Write the equation directly as Python expressions. 

# 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
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)

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
A function has a name, receives inputs (parameters), and produces outputs.
Write it once, fix it once.

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. 
import numpy as np

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. ]
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.

Step 3 — Test the Function

Exercise: run this in your notebook

# 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

4. Collimator and Camera

Collimator — Position → Angle

Role: placed at distance fcoll 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}} \] (source above axis → beam going downward: minus sign)
This is the thin lens formula applied to a point source at the focal plane: θin = 0, so θout = 0 − x/f = −x/f.
focal plane f x θ = −x/f collimator (f)

A point source at position x in the focal plane → all rays exit as a parallel beam at angle θ = −x/f.

Camera — Angle → Position

Role: re-focuses a parallel beam onto the detector.

A collimated beam at angle θ converges to a point at: \[ x = -f_\text{cam}\,\theta \] (positive angle → negative position: image is flipped)
Collimator + Camera together: \[ x_\text{det} = \frac{f_\text{cam}}{f_\text{coll}}\,x_\text{focal} \] The image is magnified by γ = fcam/fcoll.
camera (f) f detector θ x x = −f·θ

Three parallel rays at angle θ enter the camera and converge to a single point at x = −f·θ on the detector.

The Grism acts between them.
With the grism: θ → θ + δ(λ) → camera → dispersed image.
The extra δ(λ) shifts the image by −fcam·δ(λ), a different amount for each wavelength.

Collimator & Camera — Named Functions

Python concept: one function, one clear responsibility

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
Each function takes one physical quantity as input and returns one.
Named, single-purpose functions make the chain explicit.

Collimator + Camera — Test

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 γ = f_cam / f_coll = 0.3/0.5 = 0.6
# x_det = x_focal * 0.6 = 0.60 mm  ✓

5. Telescope and Classes

Telescope

Role: convert a sky angle into a position in the focal plane.
The telescope is pointed at a centre α0 (the optical axis). A source at sky angle α produces an angular offset α − α0: \[ x_\text{focal} = -f_\text{tel}\,(\alpha - \alpha_0) \] Under the paraxial approximation, α0 is simply subtracted before applying the focal length.
This is the same formula as the Camera: \( x = -f \cdot \theta \) with \(\theta = \alpha - \alpha_0\).

Euclid: ftel = 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

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
__init__ runs when the object is created — stores parameters.
self is the object itself — it carries its own data.
trace() is a method: a function that belongs to the class.

Telescope — Inheriting from Camera

Python concept: inheritance + extending with a new parameter

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)
super().__init__(f) — calls Camera.__init__
Inheritance lets us extend a class without rewriting it.

Collimator as a Class

Every element has the same interface: .trace(input)

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
Every element follows the same interface: .forward(input).
Reading the instrument chain becomes as natural as reading the pipeline diagram.

6. Refractive Index

Light in a Medium

Refractive index n
Light slows down in a transparent medium: \[ v = \frac{c}{n} \] At each interface, the ray bends (Snell's law).
Dispersion: n = n(λ)
n depends on wavelength.
Blue (short λ) slows more than red (long λ).
→ different wavelengths exit at different angles.
This is the physical origin of spectral dispersion.
glass n(λ) white light λ short (blue) λ long (red) n large → bends more

Short wavelengths (large n) refract more than long wavelengths — the physical origin of dispersion.

Cauchy model (1836): \[ n(\lambda) = n_0 + \frac{k}{\lambda^2} \] λ in µm; n0 and k depend on the glass type.

Material — Implementation

Python concept: class with a method that computes a value

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 λ)

A Quick Word on matplotlib

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.
Core functions used here
  • plt.figure(figsize=(w,h)) — new figure
  • plt.plot(x, y) — line plot
  • plt.scatter(x, y, c=...) — scatter plot
  • plt.xlabel / ylabel / title — labels
  • plt.legend() — legend
  • plt.show() — display
Install
pip install matplotlib

Documentation
matplotlib.org/stable

Quick intro
pyplot tutorial

Material — Visualization

import matplotlib.pyplot as plt

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()
n decreases with λ (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.

7. The Prism

Prism — Physics

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 λ, different wavelengths are deflected by different amounts — the dispersive action of the prism.
prism n(λ) A white 0.9 µm 2.0 µm

A thin prism of apex angle A deviates each wavelength by δ(λ) = A·(n(λ)−1).

Prism — Implementation

Python concept: Prism has-a Material (composition)

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°
print(np.degrees(prism.deviation(1.0)))   # deviation at 1 µm [deg]
self.material = material — the Prism stores a reference to a Material object.
This is composition: one object contains another.

8. The Diffraction Grating

Grating — Physics

Grating equation: \[ d\,(\sin\theta_m - \sin\theta_i) = m\,\lambda \] · d = groove spacing [µm]
· m = diffraction order
· θi = incident angle, θm = diffracted angle
Paraxial, normal incidence (θi ≈ 0): \[ \delta_\text{grating}(\lambda) \approx \frac{m\,\lambda}{d} \] Unlike the prism, grating dispersion is linear in λ.
d incident (θᵢ = 0) m = 0 m=+1 λ=0.9µm m=+1 λ=2.0µm

Transmission grating: light passes through and is diffracted at angle δ = mλ/d below the surface. Long wavelengths diffract more.

Grating — Implementation

Exercise: implement this class

class Grating:
    """
    Diffraction grating: angular deviation δ(λ) = m * λ / d.

    Parameters
    ----------
    m : int   — diffraction order
    d : float — groove spacing [µm]
    """

    def __init__(self, m, d):
        self.m = m
        self.d = d

    def deviation(self, wavelength):
        """Angular deviation [rad]. wavelength in µm."""
        return self.m * wavelength / self.d


# Test:
grating = Grating(m=1, d=13.7)
d1 = np.degrees(grating.deviation(1.0))
d2 = np.degrees(grating.deviation(2.0))
print(f"δ(1.0 µm) = {d1:.3f}°,  δ(2.0 µm) = {d2:.3f}°")
# Grating dispersion is linear: δ(2µm) = 2 × δ(1µm)

9. The Grism

What is a Grism?

A grism = a grating bonded to a prism.

The prism is chosen so that at the central wavelength λ0, prism and grating deviations cancel: \[ \delta_\text{prism}(\lambda_0) + \delta_\text{grating}(\lambda_0) = 0 \] → the beam passes straight through at λ0.
Other λ are deflected proportionally to (λ − λ0).
Total deviation: \[ \delta_\text{grism}(\lambda) = \delta_\text{prism}(\lambda) + \delta_\text{grating}(\lambda) \]
prism grating λ₁, λ₂, λ₃ λ₀ (δ=0) λ < λ₀ λ > λ₀ A

Grating etched on the hypotenuse of the prism. At λ₀ prism and grating deviations cancel — the beam passes straight through.

Grism — Implementation

Python concept: Grism has-a Prism and has-a Grating

class Grism:
    """
    Grism: prism + grating.
    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 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: x = −f·θ.
Same code, different name.

→ inherit, write pass.
Composition — "has-a"

class Grism(prism, grating)

A Grism HAS a Prism and HAS a Grating.
Its behavior combines two components.

→ store both, sum their deviations.
Rule: use inheritance when the new class specializes an existing one. Use composition when the new class is built from existing components.

Grism — Dispersion Curve

glass   = Material(n0=1.44, k=0.004)
prism   = Prism(material=glass, A=0.05)
grating = Grating(m=1, d=13.7)
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()
# The wavelength where δ = 0 is the blaze wavelength λ₀.

10. Full Instrument

Detector — Position to Pixel

The last element: converts physical units to pixels

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

The Complete Chain

Every element now has the same interface: .trace(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 — Exercise

# Build the instrument
tel   = Telescope(f=1.2)
coll  = Collimator(f=0.5)
glass = Material(n0=1.44, k=0.004)
grism = Grism(Prism(glass, A=0.05), Grating(m=1, d=13.7))
cam   = Camera(f=0.3)
det   = Detector(pixel_size=18e-6, x0=0.0)

# Source at field center (alpha = 0), vary wavelength
alpha      = 0.0
wavelength = np.linspace(0.9, 2.0, 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("Dispersion curve — α = 0 (field center)")
plt.tight_layout(); plt.show()

What the Model Gives Us

We have a 1D forward model:
given sky angle α and λ → predict detector pixel x.
\[ x_\text{pix} = f\!\left(\alpha,\, \lambda\right) \] The same chain applied to δ gives y independently.
Section 11 shows how to combine them into a single 2D model.
Next (Intro_model):
Compare with real NISP calibration data. Fit the parameters (ftel, A, d, …) to minimize residuals.
This simplified model neglects geometric distortions and lateral color (present in spectrogrism as higher-order terms).
The residuals will be corrected later by machine learning.

11. Extension to 2D

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 α and δ:
\[ 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.

Naïve Solution — Explicit x and y

Since everything couples, every step takes both coordinates as input

def trace_instrument_2d(alpha, delta, wavelength,
                         tel, coll, grism, cam, det):
    # telescope: sky angles → focal plane
    x_foc, y_foc     = tel.forward_x(alpha, delta), \
                        tel.forward_y(alpha, delta)

    # collimator: focal position → beam direction (both axes mix)
    theta_x           = coll.forward_x(x_foc, y_foc)
    theta_y           = coll.forward_y(x_foc, y_foc)

    # grism: adds dispersion (couples wavelength to both directions)
    theta_x, theta_y  = grism.forward_x(theta_x, theta_y, wavelength), \
                        grism.forward_y(theta_x, theta_y, wavelength)

    # camera: beam direction → detector position
    x_cam, y_cam      = cam.forward_x(theta_x, theta_y), \
                        cam.forward_y(theta_x, theta_y)

    # detector: physical position → pixel
    pixel_x           = det.forward_x(x_cam, y_cam)
    pixel_y           = det.forward_y(x_cam, y_cam)

    return pixel_x, pixel_y
Every element needs two methods, each taking two inputs. Every new element, every bug fix: done twice, in sync.

A Cleaner Design — 2D Vectors

Each quantity becomes a NumPy array of shape (2,) or (2, N)

import numpy as np

def trace_instrument_2d(sky, wavelength, tel, coll, grism, cam, det):
    """
    sky        : array [alpha, delta]  [rad]
    wavelength : float or array        [µm]
    Returns
    -------
    pixel : array [pixel_x, pixel_y]
    """
    pos_foc   = tel.forward(sky)                # (2,) → (2,)
    angle_col = coll.forward(pos_foc)           # (2,) → (2,)
    angle_gr  = grism.forward(angle_col, wavelength)   # (2,) → (2, N)
    pos_cam   = cam.forward(angle_gr)           # (2, N) → (2, N)
    return det.forward(pos_cam)                 # (2, N) → (2, N)


# Each element operates on a real 2D vector:
class Telescope:
    def __init__(self, f): self.f = f
    def forward(self, sky):
        return -self.f * sky      # broadcasts: (2,) or (2, N)

class Collimator:
    def __init__(self, f): self.f = f
    def forward(self, pos):
        return -pos / self.f
One method per element, one call per step. The chain reads exactly like the pipeline diagram. NumPy broadcasts over the wavelength axis automatically.

Grism in 2D — Rotation + Dispersion

The prism tilt is a 2D rotation; dispersion adds to one axis

def rotation_matrix(theta):
    c, s = np.cos(theta), np.sin(theta)
    return np.array([[c, -s],
                     [s,  c]])

class Grism:
    def __init__(self, prism, grating, tilt=0):
        self.prism   = prism
        self.grating = grating
        self.R       = rotation_matrix(tilt)   # prism tilt couples x,y

    def deviation(self, wavelength):
        return self.prism.deviation(wavelength) \
             + self.grating.deviation(wavelength)  # scalar or (N,)

    def forward(self, angle, wavelength):
        # rotate beam to align with prism apex
        angle = self.R @ angle                         # (2,N)
        # add dispersion along y axis (grating grooves along x)
        angle[1] += self.deviation(wavelength)         # broadcast over N
        # rotate back
        return self.R.T @ angle                        # (2,N)

Full 2D Chain — Test

# Instrument parameters (same as Section 10)
tel   = Telescope(f=1.2)
coll  = Collimator(f=0.5)
grism = Grism(prism, grating, tilt=np.radians(2))
cam   = Camera(f=0.3)
det   = Detector(pixel_size=18e-6, offset=np.array([0., 0.]))

wl = np.linspace(0.9, 2.0, 200)   # [µm]

# Grid of 4 sources: (alpha, delta) in rad
sources = {
    "on-axis"      : np.array([0.,    0.   ]),
    "α offset"     : np.array([1e-4,  0.   ]),
    "δ offset"     : np.array([0.,    1e-4 ]),
    "both offsets" : np.array([1e-4,  1e-4 ]),
}

fig, ax = plt.subplots()
for label, sky in sources.items():
    pix = trace_instrument_2d(sky, wl, tel, coll, grism, cam, det)
    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.show()
The δ-offset source traces a tilted line — not a horizontal one. This is the signature of the prism tilt coupling the two axes.

Pointing and Field Rotation

Sky position [α, δ] → focal-plane position [x, y]

The pointing is the sky position [α0, δ0] of the field centre).
Additionally the focal plane may be rotated relative to the equatorial axes by a roll angle ψ. 
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=np.zeros(2), angle=0.0):
        super().__init__(f)
        self.sky0  = np.asarray(sky0, dtype=float)
        self.angle = angle

    def forward(self, sky):
        """
        sky : (2,) array [alpha, delta]  [rad]
        Returns : (2,) focal-plane position [m]
        """
        offset = sky - self.sky0               # angular offset from pointing
        c, s = np.cos(self.angle), np.sin(self.angle)
        R    = np.array([[c, -s],
                         [s,  c]])             # roll rotation
        return super().forward(R @ offset)     # rotate, then focus


# No rotation — 1D result recovered on each axis:
tel = Telescope2D(f=24.5, sky0=np.array([0., 0.]))
print(tel.forward(np.array([1e-4, 2e-4])))   # → [-2.45e-3, -4.90e-3] m

# 45° roll — the α-offset now contributes to both x and y:
tel45 = Telescope2D(f=24.5, sky0=np.zeros(2), angle=np.radians(45))
print(tel45.forward(np.array([1e-4, 0.])))

Exercise — Build the 2D Chain

Goal: extend your notebook from Section 10 to produce a 2D detector map with multiple sources at different (α, δ) positions.
  1. Update each class (Telescope, Collimator, Camera, Detector) so that forward accepts a (2,) array and returns the same shape.
  2. Implement Grism.forward(angle, wavelength) using a rotation matrix and a per-axis dispersion offset.
  3. Define a 3×3 grid of sources: α ∈ {−1, 0, +1}·10⁻⁴ rad × δ ∈ {−1, 0, +1}·10⁻⁴ rad.
  4. Plot all 9 spectral traces. How does the prism tilt angle affect the pattern?

12. From Notebook to Library

The Problem with Notebooks Alone

The code lives inside the notebook.

In the next notebook you will need:
Telescope, Collimator, Camera, Grism, Detector, trace_instrument

Copy-pasting is:
· error-prone (which version is correct?)
· hard to maintain (a fix must be applied in every notebook)
· not reproducible
Solution: extract the code into a reusable Python module.

Step 1 — A Python Module

Concept: a .py file is an importable module

# optics.py
import numpy as np

class Material: ...
class Camera: ...
class Collimator: ...
class Telescope(Camera): ...
class Prism: ...
class Grating: ...
class Grism: ...
class Detector: ...

def trace_instrument(...): ...
# In any notebook:
from optics import Telescope, Collimator, Camera
from optics import Material, Prism, Grating, Grism, Detector
from optics import trace_instrument

tel = Telescope(f=1.2)   # unchanged

Step 2 — A Python Package

Concept: a folder with __init__.py is a package

dispcraft/
├── optics/
│   ├── __init__.py          ← marks the folder as a package
│   ├── materials.py         ← Material
│   ├── elements.py          ← Camera, Collimator, Telescope,
│   │                           Prism, Grating, Grism, Detector
│   └── instrument.py        ← trace_instrument
├── 01_intro.ipynb
├── 02_model.ipynb
└── README.md
from optics.elements import Telescope, Camera, Collimator
from optics.elements import Prism, Grating, Grism, Detector
from optics.instrument import trace_instrument

Step 3 — Docstrings and README

class Material:
    """
    Optical glass: refractive index following the Cauchy model.

    Parameters
    ----------
    n0 : float — base refractive index
    k  : float — dispersion coefficient [µm²]

    Examples
    --------
    >>> glass = Material(n0=1.44, k=0.004)
    >>> round(glass.refractive_index(1.5), 4)
    1.4418
    """
Docstrings appear with help(Material) or when hovering in an IDE.
A README.md at the root explains what the package does and how to install it.

Step 4 — Version Control with Git

# Initialize a repository
git init

# Check what has changed
git status

# Stage specific files
git add optics/elements.py

# Save a snapshot with a message
git commit -m "add Grism class with composition (Prism + Grating)"

# See the full history
git log --oneline
Commit after each logical change.
Messages should describe what was added and why.

Summary — The Software Path

  1. Notebook — interactive exploration and testing
  2. Module (.py) — extract reusable code
  3. Package (folder + __init__.py) — organize by topic
  4. Docstrings — document every class and function
  5. README — explain the package to others
  6. Git — track the full history of changes
Starting from the beginning avoids painful refactoring later and makes results reproducible.