ehpor / hcipy

A framework for performing optical propagation simulations, meant for high contrast imaging, in Python.
https://hcipy.org
MIT License
91 stars 30 forks source link

Deterministic atmospheric phase screens #173 #189

Closed Oswald522 closed 1 year ago

Oswald522 commented 1 year ago

Thanks for your outstanding works. I have read the #169 and #173 When using the layer.reset() , I cannot get the initial state(Expectation: the InfiniteAtmosphericLayer hould return the same phase screen every time.). or maybe print(layer.phase_for(1)) should be same?. the code follows:

from hcipy import *
import numpy as np
wavelength = 0.5e-6
D_tel = 1
fried_parameter = 0.1
outer_scale = 10
propagate_phase_screen = False
velocity = 10.0  # meters/sec
num_modes = 1000
pupil_grid = make_pupil_grid(64, D_tel)
aperture = make_circular_aperture(D_tel)(pupil_grid)
Cn_squared = Cn_squared_from_fried_parameter(fried_parameter, wavelength)
## add the seed
layer = InfiniteAtmosphericLayer(pupil_grid, Cn_squared, outer_scale, [velocity / np.sqrt(2), velocity / np.sqrt(2)],
                                 use_interpolation=False,seed=1)
layer.reset()
print(layer.phase_for(1))
layer.reset()
print(layer.phase_for(1))
layer.reset(make_independent_realization=True)
print(layer.phase_for(1))

and the outputs vary every time

[-2.57236972e-06 -2.70448871e-06 -2.75221671e-06 ...  4.62491109e-06
  4.63627301e-06  4.81470425e-06]
[ 2.06733563e-06  2.08069233e-06  2.17932777e-06 ... -1.82349645e-07
 -2.99504848e-07 -6.31952191e-08]
[ 6.92540580e-07  4.12439742e-07  2.36030268e-07 ... -4.45835357e-07
 -2.77815848e-07 -3.59797936e-07]

Is there anything wrong?

ehpor commented 1 year ago

Confirmed. Thanks for the code example. I'll look at this in the afternoon.

Oswald522 commented 1 year ago

I may have found something useful. The results differ due to the reset function in the hcipy/atmosphere/finite_atmospheric_layer.py file, specifically the self._noise = self.noise_factory.make_random() function. Or it could be changed to self._noise = self.noise_factory.make_random(self.rng). I am not very familiar with hcipy, so I hope I am not giving you inaccurate information.