esbo_ds/ESBO-ETC
esbo_ds/ESBO-ETC
11
0
Fork 1
Code Issues 1 Pull Requests Releases 1 Wiki Activity

Bugfix: scale integral of convolution

Browse Source
This commit is contained in:
Lukas Klass 2020-05-14 15:02:44 +02:00
parent 0d8b3949a0
commit 10038014f8
2 changed files with 11 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
esbo_etc/classes/psf/Airy.py
View File

@ -37,6 +37,7 @@ class Airy(IPSF):
self.__d_aperture = d_aperture self.__d_aperture = d_aperture
self.__osf = osf self.__osf = osf
self.__pixel_size = pixel_size self.__pixel_size = pixel_size
self.__psf_jitter = None
def calcReducedObservationAngle(self, contained_energy: Union[str, int, float], def calcReducedObservationAngle(self, contained_energy: Union[str, int, float],
jitter_sigma: u.Quantity = None, obstruction: float = 0.0) -> u.Quantity: jitter_sigma: u.Quantity = None, obstruction: float = 0.0) -> u.Quantity:
@ -98,7 +99,7 @@ class Airy(IPSF):
reduced_observation_angle_pixel = (self.__pixel_size / ( reduced_observation_angle_pixel = (self.__pixel_size / (
self.__f_number * self.__d_aperture) * self.__d_aperture / self.__wl).decompose() self.__f_number * self.__d_aperture) * self.__d_aperture / self.__wl).decompose()
# Calculate the width of each grid element # Calculate the width of each grid element
dx = reduced_observation_angle_pixel / self.__osf dx = reduced_observation_angle_pixel.value / self.__osf
# Calculate the necessary number of points on the grid # Calculate the necessary number of points on the grid
n_points = np.ceil(grid_width / dx).value n_points = np.ceil(grid_width / dx).value
# Calculate the corresponding x-coordinates of each grid element # Calculate the corresponding x-coordinates of each grid element
@ -118,14 +119,18 @@ class Airy(IPSF):
psf = fftconvolve(np.pad(psf, int(n_points), mode="constant", constant_values=0), kernel, mode="same") psf = fftconvolve(np.pad(psf, int(n_points), mode="constant", constant_values=0), kernel, mode="same")
# Reduce the PSF to the positive x-domain # Reduce the PSF to the positive x-domain
psf = psf[int((psf.shape[0] - 1) / 2):] psf = psf[int((psf.shape[0] - 1) / 2):]
# Scale the integral of the disturbed PSF equal to the undisturbed PSF
psf = psf / (np.sum(psf * np.arange(psf.shape[0]) * dx) * dx * 2 * np.pi) * total
self.__psf_jitter = psf
if isinstance(contained_energy, str) and contained_energy.lower() == "fwhm": if isinstance(contained_energy, str) and contained_energy.lower() == "fwhm":
reduced_observation_angle = np.argmax(psf < psf[0] / 2) * reduced_observation_angle_pixel.value /\ reduced_observation_angle = np.argmax(psf < psf[0] / 2) * reduced_observation_angle_pixel.value / \
self.__osf * 2 self.__osf * 2
else: else:
# Calculate the rolling integral of the PSF # Calculate the rolling integral of the PSF
psf_int = np.cumsum(psf * np.arange(psf.shape[0])) * reduced_observation_angle_pixel.value / self.__osf psf_int = np.cumsum(psf * np.arange(psf.shape[0]) * dx) * dx * 2 * np.pi
# Scale the integral of the disturbed PSF equal to the undisturbed PSF # Scale the integral of the disturbed PSF equal to the undisturbed PSF
psf_int = psf_int / psf_int[-1] * total / (4 / np.pi) psf_int = psf_int / (4 / np.pi) * (1 - obstruction) ** 2
# Calculate the reduced observation angle # Calculate the reduced observation angle
reduced_observation_angle = np.argmax( reduced_observation_angle = np.argmax(
psf_int > contained_energy) * reduced_observation_angle_pixel.value / self.__osf * 2 psf_int > contained_energy) * reduced_observation_angle_pixel.value / self.__osf * 2

4
tests/psf/test_Airy.py
View File

@ -29,5 +29,5 @@ class TestAiry(TestCase):
# Jitter, obstructed # Jitter, obstructed
self.assertAlmostEqual(self.airy.calcReducedObservationAngle("peak", 1 * u.arcsec, 0.04).value, 0.0) self.assertAlmostEqual(self.airy.calcReducedObservationAngle("peak", 1 * u.arcsec, 0.04).value, 0.0)
self.assertAlmostEqual(self.airy.calcReducedObservationAngle("fwhm", 1 * u.arcsec, 0.04).value, 1.725) self.assertAlmostEqual(self.airy.calcReducedObservationAngle("fwhm", 1 * u.arcsec, 0.04).value, 1.725)
self.assertAlmostEqual(self.airy.calcReducedObservationAngle("min", 1 * u.arcsec, 0.04).value, 2.875) self.assertAlmostEqual(self.airy.calcReducedObservationAngle("min", 1 * u.arcsec, 0.04).value, 3.325)
self.assertAlmostEqual(self.airy.calcReducedObservationAngle(80., 1 * u.arcsec, 0.04).value, 3.1) self.assertAlmostEqual(self.airy.calcReducedObservationAngle(80., 1 * u.arcsec, 0.04).value, 3.7)