Skip to content

Add PSF model based on Zernike coefficients - #3056

Open
mexanick wants to merge 14 commits into
mainfrom
zernike-psf-model
Open

Add PSF model based on Zernike coefficients#3056
mexanick wants to merge 14 commits into
mainfrom
zernike-psf-model

Conversation

@mexanick

Copy link
Copy Markdown
Contributor

Summary

Adds ZernikePSFModel, to ctapipe.instrument.optics alongside the existing ComaPSFModel.

Unlike ComaPSFModel's parametric radial/polar scaling approach, this model reconstructs the optical wavefront directly from Zernike coefficients (using Noll indexing) and propagates it through Fraunhofer diffraction (FFT of the aberrated pupil).

Implementation details

  • Builds a circular pupil grid (pupil_size, pupil_diameter_fraction) with a soft-edged aperture (pupil_edge_softness) to avoid FFT ringing.
  • Reconstructs the wavefront from per-telescope Zernike coefficients z2z11 (tilt, defocus, astigmatism, coma, trefoil, spherical), using the zernike package's RZern for polynomial evaluation.
  • Adds off-axis aberration growth: coma grows linearly with field angle θ (z7_theta, z8_theta), astigmatism grows quadratically with θ² (z5_theta2, z6_theta2).
  • Computes the focal-plane intensity via FFT of the complex pupil function, averaged over wavelength_samples wavelengths between wavelength_min and wavelength_max, weighted by a cherenkov_spectrum_index power law (dN/dλ ∝ λ^-index).
  • Applies Gaussian smoothing (focal_plane_smoothing_sigma_pix) to the resulting intensity to approximate pixel/light-guide response.
  • Exposes the standard pdf(tel_id, lon, lat, lon0, lat0) interface shared with ComaPSFModel, interpolating the precomputed PSF grid via scipy.ndimage.map_coordinates and normalizing to unit integral.
image

@mexanick
mexanick requested review from Voutsi, kosack and maxnoe July 17, 2026 15:57
@mexanick
mexanick marked this pull request as ready for review July 17, 2026 15:58

@kosack kosack left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good, but a couple of small requests:

  • Can you add a citation to the docs to where the "Zernike" coefficients are defined?
  • Just a question: you make all the parameters generic names like "z2-z6", but they could have nicer names, e.g. z2 could be tilt_x, etc. Maybe that is because of how they are defined in the reference paper? If so, then adding the citation is fine, but otherwise, why not use human-readable names?

@mexanick

Copy link
Copy Markdown
Contributor Author
  • Can you add a citation to the docs to where the "Zernike" coefficients are defined?

I've added a reference to Wikipedia. The original Zernike article is in German, and Noll's article is behind a paywall, but Wikipedia provides all the necessary information.

  • Just a question: you make all the parameters generic names like "z2-z6", but they could have nicer names, e.g. z2 could be tilt_x, etc. Maybe that is because of how they are defined in the reference paper? If so, then adding the citation is fine, but otherwise, why not use human-readable names?

I would like to keep them as is (Zi), as this is the industry standard. For example, Ansys Zemax uses this notation, which makes it easy to copy or cross-check coefficient values directly from ray-tracing software. The docstrings for each parameter already contain the human-readable descriptions, so I think it's a good compromise between readability and standardization.

@mexanick
mexanick requested a review from kosack July 27, 2026 08:24
kosack
kosack previously approved these changes Jul 27, 2026
Comment thread src/ctapipe/instrument/optics.py Outdated
self.z5.tel[tel_id].to_value(u.m) + z5_theta2_m_per_deg2 * theta2,
self.z6.tel[tel_id].to_value(u.m) + z6_theta2_m_per_deg2 * theta2,
self.z7.tel[tel_id].to_value(u.m) + coma_x,
self.z8.tel[tel_id].to_value(u.m) + coma_y,

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

According to wikipedia link, the Noll coefficient z7 is the vertical and z8 the horizontal. I think then you need to swap the comas, I mean add coma_x to z8 and coma_y to z7

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was an artifact of some experimenting with the internal to the model coordinate transforms, I've simplified it.

Comment thread src/ctapipe/instrument/optics.py Outdated
phase_safe[mask] = phase[mask]
pupil = aperture * np.exp(1j * phase_safe)
field = fft2(pupil)
intensity += weight * fftshift(np.abs(field) ** 2)

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The angular coordinates depend on the wavelength. I think here we ignore this dependence, it should be accounted for before the final summation.

@Voutsi

Voutsi commented Jul 31, 2026

Copy link
Copy Markdown

Hi @mexanick , thanks a lot for this. I asked for one more thing to be clarified. Other than that looks good to me.

@mexanick

Copy link
Copy Markdown
Contributor Author

Hi @mexanick , thanks a lot for this. I asked for one more thing to be clarified. Other than that looks good to me.

Thanks, it is fixed now. I also found and fixed an error in the astigmatism application. An updated PSF map:
image

@mexanick
mexanick requested a review from Voutsi July 31, 2026 17:14
Comment thread src/ctapipe/instrument/optics.py Outdated
Comment thread src/ctapipe/instrument/optics.py Outdated
Comment thread src/ctapipe/instrument/optics.py
Comment thread pyproject.toml Outdated

@maxnoe maxnoe left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, this looks very good.

Aside from the comment by @Voutsi, I'd just request that we add the new package only as an optional dependency.

@mexanick
mexanick force-pushed the zernike-psf-model branch from 277728b to 663d9ed Compare August 3, 2026 16:03
@mexanick
mexanick requested a review from maxnoe August 3, 2026 16:06
Comment thread docs/conf.py
@ctao-sonarqube

ctao-sonarqube Bot commented Aug 4, 2026

Copy link
Copy Markdown

Comment thread src/ctapipe/instrument/optics.py
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants