How to use the lenstronomy.Util.constants function in lenstronomy

anisotropy_model = 'OM'
        kwargs_numerics_galkin = {'interpol_grid_num': 500, 'log_integration': True,
                                  'max_integrate': 10, 'min_integrate': 0.001}
        self.td_cosmo.kinematics_modeling_settings(anisotropy_model, kwargs_numerics_galkin, analytic_kinematics=True,
                                             Hernquist_approx=False, MGE_light=False, MGE_mass=False)

        J = self.td_cosmo.velocity_dispersion_dimension_less(self.kwargs_lens, kwargs_lens_light, kwargs_anisotropy, r_eff=r_eff,
                                           theta_E=self.kwargs_lens[0]['theta_E'], gamma=2)

        J_map = self.td_cosmo.velocity_dispersion_map_dimension_less(self.kwargs_lens, kwargs_lens_light,
                                                                     kwargs_anisotropy, r_eff=r_eff,
                                                                     theta_E=self.kwargs_lens[0]['theta_E'], gamma=2)
        assert len(J_map) == 1
        npt.assert_almost_equal(J_map[0]/J, 1, decimal=1)
        sigma_v2 = J * D_s/D_ds * const.c ** 2
        sigma_v = np.sqrt(sigma_v2) / 1000.  # convert to [km/s]
        print(sigma_v, 'test sigma_v')
        Ds_Dds = self.td_cosmo.ds_dds_from_kinematics(sigma_v, J, kappa_s=0, kappa_ds=0)
        npt.assert_almost_equal(Ds_Dds, D_s/D_ds)

        # now we perform a mass-sheet transform in the observables but leave the models identical with a convergence correction
        kappa_s = 0.5
        dt_list = self.td_cosmo.time_delays(self.kwargs_lens, self.kwargs_ps, kappa_ext=kappa_s)
        sigma_v_kappa = sigma_v * np.sqrt(1-kappa_s)
        dt = dt_list[0] - dt_list[1]
        D_dt_infered, D_d_infered = self.td_cosmo.ddt_dd_from_time_delay_and_kinematics(d_fermat_model=d_fermat, dt_measured=dt,
                                                                                        sigma_v_measured=sigma_v_kappa, J=J, kappa_s=kappa_s,
                                                                                        kappa_ds=0, kappa_d=0)
        npt.assert_almost_equal(D_dt_infered, D_dt, decimal=6)
        npt.assert_almost_equal(D_d_infered, D_d, decimal=6)

def test_sis_travel_time(self):
        z_source = 1.5
        z_lens = 0.5
        lens_model_list = ['SIS']
        redshift_list = [z_lens]
        lensModelMutli = MultiPlane(z_source=z_source, lens_model_list=lens_model_list, lens_redshift_list=redshift_list)
        lensModel = LensModel(lens_model_list=lens_model_list)
        kwargs_lens = [{'theta_E': 1., 'center_x': 0, 'center_y': 0}]
        dt = lensModelMutli.arrival_time(1., 0., kwargs_lens)
        Dt = lensModelMutli._multi_plane_base._cosmo_bkg.ddt(z_lens=z_lens, z_source=z_source)
        fermat_pot = lensModel.fermat_potential(1, 0., kwargs_lens)
        dt_simple = const.delay_arcsec2days(fermat_pot, Dt)
        print(dt, dt_simple)
        npt.assert_almost_equal(dt, dt_simple, decimal=8)

def epsilon_crit(self):
        """
        returns the critical projected mass density in units of M_sun/Mpc^2 (physical units)
        """
        if not hasattr(self,'Epsilon_Crit'):
            const_SI = constants.c ** 2 / (4 * np.pi * constants.G)  #c^2/(4*pi*G) in units of [kg/m]
            conversion = constants.Mpc / constants.M_sun  # converts [kg/m] to [M_sun/Mpc]
            const = const_SI*conversion   #c^2/(4*pi*G) in units of [M_sun/Mpc]
            self.Epsilon_Crit = self.dist_OS/(self.dist_OL*self.dist_LS) * const #[M_sun/Mpc^2]
        return self.Epsilon_Crit

