diff --git a/galsim/image.py b/galsim/image.py
index 747c97a46b..a6d20b95bc 100644
--- a/galsim/image.py
+++ b/galsim/image.py
@@ -1590,13 +1590,18 @@ def calculateFWHM(self, center=None, Imax=0.):
 
         # Find the first value with I < 0.5 * Imax
         k = np.argmax(data < 0.5 * Imax)
-        Ik = data[k] / Imax
+        # If there are duplicate rsq values, argmax won't be deterministic across numpy versions.
+        # To achieve consistent results, we take the average of all data values with the same rsq.
+        k2 = np.searchsorted(rsqf, rsqf[k], side='right')
+        Ik = np.mean(data[k:k2]) / Imax
 
         # Interpolate (linearly) between this and the previous value.
         if k == 0:
             rsqhm = rsqf[0] * (0.5 / Ik)
         else:
-            Ikm1 = data[k-1] / Imax
+            # Similarly, use the mean data for all rsq equal to rsqf[k-1].
+            k1 = np.searchsorted(rsqf, rsqf[k-1], side='left')
+            Ikm1 = np.mean(data[k1:k]) / Imax
             rsqhm = (rsqf[k-1] * (Ik-0.5) + rsqf[k] * (0.5-Ikm1)) / (Ik-Ikm1)
 
         # This has all been done in pixels.  So normalize according to the pixel scale.
@@ -1789,7 +1794,10 @@ def __pow__(self, other):
     def __ipow__(self, other):
         if not isinstance(other, int) and not isinstance(other, float):
             raise TypeError("Can only raise an image to a float or int power!")
-        self.array[:,:] **= other
+        if not self.isinteger or isinstance(other, int):
+            self.array[:,:] **= other
+        else:
+            self.array[:,:] = self._safe_cast(self.array ** other)
         return self
 
     def __neg__(self):
diff --git a/tests/test_calc.py b/tests/test_calc.py
index a80e82c71f..655e8b683f 100644
--- a/tests/test_calc.py
+++ b/tests/test_calc.py
@@ -369,5 +369,30 @@ def test_fwhm():
                                    err_msg="non-square image.calculateFWHM is not accurate.")
 
 
+def test_fwhm_deterministic_fwhm():
+    # This test was proposed by Wiliam Lucas in issue #1336.
+    # Numpy changed some details about how argsort works, so different numpy versions had
+    # given different answers to the following calculation.
+
+    gal = galsim.Gaussian(flux=1.e5, sigma=5.)
+    rng = galsim.BaseDeviate(1234)
+    img = gal.drawImage(method='phot', rng=rng, nx=50, ny=50)
+    fwhm = img.calculateFWHM()
+    # This value used to be different on systems with different numpy versions.
+    # Note: this only affects situations where the center is a pixel center or pixel corner
+    # (or some other very special locations) such that the rsq values include repeated values.
+    # Then argsort could return a different data index for what pixel was the first one past
+    # the half max value.  For symmetric profiles, even this probalby wasn't a problem, since all
+    # data values at that radius have the same value.  But for photon shooting, there is photon
+    # noise, which means the values are not identical, leading to indeterminacy in the FWHM value.
+    print(fwhm)
+
+    # The code was changed in GalSim version 2.7 to be insensitive to the details of what argsort
+    # does with repeated values.  Specifically, we now take the mean of all data values with the
+    # same rsq value.
+    # Thus, this value should be consistent across any algorithmic changes numpy might make to
+    # how argsort works in detail.
+    assert np.isclose(fwhm, 13.715780373024351)
+
 if __name__ == "__main__":
     runtests(__file__)