Polychromatic PSF Simulation¶

What you will learn:

Why a polychromatic (broadband) PSF differs from a monochromatic one
How stellar black-body spectra are modelled in PSFCraft
How to compute a flux-weighted polychromatic PSF using build_polychromatic_star
How stellar temperature affects the PSF through chromatic broadening

Prerequisites: 04_wavefront_aberrations.ipynb

Background¶

Real astronomical sources emit light across a broad spectral range.
The PSF at wavelength $\lambda$ scales as $\lambda / D$, so different wavelengths produce PSFs of different sizes.

A polychromatic PSF is the flux-weighted sum of monochromatic PSFs:

$$\mathrm{PSF}_\mathrm{poly}(\vec{x}) = \frac{\sum_i w_i \, \mathrm{PSF}(\vec{x}, \lambda_i)}{\sum_i w_i}$$

where the weights $w_i \propto F_\nu(\lambda_i)$ are the stellar spectral flux values sampled across the filter bandpass.

PSFCraft uses the pysynphot BlackBody model to generate realistic stellar spectra.

In [1]:

Copied!





%matplotlib inline
import psfcraft
from psfcraft.utils import build_polychromatic_star
from psfcraft import constants as csts
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import psfcraft
from psfcraft.utils import build_polychromatic_star
from psfcraft import constants as csts
import numpy as np
import matplotlib.pyplot as plt

pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.

1. Inspecting Stellar Black-Body Spectra¶

PSFCraft's constants module defines standard photometric filters (Y, J, H bands) and reference stellar temperatures. Let's visualise the spectral content.

In [2]:

Copied!





def planck_law(wavelength_m, temperature_K):
    """Spectral radiance of a black body (W sr⁻¹ m⁻² m⁻¹)."""
    from astropy.constants import c, h, k_B
    c_, h_, k_ = c.value, h.value, k_B.value
    return 2 * h_ * c_**2 / (wavelength_m**5 * (np.exp(h_ * c_ / (wavelength_m * k_ * temperature_K)) - 1))

# Wavelength axis
wl_m = np.linspace(0.3e-6, 2.5e-6, 500)

# Stellar types and their approximate temperatures
stellar_types = {
    "O  (40 000 K)": 40_000,
    "A  ( 8 200 K)": 8_200,
    "G  ( 5 800 K)": 5_800,
    "M  ( 3 300 K)": 3_300,
}
colors = ["royalblue", "skyblue", "gold", "firebrick"]

fig, ax = plt.subplots(figsize=(10, 4))
for (label, T), color in zip(stellar_types.items(), colors):
    flux = planck_law(wl_m, T)
    ax.plot(wl_m * 1e6, flux / flux.max(), label=label, color=color, linewidth=2)

# Highlight photometric bands
band_info = [
    ("Y", csts.FILTER_Y_CUTS, "khaki"),
    ("J", csts.FILTER_J_CUTS, "orange"),
    ("H", csts.FILTER_H_CUTS, "salmon"),
]
for name, (lo, hi), color in band_info:
    ax.axvspan(lo * 1e-3, hi * 1e-3, alpha=0.25, color=color, label=f"{name} band")
    ax.text((lo + hi) / 2 * 1e-3, 0.05, name, ha="center", fontsize=11)

ax.set_xlabel("Wavelength (µm)")
ax.set_ylabel("Normalised flux")
ax.set_title("Black-body spectra and NIR filter bands")
ax.legend(loc="upper right", fontsize=9)
ax.set_xlim(0.3, 2.5)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
def planck_law(wavelength_m, temperature_K):
    """Spectral radiance of a black body (W sr⁻¹ m⁻² m⁻¹)."""
    from astropy.constants import c, h, k_B
    c_, h_, k_ = c.value, h.value, k_B.value
    return 2 * h_ * c_**2 / (wavelength_m**5 * (np.exp(h_ * c_ / (wavelength_m * k_ * temperature_K)) - 1))

# Wavelength axis
wl_m = np.linspace(0.3e-6, 2.5e-6, 500)

# Stellar types and their approximate temperatures
stellar_types = {
    "O  (40 000 K)": 40_000,
    "A  ( 8 200 K)": 8_200,
    "G  ( 5 800 K)": 5_800,
    "M  ( 3 300 K)": 3_300,
}
colors = ["royalblue", "skyblue", "gold", "firebrick"]

