Add FQPM coronagraph

[ivalaginja]

Fixes #164.

I have implemented the restructuring suggested by @syhaffert in #164. Notes:

• I replaced np.sign() in creating the FQPM with a subtraction of two orthogonal heaviside functions which gives a better attenuation.
• The location of the new classes was chosen by me arbitrarily and I'll be happy to change that.
• The default parameters for the VortexCoronagraph don't work as well for the FQPMCoronagraph, so I adapted them. There might still be parameter combinations that work better than what I put in here.
• The FQPMCoronagraph needs a focal-plane grid that does not contain the origin (centered between pixels instead of on a pixel) to work well. I have added this tentatively in commit bf12a81 but this breaks the VortexCoronagraph. Waiting on some advice from @syhaffert to find a solution for flexible focal-plane grids - not sure if a simple if statement inside of MultiScaleCoronagraph would suffice... for example to test a hard-coded boolean attribute inside of the subclasses? Something like self.grid_contains_origin = True/False?
• Missing tests, missing documentation.

[codecov]

Codecov Report

Merging #165 (e6feab9) into master (3249f92) will increase coverage by 0.04%.
The diff coverage is 94.81%.

@@            Coverage Diff             @@
##           master     #165      +/-   ##
==========================================
+ Coverage   81.85%   81.89%   +0.04%     
==========================================
  Files          97       99       +2     
  Lines        7328     7347      +19     
==========================================
+ Hits         5998     6017      +19     
  Misses       1330     1330              

Files	Coverage Δ
hcipy/coronagraphy/__init__.py	100.00% <100.00%> (ø)
hcipy/coronagraphy/fqpm.py	100.00% <100.00%> (ø)
hcipy/coronagraphy/multi_scale.py	95.71% <95.71%> (ø)
hcipy/coronagraphy/vortex.py	94.36% <93.98%> (-0.19%)	⬇️

📣 Codecov offers a browser extension for seamless coverage viewing on GitHub. Try it in Chrome or Firefox today!

[ivalaginja]

Using a test like already exists for the VortexCoronagraph, we can see that this FQPM doesn't reach quite the same result. Below the test, printing the resulting values instead of using assertions:

def test_fqpm_coronagraph():
	pupil_grid = make_pupil_grid(256)
	focal_grid = make_focal_grid(4, 32)
	prop = FraunhoferPropagator(pupil_grid, focal_grid)

	aperture = make_circular_aperture(1)
	aperture = evaluate_supersampled(aperture, pupil_grid, 8)

	lyot = make_circular_aperture(0.99)
	lyot = evaluate_supersampled(lyot, pupil_grid, 8) > 1 - 1e-5

	fqpm = FQPMCoronagraph(pupil_grid, q=4, scaling_factor=2, window_size=256)

	wf = Wavefront(aperture)
	wf.total_power = 1

	img_ref = prop(wf)

	wf = fqpm(wf)
	wf.electric_field *= lyot
	img = prop(wf)

        # assert img.total_power < 1e-6
        # assert img.intensity.max() / img_ref.intensity.max() < 1e-8

	print(img.total_power)
	print(img.intensity.max() / img_ref.intensity.max())

test_fqpm_coronagraph()

And the output - note that this uses make_pupil_grid() to create a grid for the phase mask (FQPM) in the focal plane:

Total power: 1.3609898015020445e-05
Normalized max intensity: 1.8636041696954595e-07

The test for the Vortex found here:

hcipy/tests/test_coronagraphy.py

Line 4 in 3a6c42f

def test_vortex_coronagraph():

returns the following results - note that this now uses make_focal_grid() to create a grid for the phase masks (Vortex of different charges) in the focal plane:

Charge 2
Total power: 1.7280129555267325e-07
Normalized max intensity: 6.274304397647414e-10

Charge 4
Total power: 4.65214713034196e-07
Normalized max intensity: 2.678918513019481e-09

Charge 6
Total power: 7.743935824913636e-07
Normalized max intensity: 6.76816821429226e-09

Charge 8
Total power: 7.590527468168938e-07
Normalized max intensity: 5.038705488143835e-09

[ivalaginja]

@syhaffert any thoughts on this?

[syhaffert]

