def get_integral_powerlaw_model(Emin, Emax, k0, index, E0=1.0):
    """
    Compute the enery integral :
    \int_Emin^Emax f(E) dE
    for f(E) a power law

    Parameters
    ----------
    - Emin (GeV): the lower bound
    - Emax (GeV): the upper bound
    - k0 : normalization
    - E0 : pivot energy (GeV)
    - index : spectral index

    Outputs
    --------
    - The integrated function
    
    """

    if Emin > Emax:
        raise TypeError("Emin is larger than Emax")

    output = k0 * E0**(1-index) / (1-index) * ( (Emax/E0)**(1-index) - (Emin/E0)**(1-index))
    
    return output