fig, ax = plt.subplots(figsize=(10, 4))
for (label, T), color in zip(stellar_types.items(), colors):
    flux = planck_law(wl_m, T)
    ax.plot(wl_m * 1e6, flux / flux.max(), label=label, color=color, linewidth=2)

# Highlight photometric bands
band_info = [
    ("Y", csts.FILTER_Y_CUTS, "khaki"),
    ("J", csts.FILTER_J_CUTS, "orange"),
    ("H", csts.FILTER_H_CUTS, "salmon"),
]
for name, (lo, hi), color in band_info:
    ax.axvspan(lo * 1e-3, hi * 1e-3, alpha=0.25, color=color, label=f"{name} band")
    ax.text((lo + hi) / 2 * 1e-3, 0.05, name, ha="center", fontsize=11)

ax.set_xlabel("Wavelength (µm)")
ax.set_ylabel("Normalised flux")
ax.set_title("Black-body spectra and NIR filter bands")
ax.legend(loc="upper right", fontsize=9)
ax.set_xlim(0.3, 2.5)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

No description has been provided for this image

2. Building a Polychromatic Source Dictionary¶

build_polychromatic_star(temperature, filter_name, sampling) returns a dict with keys:

wavelengths — array of wavelengths in metres
weights — corresponding spectral flux values (Planck black-body at temperature K)

This dictionary is passed directly to calc_psf(source=...).

In [3]:

Copied!





# Reload psfcraft.utils so the kernel picks up the updated build_polychromatic_star
import importlib, psfcraft.utils
from psfcraft.utils import build_polychromatic_star

# Build a G-star source (T = 5800 K) in the J band, sampled at 5 wavelengths
source = build_polychromatic_star(temperature=5800, filter_name="J", sampling=5)

print("Source wavelengths (µm):", np.round(source["wavelengths"] * 1e6, 3))
print("Source weights (normalised):", np.round(source["weights"] / source["weights"].max(), 3))

# Plot the sampled spectrum
fig, ax = plt.subplots(figsize=(7, 3))
ax.stem(source["wavelengths"] * 1e6, source["weights"] / source["weights"].max(),
        linefmt="C0-", markerfmt="C0o", basefmt="k-")
lo, hi = csts.FILTER_J_CUTS
ax.axvspan(lo * 1e-3, hi * 1e-3, alpha=0.15, color="orange", label="J band")
ax.set_xlabel("Wavelength (µm)")
ax.set_ylabel("Normalised flux")
ax.set_title("G-star (5800 K) sampled in the J band — 5 wavelengths")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Reload psfcraft.utils so the kernel picks up the updated build_polychromatic_star
import importlib, psfcraft.utils
from psfcraft.utils import build_polychromatic_star

# Build a G-star source (T = 5800 K) in the J band, sampled at 5 wavelengths
source = build_polychromatic_star(temperature=5800, filter_name="J", sampling=5)

print("Source wavelengths (µm):", np.round(source["wavelengths"] * 1e6, 3))
print("Source weights (normalised):", np.round(source["weights"] / source["weights"].max(), 3))

# Plot the sampled spectrum
fig, ax = plt.subplots(figsize=(7, 3))
ax.stem(source["wavelengths"] * 1e6, source["weights"] / source["weights"].max(),
        linefmt="C0-", markerfmt="C0o", basefmt="k-")
lo, hi = csts.FILTER_J_CUTS
ax.axvspan(lo * 1e-3, hi * 1e-3, alpha=0.15, color="orange", label="J band")
ax.set_xlabel("Wavelength (µm)")
ax.set_ylabel("Normalised flux")
ax.set_title("G-star (5800 K) sampled in the J band — 5 wavelengths")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Source wavelengths (µm): [1.168 1.268 1.368 1.467 1.567]
Source weights (normalised): [1.    0.804 0.652 0.532 0.438]

3. Computing a Polychromatic PSF¶

Pass the source dictionary to calc_psf(source=...) instead of monochromatic=. POPPY will compute one PSF per wavelength sample and combine them with the given weights.

In [4]:

Copied!





tel = psfcraft.NewtonianTelescope(version="1_3")
tel.pixelscale = 0.0298   # Euclid pixel scale

# Monochromatic reference (J-band centre)
wl_centre = np.mean(source["wavelengths"])
psf_mono  = tel.calc_psf(monochromatic=wl_centre, fov_pixels=64)

