def ccm_extinction(wave, ebv, r_v=3.1):
    """
    The extinction (A_lambda) for given wavelength (or vector of wavelengths)
    from the CCM 1989 (and O'Donnell 1994) parameterization. Returns an
    array of extinction values for each wavelength in 'wave'

    INPUTS:
    wave - array of wavelengths (in Angstroms)
    ebv - colour excess E(B-V) float. If a negative ebv is supplied
          fluxes will be reddened rather than dereddened

    OPTIONAL INPUT:
    r_v - float specifying the ratio of total selective
          extinction R(V) = A(V)/E(B-V). If not specified,
          then r_v = 3.1

    OUTPUTS:
    funred - unreddened calibrated flux array, same number of
             elements as wave

    NOTES:
    1. This function was converted from the IDL Astrolib procedure
       last updated in April 1998. All notes from that function
       (provided below) are relevant to this function

    2. (From IDL:) The CCM curve shows good agreement with the Savage & Mathis (1979)
       ultraviolet curve shortward of 1400 A, but is probably
       preferable between 1200 and 1400 A.

    3. (From IDL:) Many sightlines with peculiar ultraviolet interstellar extinction
       can be represented with a CCM curve, if the proper value of
       R(V) is supplied.

    4. (From IDL:) Curve is extrapolated between 912 and 1000 A as suggested by
       Longo et al. (1989, ApJ, 339,474)

    5. (From IDL:) Use the 4 parameter calling sequence if you wish to save the
       original flux vector.

    6. (From IDL:) Valencic et al. (2004, ApJ, 616, 912) revise the ultraviolet CCM
       curve (3.3 -- 8.0 um-1).    But since their revised curve does
       not connect smoothly with longer and shorter wavelengths, it is
       not included here.

    7. For the optical/NIR transformation, the coefficients from
       O'Donnell (1994) are used

    >>> ccm_unred([1000, 2000, 3000], [1, 1, 1], 2 )
    array([9.7976e+012, 1.12064e+07, 32287.1])
    """
    import numpy as np
    scalar = not np.iterable(wave)
    if scalar:
        wave = np.array([wave], float)
    else:
        wave = np.array(wave, float)
    x = 10000.0 / wave
    npts = wave.size
    a = np.zeros(npts, float)
    b = np.zeros(npts, float)
    good = np.where((x > 0.3) & (x < 1.1))
    a[good] = 0.574 * x[good] ** 1.61
    b[good] = -0.527 * x[good] ** 1.61
    good = np.where((x >= 1.1) & (x < 3.3))
    y = x[good] - 1.82
    c1 = np.array([1.0, 0.104, -0.609, 0.701, 1.137, -1.718, -0.827, 1.647,
        -0.505])
    c2 = np.array([0.0, 1.952, 2.908, -3.989, -7.985, 11.102, 5.491, -
        10.805, 3.347])
    a[good] = np.polyval(c1[::-1], y)
    b[good] = np.polyval(c2[::-1], y)
    good = np.where((x >= 3.3) & (x < 8))
    y = x[good]
    F_a = np.zeros(np.size(good), float)
    F_b = np.zeros(np.size(good), float)
    good1 = np.where(y > 5.9)
    if np.size(good1) > 0:
        y1 = y[good1] - 5.9
        F_a[good1] = -0.04473 * y1 ** 2 - 0.009779 * y1 ** 3
        F_b[good1] = 0.213 * y1 ** 2 + 0.1207 * y1 ** 3
    a[good] = 1.752 - 0.316 * y - 0.104 / ((y - 4.67) ** 2 + 0.341) + F_a
    b[good] = -3.09 + 1.825 * y + 1.206 / ((y - 4.62) ** 2 + 0.263) + F_b
    good = np.where((x >= 8) & (x <= 11))
    y = x[good] - 8.0
    c1 = [-1.073, -0.628, 0.137, -0.07]
    c2 = [13.67, 4.257, -0.42, 0.374]
    a[good] = np.polyval(c1[::-1], y)
    b[good] = np.polyval(c2[::-1], y)
    a_v = r_v * ebv
    a_lambda = a_v * (a + b / r_v)
    if scalar:
        a_lambda = a_lambda[0]
    return a_lambda
