import copy


def continuous_injection_model(energy_gev, k0, index, Ebreak,
                               E0=1.0):
    """
    Compute a CI spectrum. See Pacholczyk 1970,
    and Turner 2017

    Parameters
    ----------
    - energy_GeV: scalar or vector
    - k0 : normalization
    - index : spectral index
    - E0 : pivot energy (GeV)
    - Ebreak : break energy (GeV)

    Outputs
    --------
    - spectrum
    """

    energy_gev_1 = copy.copy(energy_gev)
    energy_gev_2 = copy.copy(energy_gev)
    w1 = energy_gev >= Ebreak
    w2 = energy_gev < Ebreak
    energy_gev_2[w1] = Ebreak/2.0

    f1 = k0 * (energy_gev_1/E0)**(-index-1)
    f2 = 1 - (1-energy_gev_2/Ebreak)**(index-1)
    S = energy_gev*0
    S[w1] = f1[w1]
    S[w2] = (f1*f2)[w2]
    
    return S