# Polychromatic
psf_poly  = tel.calc_psf(source=source, fov_pixels=64)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
psfcraft.display_psf(psf_mono, ax=axes[0],
                     title=f"Monochromatic  λ = {wl_centre*1e6:.2f} µm",
                     colorbar=False)
psfcraft.display_psf(psf_poly, ax=axes[1],
                     title="Polychromatic (J band, G star)",
                     colorbar=False)
plt.suptitle("Monochromatic vs Polychromatic PSF", fontsize=13)
plt.tight_layout()
plt.show()
tel = psfcraft.NewtonianTelescope(version="1_3")
tel.pixelscale = 0.0298   # Euclid pixel scale

# Monochromatic reference (J-band centre)
wl_centre = np.mean(source["wavelengths"])
psf_mono  = tel.calc_psf(monochromatic=wl_centre, fov_pixels=64)

# Polychromatic
psf_poly  = tel.calc_psf(source=source, fov_pixels=64)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
psfcraft.display_psf(psf_mono, ax=axes[0],
                     title=f"Monochromatic  λ = {wl_centre*1e6:.2f} µm",
                     colorbar=False)
psfcraft.display_psf(psf_poly, ax=axes[1],
                     title="Polychromatic (J band, G star)",
                     colorbar=False)
plt.suptitle("Monochromatic vs Polychromatic PSF", fontsize=13)
plt.tight_layout()
plt.show()

pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.

4. Effect of Stellar Temperature on the PSF¶

Hotter stars emit more flux at shorter wavelengths → their PSF is slightly sharper (shorter-λ component dominates the weighted sum).
Cooler stars are redder → the PSF is broader.

In [5]:

Copied!





temperatures = {
    "O  (40 000 K)": 40_000,
    "A  ( 8 200 K)": 8_200,
    "G  ( 5 800 K)": 5_800,
    "M  ( 3 300 K)": 3_300,
}
filter_name = "J"

tel = psfcraft.NewtonianTelescope(version="1_3")
tel.pixelscale = 0.0298

fig, axes = plt.subplots(1, len(temperatures), figsize=(14, 4))
for ax, (label, T) in zip(axes, temperatures.items()):
    src = build_polychromatic_star(temperature=T, filter_name=filter_name, sampling=7)
    psf = tel.calc_psf(source=src, fov_pixels=64)
    psfcraft.display_psf(psf, ax=ax, title=label, colorbar=False)

plt.suptitle(f"Polychromatic PSF — {filter_name}-band, varying stellar temperature", fontsize=12)
plt.tight_layout()
plt.show()
temperatures = {
    "O  (40 000 K)": 40_000,
    "A  ( 8 200 K)": 8_200,
    "G  ( 5 800 K)": 5_800,
    "M  ( 3 300 K)": 3_300,
}
filter_name = "J"

tel = psfcraft.NewtonianTelescope(version="1_3")
tel.pixelscale = 0.0298

fig, axes = plt.subplots(1, len(temperatures), figsize=(14, 4))
for ax, (label, T) in zip(axes, temperatures.items()):
    src = build_polychromatic_star(temperature=T, filter_name=filter_name, sampling=7)
    psf = tel.calc_psf(source=src, fov_pixels=64)
    psfcraft.display_psf(psf, ax=ax, title=label, colorbar=False)

plt.suptitle(f"Polychromatic PSF — {filter_name}-band, varying stellar temperature", fontsize=12)
plt.tight_layout()
plt.show()

pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.
pysynphot is unavailable (`np.alltrue` was removed in the NumPy 2.0 release. Use `np.all` instead.). Polychromatic sources will use the built-in Planck-law implementation.

Key Takeaways¶

A polychromatic PSF is the flux-weighted average of monochromatic PSFs across the bandpass.
Use build_polychromatic_star(T, filter_name, sampling) to build the source dictionary (requires pysynphot).
Pass source=... to calc_psf() instead of monochromatic= to trigger polychromatic mode.
Hotter (bluer) stars produce slightly sharper PSFs in a given band; cooler (redder) stars produce broader ones.
More wavelength samples → better accuracy but longer computation time; 5–10 samples is usually sufficient.

Next: 07_psf_generation_pipeline.ipynb — build a full batch PSF dataset for ML or analysis pipelines.