from scipy.interpolate import interp1d

def shift_spectrum(wavelengths, intensities, vel):
    c = 3e8
    shiftedWvl = wavelengths * (1.0 + vel / c)
    shiftedIntens = intensities / (1.0 + vel / c)
    interpolatedShift = interp1d(shiftedWvl, shiftedIntens, bounds_error=False, fill_value=(intensities[0], intensities[-1]))(wavelengths)
    return interpolatedShift


def cost_fn(v, wavelengths, intensities, template):
    return np.sum((template - shift_spectrum(wavelengths, intensities, v))**2)