Don't error if two photon_arrays being convolved have each set ancill…

…ary arrays

CHANGELOG.rst

@@ -64,3 +64,5 @@ Bug Fixes
 - Fixed some errors when drawing ChromaticTransformation objects with photon shooting. (#1229)
 - Fixed the flux drawn by ChromaticConvolution with photon shooting when poisson_flux=True. (#1229)
 - Fixed a slight inaccuracy in the FFT phase shifts for single-precision images. (#1231, #1234)
+- Fixed a bug that prevented a convolution of two PhaseScreenPSF objects from being drawn with
+  photon shooting. (#1242)

galsim/photon_array.py

@@ -465,44 +465,28 @@ def assignAt(self, istart, rhs):
 
     def convolve(self, rhs, rng=None):
         """Convolve this `PhotonArray` with another.
+
+        ..note::
+
+            If both self and rhs have wavelengths, angles, pupil coordinates or times assigned,
+            then the values from the first array (i.e. self) take precedence.
         """
         if rhs.size() != self.size():
             raise GalSimIncompatibleValuesError("PhotonArray.convolve with unequal size arrays",
                                                 self_pa=self, rhs=rhs)
-        if rhs.hasAllocatedAngles():
-            if self.hasAllocatedAngles():
-                raise GalSimIncompatibleValuesError(
-                    "PhotonArray.convolve with doubly assigned angles"
-                )
-            else:
-                self.dxdz = rhs.dxdz
-                self.dydz = rhs.dydz
+        if rhs.hasAllocatedAngles() and not self.hasAllocatedAngles():
+            self.dxdz = rhs.dxdz
+            self.dydz = rhs.dydz
 
-        if rhs.hasAllocatedWavelengths():
-            if self.hasAllocatedWavelengths() and self._wave is not rhs._wave:
-                raise GalSimIncompatibleValuesError(
-                    "PhotonArray.convolve with doubly assigned wavelengths"
-                )
-            else:
-                self.wavelength = rhs.wavelength
+        if rhs.hasAllocatedWavelengths() and not self.hasAllocatedWavelengths():
+            self.wavelength = rhs.wavelength
 
-        if rhs.hasAllocatedPupil():
-            if self.hasAllocatedPupil() and (
-                    self._pupil_u is not rhs._pupil_u or self._pupil_v is not rhs._pupil_v):
-                raise GalSimIncompatibleValuesError(
-                    "PhotonArray.convolve with doubly assigned pupil coordinates"
-                )
-            else:
-                self.pupil_u = rhs.pupil_u
-                self.pupil_v = rhs.pupil_v
+        if rhs.hasAllocatedPupil() and not self.hasAllocatedPupil():
+            self.pupil_u = rhs.pupil_u
+            self.pupil_v = rhs.pupil_v
 
-        if rhs.hasAllocatedTimes():
-            if self.hasAllocatedTimes() and self._time is not rhs._time:
-                raise GalSimIncompatibleValuesError(
-                    "PhotonArray.convolve with doubly assigned time stamps"
-                )
-            else:
-                self.time = rhs.time
+        if rhs.hasAllocatedTimes() and not self.hasAllocatedTimes():
+            self.time = rhs.time
 
         rng = BaseDeviate(rng)
         self._pa.convolve(rhs._pa, rng._rng)

tests/test_phase_psf.py

@@ -1521,6 +1521,7 @@ def test_uv_persistence():
         check_pickle(photons)
 
 
+@timer
 def test_t_persistence():
     rng = galsim.BaseDeviate(10)
     atm = galsim.Atmosphere(screen_size=10.0, altitude=[0,1,2,3], r0_500=0.2, rng=rng)
@@ -1531,6 +1532,21 @@ def test_t_persistence():
     assert np.min(photons.time) > 10.0
     assert np.max(photons.time) < 25.0
 
+@timer
+def test_convolve_phasepsf():
+    # This snippet didn't use to be allowed, since the two psfs populate different
+    # pupil angles.
+    star = galsim.DeltaFunction(flux=10)
+    psf1 = galsim.Atmosphere(screen_size=10).makePSF(lam=700, diam=8)
+    psf2 = galsim.OpticalPSF(lam=700, diam=8, geometric_shooting=True)
+    sed = galsim.SED('vega.txt', wave_type='nm', flux_type='fphotons')
+    obj = galsim.Convolve(star * sed, psf1, psf2)
+    bandpass = galsim.Bandpass('LSST_r.dat', wave_type='nm')
+    obj = obj.withFlux(10, bandpass)
+    im = obj.drawImage(bandpass, method='phot', n_photons=10)
+
+    # The main thing is that it works.  But check that flux makes sense.
+    assert im.array.sum() == 10
 
 if __name__ == "__main__":
     testfns = [v for k, v in vars().items() if k[:5] == 'test_' and callable(v)]

tests/test_photon_array.py

@@ -365,8 +365,15 @@ def test_convolve():
         np.testing.assert_array_equal(getattr(pa2, attr), data)
 
         # both have data now...
-        with assert_raises(galsim.GalSimIncompatibleValuesError):
-            pa1.convolve(pa2)
+        pa1.convolve(pa2)
+        np.testing.assert_array_equal(getattr(pa1, attr), data)
+        np.testing.assert_array_equal(getattr(pa2, attr), data)
+
+        # If the second one has different data, the first takes precedence.
+        setattr(pa2, attr, data * 2)
+        pa1.convolve(pa2)
+        np.testing.assert_array_equal(getattr(pa1, attr), data)
+        np.testing.assert_array_equal(getattr(pa2, attr), 2*data)
 
 
 @timer