spacetelescope / PASTIS

Algorithm for analytical contrast predictions of coronagraphs on segmented telescopes
BSD 3-Clause "New" or "Revised" License
8 stars 3 forks source link

Investigate regression test failure due to hcipy 0.5.0 upgrade #140

Open ivalaginja opened 2 years ago

ivalaginja commented 2 years ago

See #137

ivalaginja commented 2 years ago

See also https://github.com/spacetelescope/PASTIS/pull/130#issuecomment-1227997114

ivalaginja commented 1 year ago

Trying things to figure out what makes the tests fail with hcipy>0.4.0:

I found a difference in the produced efields. Here some differences between the efield in the final focal plane as calculated by calc_psf() of the telescope class - it's the imaginary part of the efield in this plane that changed significantly (bottom right plot). The top two plots are intensities, the bottom ones real and imaginary part of the efield. This was done with the normalization from the simulator, normalizing the total energy to one photon:

Screenshot 2023-05-10 at 00 10 46

The difference is also visible in the focal plane of the FPM:

Screenshot 2023-05-10 at 00 14 38

And in the Lyot plane:

Screenshot 2023-05-10 at 00 16 26

I repeated the data above without the normalization to make sure it's not that (and also centered the colorbar on zero): Screenshot 2023-05-10 at 10 54 51 Screenshot 2023-05-10 at 11 09 38 Screenshot 2023-05-10 at 11 09 54

The simulator uses a FraunhoferPropagator to get to the FPM plane and this is used only for illustration purposes. The main propagation path of the simulator uses a LyotCoronagraph to get from the compound pupil plane (all conjugate pupils treated like they were one) to get to the Lyot plane. The final propagation to the camera plane is then done with a FaunhoferPropagator.

The below example of the FraunhoferPropagator from one of the hcipy tutorials yields exactly the same result in hcipy 0.4.0 and 0.5.2. so that alone is not it, it must be something in the exact setup of the simulator in PASTIS.

wavelength = 500e-9 # meter
pupil_diameter = 15
pupil_grid = hcipy.make_pupil_grid(256, 1.2 * pupil_diameter)
aperture_luvoir = hcipy.evaluate_supersampled(hcipy.make_luvoir_a_aperture(), pupil_grid, 8)
wf_luvoir = hcipy.Wavefront(aperture_luvoir, wavelength)

focal_length = 0.5 # meter
spatial_resolution = focal_length / pupil_diameter * wavelength
focal_grid = hcipy.make_focal_grid(8, 12, spatial_resolution=spatial_resolution)

fraunhofer = hcipy.FraunhoferPropagator(pupil_grid, focal_grid, focal_length=focal_length)
img_luvoir = fraunhofer(wf_luvoir)

I repeated the same thing with reading the aperture file from disk and the difference between the two versions is again zero.

My next best guess is to have a closer look at the grids being used in the pupil and focal plane(s), and at the apodizer and how it's being sampled.

ivalaginja commented 1 year ago

@ehpor considering the plots above (edit: or below), can you tell me if anything changed in the propagators when going to hcipy 0.5.0 that might be the reason between the differences I am seeing in my simulator?

ivalaginja commented 1 year ago

I think I managed to extract the bit of code that is causing trouble. Now I just have to figure out why.

import hcipy

pxsize = 1024
diameter = 7.258823341
wvln = 500e-9
imlamD = 1.2 * 32
fpm_rad = 3.5
fpm_px = 150
sampling = 4

pupil_grid = hcipy.make_pupil_grid(dims=pxsize, diameter=diameter)
lam_over_d = wvln / diameter

aper_fname = '/Users/ilaginja/repos/PASTIS/data/SCDA/2-Hex/masks/TelAp_LUVex_02-Hex_gy_ovsamp04_N1024.fits'
pup_read = hcipy.read_fits(aper_fname)
aperture = hcipy.Field(pup_read.ravel(), pupil_grid)

# Create a focal plane mask
samp_foc = fpm_px / (fpm_rad * 2)
focal_grid_fpm = hcipy.make_focal_grid_from_pupil_grid(pupil_grid=pupil_grid, q=samp_foc, num_airy=fpm_rad, wavelength=wvln)
fpm = 1 - hcipy.circular_aperture(2 * fpm_rad * lam_over_d)(focal_grid_fpm)

# Create a focal plane grid for the detector
focal_det = hcipy.make_focal_grid_from_pupil_grid(pupil_grid=pupil_grid, q=sampling, num_airy=imlamD, wavelength=wvln)

prop_method = hcipy.FraunhoferPropagator(pupil_grid, focal_det)

wf_active_pupil = hcipy.Wavefront(aperture, wavelength=wvln)
wf_before_fpm = prop_method(wf_active_pupil)

It looks like I just have a shift somewhere, but I am not quite sure how this would have changed from one hcipy version to another: Screenshot 2023-05-10 at 12 03 11

@ehpor I can see that hcipy 0.5.0 introduced some changes in the shift properties of the FFT, maybe that's it?

ehpor commented 1 year ago

@ivalaginja Can you check if make_focal_grid_from_pupil_grid() returns the exact same grid for both HCIPy 0.5.2dev and HCIPy 0.4.0. You should only compare images directly if they are defined on the same grid.

Also, can you check that this is the case for HCIPy 0.5.0 too, not just HCIPy 0.5.2dev? (To narrow down the possible changes.)

ivalaginja commented 1 year ago

Will do. The images are indeed not the same so I doubt the grid is, so much I know. The change made all tests explode though - will report back when I get my hands on this again. Good call to actually compare to 0.5.0.