Merge pull request astropy#15224 from nstarman/auto-backport-of-pr-15…

…049-on-v5.3.x

Backport PR astropy#15049: fix sign of w0wzCDM inv_efuncs pyx (v5.3.x)

astropy/cosmology/flrw/scalar_inv_efuncs.pyx

@@ -206,14 +206,14 @@ def fwpwacdm_inv_efunc(double z, double Om0, double Ode0,
 def w0wzcdm_inv_efunc_norel(double z, double Om0, double Ode0, double Ok0,
     double w0, double wz):
   cdef double opz = 1.0 + z
-  cdef Odescl = opz**(3. * (1. + w0 - wz)) * exp(-3. * wz * z)
+  cdef Odescl = opz**(3. * (1. + w0 - wz)) * exp(3. * wz * z)
   return pow(opz**2 * (opz * Om0 + Ok0) + Ode0 * Odescl, -0.5)
 
 # Massless neutrinos
 def w0wzcdm_inv_efunc_nomnu(double z, double Om0, double Ode0, double Ok0,
     double Or0, double w0, double wz):
   cdef double opz = 1.0 + z
-  cdef Odescl = opz**(3. * (1. + w0 - wz)) * exp(-3. * wz * z)
+  cdef Odescl = opz**(3. * (1. + w0 - wz)) * exp(3. * wz * z)
   return pow((((opz * Or0 + Om0) * opz) + Ok0) * opz**2 +
           Ode0 * Odescl, -0.5)
 
@@ -224,22 +224,22 @@ def w0wzcdm_inv_efunc(double z, double Om0, double Ode0, double Ok0,
 
   cdef double opz = 1.0 + z
   cdef double Or0 = Ogamma0 * (1.0 + nufunc(opz, NeffPerNu, nmasslessnu, nu_y))
-  cdef Odescl = opz**(3. * (1. + w0 - wz)) * exp(-3. * wz * z)
+  cdef Odescl = opz**(3. * (1. + w0 - wz)) * exp(3. * wz * z)
   return pow((((opz * Or0 + Om0) * opz) + Ok0) * opz**2 + Ode0 * Odescl, -0.5)
 
 ######## Flatw0wzCDM
 # No relativistic species
 def fw0wzcdm_inv_efunc_norel(double z, double Om0, double Ode0,
     double w0, double wz):
   cdef double opz = 1.0 + z
-  cdef Odescl = opz**(3. * (1. + w0 - wz)) * exp(-3. * wz * z)
+  cdef Odescl = opz**(3. * (1. + w0 - wz)) * exp(3. * wz * z)
   return pow(opz**3 * Om0 + Ode0 * Odescl, -0.5)
 
 # Massless neutrinos
 def fw0wzcdm_inv_efunc_nomnu(double z, double Om0, double Ode0,
     double Or0, double w0, double wz):
   cdef double opz = 1.0 + z
-  cdef Odescl = opz**(3. * (1. + w0 - wz)) * exp(-3. * wz * z)
+  cdef Odescl = opz**(3. * (1. + w0 - wz)) * exp(3. * wz * z)
   return pow((opz * Or0 + Om0) * opz**3 + Ode0 * Odescl, -0.5)
 
 # With massive neutrinos
@@ -249,7 +249,7 @@ def fw0wzcdm_inv_efunc(double z, double Om0, double Ode0,
 
   cdef double opz = 1.0 + z
   cdef double Or0 = Ogamma0 * (1.0 + nufunc(opz, NeffPerNu, nmasslessnu, nu_y))
-  cdef Odescl = opz**(3. * (1. + w0 - wz)) * exp(-3. * wz * z)
+  cdef Odescl = opz**(3. * (1. + w0 - wz)) * exp(3. * wz * z)
   return pow((opz * Or0 + Om0) * opz**3 + Ode0 * Odescl, -0.5)
 
 ######## Neutrino relative density function

astropy/cosmology/flrw/tests/test_w0wzcdm.py

@@ -132,6 +132,14 @@ def test_Otot_overflow(self, cosmo):
         with pytest.warns(RuntimeWarning, match="overflow encountered in exp"):
             cosmo.Otot(1e3)
 
+    # ===============================================================
+    # I/O Tests
+
+    @pytest.mark.filterwarnings("ignore:overflow encountered")
+    def test_toformat_model(self, cosmo, to_format, method_name):
+        """Test cosmology -> astropy.model."""
+        super().test_toformat_model(cosmo, to_format, method_name)
+
     # ===============================================================
     # Usage Tests
 
@@ -142,17 +150,17 @@ def test_Otot_overflow(self, cosmo):
             (  # no relativistic species
                 (75.0, 0.3, 0.6),
                 {},
-                [3051.68786716, 4756.17714818, 5822.38084257, 6562.70873734] * u.Mpc,
+                [2934.20187523, 4559.94636182, 5590.71080419, 6312.66783729] * u.Mpc,
             ),
             (  # massless neutrinos
                 (75.0, 0.25, 0.5),
                 {"Tcmb0": 3.0, "Neff": 3, "m_nu": 0 * u.eV},
-                [2997.8115653, 4686.45599916, 5764.54388557, 6524.17408738] * u.Mpc,
+                [2904.47062713, 4528.59073707, 5575.95892989, 6318.98689566] * u.Mpc,
             ),
             (  # massive neutrinos
                 (75.0, 0.25, 0.5),
                 {"Tcmb0": 3.0, "Neff": 4, "m_nu": 5 * u.eV},
-                [2676.73467639, 3940.57967585, 4686.90810278, 5191.54178243] * u.Mpc,
+                [2613.84726408, 3849.66574595, 4585.51172509, 5085.16795412] * u.Mpc,
             ),
         ],
     )
@@ -166,6 +174,23 @@ def test_comoving_distance_example(self, cosmo_cls, args, kwargs, expected):
             cosmo_cls, args, {**COMOVING_DISTANCE_EXAMPLE_KWARGS, **kwargs}, expected
         )
 
