import common
import numpy as np
from scipy.integrate import cumulative_trapezoid


def d_L(N, Omega_m0, type, Omega_r0=0, Omega_phi0=0, omega_phi_values=0):
    """
    Compute dimensionless luminosity distance: H_0 * d_L / c
    
    Formula: d_L(N) = exp(-N) * กา_0^N dN' [-exp(-N') / (H(N')/H_0)]
    """
    N = N[::-1]
    if type == 'quintessence':
        omega_phi_values = omega_phi_values[::-1]

    integrand = np.exp(-N) / common.H_over_H0(N, Omega_m0, type, Omega_r0, Omega_phi0,
                                                omega_phi_values)
    
    integral_from_start = -cumulative_trapezoid(integrand, N, initial=0)
    
    return (np.exp(-N) * integral_from_start)[::-1]