Minor PSF refactor and more logs (#901)

Author: teutoburg

scopesim/effects/psfs/psf_base.py

@@ -49,7 +49,7 @@ def __init__(self, **kwargs):
             "wave_key": "WAVE0",
             "normalise_kernel": True,
             "rounded_edges": True,
-            "rotational_blur_angle": 0,
+            "rotational_blur_angle": 0*u.deg,
         }
         self.meta.update(params)
         self.meta.update(kwargs)
@@ -76,7 +76,7 @@ def apply_to(self, obj, **kwargs):
 
                 # apply rotational blur for field-tracking observations
                 rot_blur_angle = self.meta["rotational_blur_angle"]
-                if abs(rot_blur_angle) > 0:
+                if abs(rot_blur_angle << u.deg) > 0*u.deg:
                     # makes a copy of kernel
                     kernel = rotational_blur(kernel, rot_blur_angle)
 
@@ -150,35 +150,41 @@ def plot(self, obj=None, **kwargs):
         return fig
 
 
-def rotational_blur(image, angle):
+@u.quantity_input
+def rotational_blur(image, angle: u.Quantity[u.deg]):
     """
     Rotate and coadd an image over a given angle to imitate a blur.
 
     Parameters
     ----------
-    image : array
-        Image to blur
-    angle : float
-        [deg] Angle over which the image should be rotationally blurred
+    image : array-like
+        Image to blur.
+    angle : u.Quantity["angle"]
+        Angle over which the image should be rotationally blurred.
+
+    .. versionchanged:: PLACEHOLDER_NEXT_RELEASE_VERSION
+
+       Require `angle` to be a Quantity.
 
     Returns
     -------
-    image_rot : array
+    image_rot : np.ndarray
         Blurred image
 
     """
     image_rot = np.copy(image)
 
-    n_angles = 1
-    rad_to_deg = 57.29578
-    edge_pixel_unit_angle = np.arctan2(1, (image.shape[0] // 2)) * rad_to_deg
-    while abs(angle) > edge_pixel_unit_angle and n_angles < 25:
-        angle /= 2.
-        image_rot += rotate(image_rot, angle, reshape=False, order=1)
+    edge_pixel_unit_angle = np.arctan2(1, (image.shape[0] // 2)) * u.rad
+    n_steps = np.ceil(np.log2(abs(angle) / edge_pixel_unit_angle))
+    n_steps = int(min(n_steps, 8))  # avoid overrun
+
+    current_angle = angle.copy()
+    for _ in range(n_steps):
+        current_angle /= 2.
+        image_rot += rotate(image_rot, current_angle, reshape=False, order=3)
         # each time kernel is rotated and added, the frame total doubles
-        n_angles *= 2
 
-    return image_rot / n_angles
+    return image_rot / image_rot.sum() * image.sum()
 
 
 def get_bkg_level(obj, bg_w):

scopesim/effects/psfs/semianalytical.py

@@ -11,10 +11,18 @@
 
 from .. import ter_curves_utils as tu
 from ...optics.fov import FieldOfView
-from ...utils import (figure_factory, figure_grid_factory, from_currsys,
-                      quantify, check_keys)
+from ...utils import (
+    figure_factory,
+    figure_grid_factory,
+    from_currsys,
+    quantify,
+    check_keys,
+    get_logger,
+)
 from . import PSF
 
+logger = get_logger(__name__)
+
 
 class SemiAnalyticalPSF(PSF):
     """Base class for semianalytical PSFs."""
@@ -115,9 +123,11 @@ def get_kernel(self, fov):
 
         n = self.meta["psf_side_length"]
         wave = self.wavelength
+        logger.debug("    start AnisoCADO")
         self._psf_object = aniso.AnalyticalScaoPsf(pixelSize=pixel_scale,
                                                    N=n, wavelength=wave,
                                                    nmRms=self.nmRms)
+        logger.debug("    done AnisoCADO")
         if np.any(self.meta["offset"]):
             self._psf_object.shift_off_axis(self.meta["offset"][0],
                                             self.meta["offset"][1])

scopesim/tests/tests_effects/test_PSF.py

@@ -3,6 +3,7 @@
 from pytest import approx
 import numpy as np
 from matplotlib import pyplot as plt
+from astropy import units as u
 
 from scopesim.effects import PSF
 from scopesim.effects.psfs.psf_base import get_bkg_level
@@ -56,7 +57,7 @@ def test_rescale_produces_odd_size_kernel(self, nkern, scale_factor):
         assert n_x % 2 == 1
 
 class TestRotationBlur:
-    @pytest.mark.parametrize("angle", ([1, 5, 15, 60]))
+    @pytest.mark.parametrize("angle", ([1, 5, 15, 60]*u.deg))
     def test_returns_rotated_kernel_array_has_same_sum(self, angle):
         # Without blur
         implane = basic_image_plane()
@@ -169,7 +170,7 @@ def test_convolves_with_3D_cube(self):
             plt.subplot(231)
             plt.imshow(implane.data[1, :, :])
 
-        psf = PSF(rotational_blur_angle=15, bkg_width=5, rounded_edges=False)
+        psf = PSF(rotational_blur_angle=15*u.deg, bkg_width=5, rounded_edges=False)
         psf.kernel = basic_kernel(n=63)
         implane = psf.apply_to(implane)