text: import numpy as np
  from scipy.integrate import quad
  import astropy.units as u
  from utils import _get_cosmology
  
  
  class Cosmology:
      def __init__(self, cosmo_name='Planck18'):
          self.cosmo_name = cosmo_name
          self.cosmology = self._get_cosmology()
          self.H0 = self.cosmology.H0
          self.h = self.cosmology.h
          self.t_H = self.cosmology.hubble_time
          self.D_H = self.cosmology.hubble_distance
          self.Omega_M = self.cosmology.Om0
          self.Omega_l = self.cosmology.Ode0
          self.Omega_k = 1 - self.Omega_M - self.Omega_l
  
      def _get_cosmology(self):
          """
          Import a built-in cosmology from Astropy
          """
          if self.cosmo_name == 'WMAP5':
              from astropy.cosmology import WMAP5 as cosmo
          elif self.cosmo_name == 'WMAP7':
              from astropy.cosmology import WMAP7 as cosmo
          elif self.cosmo_name == 'WMAP9':
              from astropy.cosmology import WMAP9 as cosmo
          elif self.cosmo_name == 'Planck13':
              from astropy.cosmology import Planck13 as cosmo
          elif self.cosmo_name == 'Planck15':
              from astropy.cosmology import Planck15 as cosmo
          elif self.cosmo_name == 'Planck18':
              from astropy.cosmology import Planck18 as cosmo
          else:
              raise ValueError(f'Cosmology "{self.cosmo_name}" not recognized,' +
                               ' must be "WMAP5", "WMAP7", "WMAP9", ' +
                               '"Planck13", "Planck15", or "Planck18"')
  
          return cosmo
  
      def _E(self, z):
          """
          A convenient function proportional to the time derivative of the
          logarithm of the scale factor, E(z) = H(z) / H0
  
          Parameters
          ----------
          z : int or float
              The redshift at which to calculate the comoving radial distance
  
          Returns
          -------
          E_z : float
              The value of the function E(z)
          """
          return np.sqrt(self.Omega_M * (1 + z)**3 +
                         self.Omega_k * (1 + z)**2 +
                         self.Omega_l)
  
      def comoving_radial_distance(self, z):
          """
          The comoving radial distance at redshift z
  
          Parameters
          ----------
          z : int or float
              The redshift at which to calculate the comoving radial distance
  
          Returns
          -------
          D_C : astropy.Quantity
              The comoving radial distance
          """
          def integrand(z):
              return self.D_H.value / self._E(z)
  
          D_C = quad(integrand, 0, z)[0] * u.Mpc
  
          return D_C
  
      def comoving_transverse_distance(self, z):
          """
          The comoving transverse distance at redshift z
  
          Parameters
          ----------
          z : int or float
              The redshift at which to calculate the comoving transverse distance
  
          Returns
          -------
          D_M : astropy.Quantity
              The comoving transverse distance
          """
          D_C = self.comoving_radial_distance(z)
  
          if self.Omega_l == 0:
              D_M = (2 * self.D_H / (self.Omega_M**2 * (1 + z)) *
                    (2 - self.Omega_M * (1 - z) - (2 - self.Omega_M) *  # noqa
                    np.sqrt(1 + self.Omega_M * z)))  # noqa
          else:
              if self.Omega_k == 0:
                  D_M = D_C
              elif self.Omega_k > 0:
                  sinh_arg = (np.sqrt(self.Omega_k) * D_C / self.D_H) * u.rad
                  D_M = self.D_H / np.sqrt(self.Omega_k) * np.sinh(sinh_arg)
              elif self.Omega_k < 0:
                  sin_arg = (np.sqrt(np.abs(self.Omega_k)) *
                             D_C / self.D_H * u.rad)
                  D_M = (self.D_H / np.sqrt(np.abs(self.Omega_k)) *
                         np.sin(sin_arg))
  
          return D_M
  
      def angular_diameter_distance(self, z):
          """
          The angular diameter distance at redshift z
  
          Parameters
          ----------
          z : int or float
              The redshift at which to calculate the angular diameter distance
  
          Returns
          -------
          D_A : astropy.Quantity
              The angular diameter distance
          """
          D_M = self.comoving_transverse_distance(z)
          D_A = D_M / (1 + z)
  
          return D_A
  
      def luminosity_distance(self, z):
          """
          The luminosity distance at redshift z
  
          Parameters
          ----------
          z : int or float
              The redshift at which to calculate the luminosity distance
  
          Returns
          -------
          D_L : astropy.Quantity
              The luminosity distance
          """
          D_L = (1 + z) * self.comoving_transverse_distance(z)
  
          return D_L
  
      def comoving_volume_in_area(self, z_min, z_max, area):
          """
          The comoving volume within a given area at redshift z
  
          Parameters
          ----------
          z_min : int or float
              The lower bound of the redshift interval in which to calculate the
              comoving volume
          z_max : int or float
              The upper bound of the redshift interval in which to calculate the
              comoving volume
          area : astropy.Quantity
              The area in which to calculate the comoving volume, must be an
              astropy.Quantity of area (e.g., arcsec^2, deg^2, sr)
  
          Returns
          -------
          V_C : astropy.Quantity
              The comoving volume in the area
          """
          def integrand(z):
              D_A = self.angular_diameter_distance(z)
              dV_C_dz = (self.D_H * (1 + z)**2 * D_A**2 * area.to(u.sr) /
                         self._E(z))
              return dV_C_dz.value
  
          V_C = quad(integrand, z_min, z_max)[0] * u.Mpc**3
  
          return V_C
  
      def convert_m_to_M(self, m, z):
          """
          Convert from apparent to absolute magnitude at redshift z
  
          Parameters
          ----------
          m : float
              The apparent magnitude
          z : int or float
              The redshift at which to convert
  
          Returns
          -------
          M : astropy.Quantity
              The absolute magnitude
          """
          D_L = self.luminosity_distance(z) * 1e6  # pc
          if hasattr(m, 'unit'):
              m = m.value
          M = (m - 5 * (np.log10(D_L.value) - 1) +
               2.5 * np.log10(1 + z)) * u.ABmag
  
          return M
  
      def convert_M_to_m(self, M, z):
          """
          Convert from absolute to apparent magnitude at redshift z
  
          Parameters
          ----------
          M : float
              The absolute magnitude
          z : int or float
              The redshift at which to convert
  
          Returns
          -------
          m : astropy.Quantity
              The apparent magnitude
          """
          D_L = self.luminosity_distance(z) * 1e6  # pc
          if hasattr(M, 'unit'):
              M = M.value
          m = (M + 5 * (np.log10(D_L.value) - 1) -
               2.5 * np.log10(1 + z)) * u.ABmag
  
          return m
  