+    @pytest.mark.skipif(not HAS_SCIPY, reason="scipy is not installed")
+    def test_comoving_distance_mathematica(self, cosmo_cls):
+        """Test with Mathematica example.
+
+        This test should be updated as the code changes.
+
+        ::
+            In[1]:= {Om0, w0, wz, H0, c}={0.3,-0.9, 0.2, 70, 299792.458};
+            c/H0 NIntegrate[1/Sqrt[Om0*(1+z)^3+(1-Om0)(1+z)^(3(1+w0-wz)) Exp[3 *wz*z]],{z, 0, 0.5}]
+            Out[1]= 1849.75
+        """
+        assert u.allclose(
+            cosmo_cls(H0=70, Om0=0.3, w0=-0.9, wz=0.2, Ode0=0.7).comoving_distance(0.5),
+            1849.75 * u.Mpc,
+            rtol=1e-4,
+        )
+
 
 class TestFlatw0wzCDM(FlatFLRWMixinTest, Testw0wzCDM):
     """Test :class:`astropy.cosmology.Flatw0wzCDM`."""
@@ -210,17 +235,17 @@ def test_Otot_overflow(self, cosmo):
             (  # no relativistic species
                 (75.0, 0.3),
                 {},
-                [3156.41804372, 4951.19475878, 6064.40591021, 6831.18710042] * u.Mpc,
+                [3004.55645039, 4694.15295565, 5760.90038238, 6504.07869144] * u.Mpc,
             ),
             (  # massless neutrinos
                 (75.0, 0.25),
                 {"Tcmb0": 3.0, "Neff": 3, "m_nu": 0 * u.eV},
-                [3268.38450997, 5205.96494068, 6419.75447923, 7257.77819438] * u.Mpc,
+                [3086.14574034, 4885.09170925, 6035.4563298, 6840.89215656] * u.Mpc,
             ),
             (  # massive neutrinos
                 (75.0, 0.25),
                 {"Tcmb0": 3.0, "Neff": 4, "m_nu": 5 * u.eV},
-                [2536.77159626, 3721.76294016, 4432.3526772, 4917.90352107] * u.Mpc,
+                [2510.44035219, 3683.87910326, 4389.97760294, 4873.33577288] * u.Mpc,
             ),
         ],
     )
@@ -234,6 +259,23 @@ def test_comoving_distance_example(self, cosmo_cls, args, kwargs, expected):
             cosmo_cls, args, {**COMOVING_DISTANCE_EXAMPLE_KWARGS, **kwargs}, expected
         )
 
+    @pytest.mark.skipif(not HAS_SCIPY, reason="scipy is not installed")
+    def test_comoving_distance_mathematica(self, cosmo_cls):
+        """Test with Mathematica example.
+
+        This test should be updated as the code changes.
+
+        ::
+            In[1]:= {Om0, w0, wz, H0, c}={0.3,-0.9, 0.2, 70, 299792.458};
+            c/H0 NIntegrate[1/Sqrt[Om0*(1+z)^3+(1-Om0)(1+z)^(3(1+w0-wz)) Exp[3 *wz*z]],{z, 0, 0.5}]
+            Out[1]= 1849.75
+        """
+        assert u.allclose(
+            cosmo_cls(H0=70, Om0=0.3, w0=-0.9, wz=0.2).comoving_distance(0.5),
+            1849.75 * u.Mpc,
+            rtol=1e-4,
+        )
+
 
 ##############################################################################
 # Miscellaneous

astropy/cosmology/flrw/w0wzcdm.py

@@ -273,7 +273,7 @@ class Flatw0wzCDM(FlatFLRWMixin, w0wzCDM):
     The comoving distance in Mpc at redshift z:
 
     >>> cosmo.comoving_distance(0.5)
-    <Quantity 1982.66012926 Mpc>
+    <Quantity 1849.74726272 Mpc>
     """
 
     def __init__(

docs/changes/cosmology/15224.bugfix.rst

@@ -0,0 +1 @@
+The exponent of ``w0wzCDM`` functions in ``inv_efunc`` has been corrected to 3, from -3.