def mass_3d_interp(self, r, kwargs, new_compute=False):
        """

        :param r: in arc seconds
        :param kwargs: lens model parameters in arc seconds
        :param new_compute: bool, if True, recomputes the interpolation
        :return: mass enclosed physical radius in kg
        """
        if not hasattr(self, '_log_mass_3d') or new_compute is True:
            r_array = np.logspace(np.log10(self._min_interpolate), np.log10(self._max_interpolate), self._interp_grid_num)
            mass_3d_array = self.model.mass_3d(r_array, kwargs)
            mass_3d_array[mass_3d_array < 10. ** (-10)] = 10. ** (-10)
            mass_dim_array = mass_3d_array * const.arcsec ** 2 * self.cosmo.dd * self.cosmo.ds \
                             / self.cosmo.dds * const.Mpc * const.c ** 2 / (4 * np.pi * const.G)
            f = interp1d(np.log(r_array), np.log(mass_dim_array/r_array), fill_value="extrapolate")
            self._log_mass_3d = f
        return np.exp(self._log_mass_3d(np.log(r))) * r

l(r) \sigma_r(r) ^ 2 =  1/f(r) \int_r^{\infty} f(s) l(s) G M(s) / s^2 ds
        where l(r) is the 3d light profile
        M(s) is the enclosed 3d mass
        f is the solution to
        d ln(f)/ d ln(r) = 2 beta(r)

        :param r: 3d radius
        :param kwargs_mass: mass model parameters (following lenstronomy lens model conventions)
        :param kwargs_light: deflector light parameters (following lenstronomy light model conventions)
        :param kwargs_anisotropy: anisotropy parameters, may vary according to anisotropy type chosen.
            We refer to the Anisotropy() class for details on the parameters.
        :return: sigma_r**2
        """
        l_r = self.lightProfile.light_3d_interp(r, kwargs_light)
        f_r = self.anisotropy_solution(r, **kwargs_anisotropy)
        return 1 / f_r / l_r * self._jeans_solution_integral(r, kwargs_mass, kwargs_light, kwargs_anisotropy) * const.G / (const.arcsec * self.cosmo.dd * const.Mpc)

def _geometrical_delay(beta_i_x, beta_i_y, beta_j_x, beta_j_y, T_i, T_j, T_ij):
        """

        :param beta_i_x: angle on the sky at plane i
        :param beta_i_y: angle on the sky at plane i
        :param beta_j_x: angle on the sky at plane j
        :param beta_j_y: angle on the sky at plane j
        :param T_i: transverse diameter distance to z_i
        :param T_j: transverse diameter distance to z_j
        :param T_ij: transverse diameter distance from z_i to z_j
        :return: excess delay relative to a straight line
        """
        d_beta_x = beta_j_x - beta_i_x
        d_beta_y = beta_j_y - beta_i_y
        tau_ij = T_i * T_j / T_ij * const.Mpc / const.c / const.day_s * const.arcsec**2
        return tau_ij * (d_beta_x ** 2 + d_beta_y ** 2) / 2

fermat potential (negative sign means earlier arrival time)
        for Multi-plane lensing, it computes the effective Fermat potential (derived from the arrival time and
        subtracted off the time-delay distance for the given cosmology). The units are given in arcsecond square.

        :param x_image: image position
        :param y_image: image position
        :param x_source: source position
        :param y_source: source position
        :param kwargs_lens: list of keyword arguments of lens model parameters matching the lens model classes
        :return: fermat potential in arcsec**2 without geometry term (second part of Eqn 1 in Suyu et al. 2013) as a list
        """
        if hasattr(self.lens_model, 'fermat_potential'):
            return self.lens_model.fermat_potential(x_image, y_image, kwargs_lens, x_source, y_source)
        elif hasattr(self.lens_model, 'arrival_time') and hasattr(self, '_lensCosmo'):
            dt = self.lens_model.arrival_time(x_image, y_image, kwargs_lens)
            fermat_pot_eff = dt * const.c / self._lensCosmo.ddt / const.Mpc * const.day_s / const.arcsec ** 2
            return fermat_pot_eff
        else:
            raise ValueError('In multi-plane lensing you need to provide a specific z_lens and z_source for which the effective Fermat potential is evaluated')

def days_D_model(self, delay_arcsec, D_dt_model):
        """
        given a delay in arcsec^2 and a Delay distance, the delay is computed in days

        :param delay_arc_sec:
        :param D_dt_model:
        :return:
        """
        D_dt = D_dt_model * const.Mpc  # eqn 7 in Suyu et al.
        return D_dt / const.c * delay_arcsec / const.day_s * const.arcsec**2  # * self.arcsec2phys_lens(1.)**2