# Linearly propagates the error of a velocity given in celestial coordinates to cartesian via the
# formulas identical to the ones in function cartesian_velocity_from_celestial.
def celestial_coords_to_cartesian_error(
            ra, dec, r,
            ra_error, dec_error, r_error,
            pmra, pmdec, rv,
            pmra_error, pmdec_error, rv_error):
    dvx_dra = -cos(dec)*sin(ra)*rv + r*cos(dec)*cos(ra)*pmra - r*sin(dec)*sin(ra)*pmdec
    dvx_ddec = -sin(dec)*cos(ra)*rv - r*sin(dec)*sin(ra)*pmra + r*cos(dec)*cos(ra)*pmdec
    dvx_dr = cos(dec)*sin(ra)*pmra + sin(dec)*cos(ra)*pmdec
    dvx_dpmra = r*cos(dec)*sin(ra)
    dvx_dpmdec = r*sin(dec)*cos(ra)
    dvx_drv = cos(dec)*cos(ra)

    dvy_dra = cos(dec)*cos(ra)*rv + r*cos(dec)*sin(ra)*pmra + r*sin(dec)*cos(ra)*pmdec
    dvy_ddec = -sin(dec)*sin(ra)*rv + r*sin(dec)*cos(ra)*pmra + r*cos(dec)*sin(ra)*pmdec
    dvy_dr = -cos(dec)*cos(ra)*pmra + sin(dec)*sin(ra)*pmdec
    dvy_dpmra = -r*cos(dec)*cos(ra)
    dvy_dpmdec = r*sin(dec)*sin(ra)
    dvy_drv = cos(dec)*sin(ra)

    dvz_dra = 0
    dvz_ddec = cos(dec)*rv + r*sin(dec)*pmdec
    dvz_dr = -cos(dec)*pmdec
    dvz_dpmra = 0
    dvz_dpmdec = -r*cos(dec)
    dvz_drv = sin(dec)

    return [sqrt((dvx_dra*ra_error)**2 + (dvx_ddec*dec_error)**2 + (dvx_dr*r_error)**2 + (dvx_dpmra*pmra_error)**2 + (dvx_dpmdec*pmdec_error)**2 + (dvx_drv*rv_error)**2),
            sqrt((dvy_dra*ra_error)**2 + (dvy_ddec*dec_error)**2 + (dvy_dr*r_error)**2 + (dvy_dpmra*pmra_error)**2 + (dvy_dpmdec*pmdec_error)**2 + (dvy_drv*rv_error)**2),
            sqrt((dvz_dra*ra_error)**2 + (dvz_ddec*dec_error)**2 + (dvz_dr*r_error)**2 + (dvz_dpmra*pmra_error)**2 + (dvz_dpmdec*pmdec_error)**2 + (dvz_drv*rv_error)**2)]

# Linearly propagates the error of the expression |v2 - v1| where v2 and v1 are velocity vectors and
# the errors are given in celestial coordinates. It first propagates the error
# from celestial to cartesian velocity and then from cartesian velocity to the |v2 -v1| expression.
def celestial_magnitude_of_velocity_difference_error(
            ra1, dec1, r1,
            ra1_error, dec1_error, r1_error,
            pmra1, pmdec1, rv1,
            pmra1_error, pmdec1_error, rv1_error,
            ra2, dec2, r2,
            ra2_error, dec2_error, r2_error,
            pmra2, pmdec2, rv2,
            pmra2_error, pmdec2_error, rv2_error):
    v1_error = celestial_coords_to_cartesian_error(
        ra1, dec1, r1,
        ra1_error, dec1_error, r1_error,
        pmra1, pmdec1, rv1,
        pmra1_error, pmdec1_error, rv1_error)

    v2_error = celestial_coords_to_cartesian_error(
        ra2, dec2, r2,
        ra2_error, dec2_error, r2_error,
        pmra2, pmdec2, rv2,
        pmra2_error, pmdec2_error, rv2_error)

    v1 = cartesian_velocity_from_celestial(ra1, dec1, r1, pmra1, pmdec1, rv1)
    v2 = cartesian_velocity_from_celestial(ra2, dec2, r2, pmra2, pmdec2, rv2)

    s = mag(sub(v2, v1))

    # note that abs(ds_vx1) == abs(ds_vx2), so we can use same expression
    # for both since squared in error propagation also, 1/s is already
    # taken outside (see return statement)
    ds_dvx = (v2[0] - v1[0])
    ds_dvy = (v2[1] - v1[1])
    ds_dvz = (v2[2] - v1[2])

    return (1/s)*sqrt((ds_dvx*v1_error[0])**2 + (ds_dvy*v1_error[1])**2 + (ds_dvz*v1_error[2])**2 +
                      (ds_dvx*v2_error[0])**2 + (ds_dvy*v2_error[1])**2 + (ds_dvz*v2_error[2])**2)