import numpy as np
from scipy.integrate import simpson


def t0H0(HoverH0, N):
    integrand = 1 / HoverH0
    return simpson(integrand, x=N)

def t0H0_lcdm(Omega_m0):
    return 2 / (3 * (1 - Omega_m0) ** 0.5) * np.arcsinh(((1 - Omega_m0) / Omega_m0) ** 0.5)