import astropy.units as u
import numpy as np
from scipy.integrate import cumulative_trapezoid

h = 0.7
G = 6.67430e-11  # m^3 kg^-1 s^-2
c = 2.99792458

def N_of_a(a, a0=1):
    """
    Returns the value of N for a given a.
    """
    return np.log(a/a0)

def a_of_z(z, a0=1):
    return a0/(1+z)

def z_of_a(a, a0=1):
    return a0/a - 1

def a_of_N(N, a0=1):
    return a0 * np.exp(N)

def omega_phi(x1, x2):
    return t

def H_over_H0(N, Omega_m0, type, Omega_r0=0, Omega_phi0=0, omega_phi_values=0):
    if type == 'quintessence':
        integrand = 3 * (1 + omega_phi_values)
        integral = cumulative_trapezoid(integrand, N, initial=0)
        return np.sqrt(Omega_m0 * np.exp(-3*N) + Omega_r0 * np.exp(-4*N) +
                        Omega_phi0 * np.exp(-integral))
    elif type == 'lcdm':
        return np.sqrt(Omega_m0 * np.exp(-3*N) + (1-Omega_m0))