I have been running through some tests and I think for the FQPM the high resolution is necessary along the edges of the focal plane mask. Not only in the middle. The multiscalecoronagraph zooms in at the center of the mask, exactly where the singularity of the vortex is. For the FQPM we should zoom in along the edges. I have made something that might work, but I still need to figure out what the correct parameters are. Its based on the KnifeEdgeCoronagraph.

[ivalaginja]

Any news @syhaffert ? (This is not pressing in any way, just checking in.)

[syhaffert]

Not yet. I still have to resolve a pixel offset somewhere in the code. Very frustrating last bug.

[ivalaginja]

Any luck with that bug?

[ivalaginja]

Happy new year! Any news on this?

[ivalaginja]

Motivate by last week's workshop, I just wanted to up this here again ^

[syhaffert]

@ehpor and me discussed a better implementation last week using adaptive sampling and multi-scale transforms. I am now on an observing run and I don't know if I have time to implement it. But usually during some slow observing nights I start a lot of hcipy coding, so I might have a first go at it.

[ivalaginja]

Any news on this? Pierre got curious about the FQPM in hcipy in light of having a new intern and also looking ahead to the CDI workshop in June.

[ivalaginja]

General update:

@ehpor helped me figure out how to feed a FQPM phase mask properly into the multi-scale class. I adapted the code accordingly and I think for the matter of adding a FQPM to hcipy, this is good to go. If there is ever the need to adapt the multi-scaling to the discontinuities of the FQPM in particular, this can be done separately. I don't think anybody needs this right at this moment.

I have double-checked that the VortexCoronagraph gives exactly the same result before and after this PR, the Lyot plane array and the final focal plane array subtract to zero exactly, respectively.

The FQPM propagation with the defaults from this PR look as follows (note the large range of lam/D in the focal plane for the display):

Code this was created with:

import hcipy
from matplotlib.colors import LogNorm
import matplotlib.pyplot as plt
import numpy as np


pup_size_px = 512    # [px]
pup_size_m = 2.0   # [m]
pupil_grid = hcipy.make_pupil_grid(pup_size_px, pup_size_m)

aperture = hcipy.evaluate_supersampled(hcipy.make_circular_aperture(1), pupil_grid, 32)
ls_undersizing = 0.89
lyot_mask = hcipy.make_circular_aperture(ls_undersizing)(pupil_grid)
lyot_stop = hcipy.Apodizer(lyot_mask)

samp = 8    # px per lam/D
imglamD = 200    # focal-plane images half size in units of lam/D
focal_px = samp * imglamD * 2
focal_grid_fpm = hcipy.make_uniform_grid(dims=(focal_px, focal_px), extent=imglamD * 2)
prop = hcipy.FraunhoferPropagator(pupil_grid, focal_grid_fpm, focal_length=1)

coro = hcipy.FQPMCoronagraph(pupil_grid)

wf = hcipy.Wavefront(aperture)
lyot_plane = coro(wf)
post_lyot_mask = lyot_stop(lyot_plane)
img = prop(post_lyot_mask).intensity

plt.figure(figsize=(18, 9))
plt.subplot(1, 2, 1)
hcipy.imshow_field(lyot_plane.intensity, cmap='inferno', norm=LogNorm())
plt.title('Lyot plane')
plt.colorbar()
plt.subplot(1, 2, 2)
hcipy.imshow_field(np.log10(img), cmap='inferno')
plt.title('Coro PSF')
plt.colorbar()

[ivalaginja]

Ready for review @syhaffert - just fyi @ehpor in case you had anything to add, since you talked both to Sebastiaan and me about it.

[ivalaginja]

Getting the same result like above after the rebase and reformatting to space indentations, so still ready for a review.

[ivalaginja]

@syhaffert @ehpor anybody have some time for a review? :)

[syhaffert]

Sorry for the late review. I think this is pretty much what I coded up too a couple of months ago. So I think it's a good approach. The only comment I have is that the documentation for the MultiScaleCoronagraph still explicitly mentions the vortex coronagraph.

[syhaffert]

I think everything looks good. Please address the documentation of the backward propagation. After that it is ready for merging.

[ivalaginja]

Sounds good, and all done!

[ivalaginja]

We're so close ^^ Anything else on this @syhaffert or are we good to go? :)

[syhaffert]

This is ready for merging.