def dphi_dop(t, profile, r0_vec, v_vec, d_hat, use_form=False, form_fun=None, interp_table=None):
    '\n\n    Returns the phase shift due to the Doppler delay for subhalos of mass, mass\n\n    TODO: add use_closest option\n\n    '
    v_mag = np.linalg.norm(v_vec, axis=1)
    r0_v = np.einsum('ij, ij -> i', r0_vec, v_vec)
    t0 = ((- r0_v) / np.square(v_mag))
    b_vec = (r0_vec + (v_vec * t0[:, np.newaxis]))
    b_mag = np.linalg.norm(b_vec, axis=1)
    tau = (b_mag / v_mag)
    b_hat = (b_vec / b_mag[:, np.newaxis])
    v_hat = (v_vec / v_mag[:, np.newaxis])
    b_d = np.dot(b_hat, d_hat)
    v_d = np.dot(v_hat, d_hat)
    x = (np.subtract.outer(t, t0) / tau)
    x0 = ((- t0) / tau)
    prefactor = ((const.yr_to_s * const.GN) / ((const.km_s_to_kpc_yr * const.c_light) * np.square(v_mag)))
    if (interp_table is None):
        bd_term = ((np.sqrt((1 + (x ** 2))) + x) - (np.sqrt((1 + (x0 ** 2))) + x0))
        vd_term = (np.arcsinh(x) - np.arcsinh(x0))
        sig = ((bd_term * b_d) - (vd_term * v_d))
        if ('M' in list(profile)):
            prefactor *= profile['M']
            if use_form:
                t_cl = np.maximum(np.minimum(t0, t[(- 1)]), 0)
                x_cl = ((t_cl - t0) / tau)
                r_cl = ((tau * v_mag) * np.sqrt((1 + (x_cl ** 2))))
                rv = (((((3 * profile['M']) / (4 * np.pi)) * (1 / 200)) * (1 / const.rho_crit)) ** (1 / 3))
                form_func = np.where((r_cl < rv), form((r_cl / rv), profile['c']), 1)
                sig = (form_func * sig)
        elif (form_fun is not None):
            t_cl = np.maximum(np.minimum(t0, t[(- 1)]), 0)
            x_cl = ((t_cl - t0) / tau)
            r_cl = ((tau * v_mag) * np.sqrt((1 + (x_cl ** 2))))
            form_func = form_fun(r_cl, profile['rs'], profile['rhos'])
            sig = (form_func * sig)
        else:
            raise ValueError('rho_s, r_s halo description currently requires custom density profile ("USE_FORMTAB")')
    else:
        y = (b_mag / profile['rs'])
        (bd_term0, vd_term0) = interp_table.bd_vd_terms(x0, y)
        y.shape = (1, (- 1))
        y = np.broadcast_to(y, x.shape)
        (bd_term, vd_term) = interp_table.bd_vd_terms(x, y)
        bd_term -= bd_term0
        vd_term -= vd_term0
        sig = ((profile['rhos'] * (profile['rs'] ** 3)) * ((bd_term * b_d) + (vd_term * v_d)))
    sig = (prefactor * sig)
    return np.sum(sig, axis=(- 1))
Returns the phase shift due to the Doppler delay for subhalos of mass, mass

TODO: add use_closest option
def dphi_dop(t, profile, r0_vec, v_vec, d_hat, use_form=False, form_fun=None, interp_table=None):
    '\n\n    Returns the phase shift due to the Doppler delay for subhalos of mass, mass\n\n    TODO: add use_closest option\n\n    '
    v_mag = np.linalg.norm(v_vec, axis=1)
    r0_v = np.einsum('ij, ij -> i', r0_vec, v_vec)
    t0 = ((- r0_v) / np.square(v_mag))
    b_vec = (r0_vec + (v_vec * t0[:, np.newaxis]))
    b_mag = np.linalg.norm(b_vec, axis=1)
    tau = (b_mag / v_mag)
    b_hat = (b_vec / b_mag[:, np.newaxis])
    v_hat = (v_vec / v_mag[:, np.newaxis])
    b_d = np.dot(b_hat, d_hat)
    v_d = np.dot(v_hat, d_hat)
    x = (np.subtract.outer(t, t0) / tau)
    x0 = ((- t0) / tau)
    prefactor = ((const.yr_to_s * const.GN) / ((const.km_s_to_kpc_yr * const.c_light) * np.square(v_mag)))
    if (interp_table is None):
        bd_term = ((np.sqrt((1 + (x ** 2))) + x) - (np.sqrt((1 + (x0 ** 2))) + x0))
        vd_term = (np.arcsinh(x) - np.arcsinh(x0))
        sig = ((bd_term * b_d) - (vd_term * v_d))
        if ('M' in list(profile)):
            prefactor *= profile['M']
            if use_form:
                t_cl = np.maximum(np.minimum(t0, t[(- 1)]), 0)
                x_cl = ((t_cl - t0) / tau)
                r_cl = ((tau * v_mag) * np.sqrt((1 + (x_cl ** 2))))
                rv = (((((3 * profile['M']) / (4 * np.pi)) * (1 / 200)) * (1 / const.rho_crit)) ** (1 / 3))
                form_func = np.where((r_cl < rv), form((r_cl / rv), profile['c']), 1)
                sig = (form_func * sig)
        elif (form_fun is not None):
            t_cl = np.maximum(np.minimum(t0, t[(- 1)]), 0)
            x_cl = ((t_cl - t0) / tau)
            r_cl = ((tau * v_mag) * np.sqrt((1 + (x_cl ** 2))))
            form_func = form_fun(r_cl, profile['rs'], profile['rhos'])
            sig = (form_func * sig)
        else:
            raise ValueError('rho_s, r_s halo description currently requires custom density profile ("USE_FORMTAB")')
    else:
        y = (b_mag / profile['rs'])
        (bd_term0, vd_term0) = interp_table.bd_vd_terms(x0, y)
        y.shape = (1, (- 1))
        y = np.broadcast_to(y, x.shape)
        (bd_term, vd_term) = interp_table.bd_vd_terms(x, y)
        bd_term -= bd_term0
        vd_term -= vd_term0
        sig = ((profile['rhos'] * (profile['rs'] ** 3)) * ((bd_term * b_d) + (vd_term * v_d)))
    sig = (prefactor * sig)
    return np.sum(sig, axis=(- 1))