Update atmPSF test

tests/test_atmPSF.py

@@ -1,23 +1,29 @@
 import unittest
 
 import numpy as np
+import galsim
 
 from desc.imsim.atmPSF import AtmosphericPSF
 
 class AtmPSF(unittest.TestCase):
     def test_r0_500(self):
-        """Test that _seeing_resid has the API fsolve wants."""
-        for wavelength in [300.0, 500.0, 1100.0]:
-            for L0 in [10.0, 25.0, 100.0]:
-                for target_seeing in [0.5, 0.7, 1.0]:
-                    r0s = [0.1, 0.2]
-                    np.testing.assert_array_equal(
-                        np.hstack([
-                            AtmosphericPSF._seeing_resid(r0s[0], wavelength, L0, target_seeing),
-                            AtmosphericPSF._seeing_resid(r0s[1], wavelength, L0, target_seeing),
-                        ]),
-                        AtmosphericPSF._seeing_resid(r0s, wavelength, L0, target_seeing)
-                    )
+        """Test that inversion of the Tokovinin fitting formula for r0_500 works."""
+        np.random.seed(57721)
+        for _ in range(100):
+            airmass = np.random.uniform(1.001, 1.5)
+            rawSeeing = np.random.uniform(0.5, 1.5)
+            band = 'ugrizy'[np.random.randint(6)]
+            rng = galsim.BaseDeviate(np.random.randint(2**32))
+            atmPSF = AtmosphericPSF(airmass, rawSeeing, band, rng, screen_size=6.4)
+
+            wlen = dict(u=365.49, g=480.03, r=622.20, i=754.06, z=868.21, y=991.66)[band]
+            targetFWHM = rawSeeing * airmass**0.6 * (wlen/500)**(-0.3)
+
+            r0_500 = atmPSF.atm.r0_500_effective
+            L0 = atmPSF.atm[0].L0
+            vkFWHM = AtmosphericPSF._vkSeeing(r0_500, wlen, L0)
+
+            np.testing.assert_allclose(targetFWHM, vkFWHM, atol=1e-3, rtol=0)
 
 
 if __name__ == '__main__':