celestial_magnitude_of_velocity_difference_error.py

# propagates error of |v2 - v1| where v2 and v1 are vectors and parameters are in celestial coordinates
def celestial_magnitude_of_velocity_difference_error(
            ra1_deg, dec1_deg, r1,
            ra1_error_deg, dec1_error_deg, r1_error,
            pmra1_deg_per_s, pmdec1_deg_per_s, rv1,
            pmra1_error_deg_per_s, pmdec1_error_deg_per_s, rv1_error,
            ra2_deg, dec2_deg, r2,
            ra2_error_deg, dec2_error_deg, r2_error,
            pmra2_deg_per_s, pmdec2_deg_per_s, rv2,
            pmra2_error_deg_per_s, pmdec2_error_deg_per_s, rv2_error):
    ra1 = ra1_deg * conv.deg_to_rad
    ra1_error = ra1_error_deg * conv.deg_to_rad
    dec1 = dec1_deg * conv.deg_to_rad
    dec1_error = dec1_error_deg * conv.deg_to_rad
    pmra1 = pmra1_deg_per_s * conv.deg_to_rad
    pmra1_error = pmra1_error_deg_per_s * conv.deg_to_rad
    pmdec1 = pmdec1_deg_per_s * conv.deg_to_rad
    pmdec1_error = pmdec1_error_deg_per_s * conv.deg_to_rad
    ra2 = ra2_deg * conv.deg_to_rad
    ra2_error = ra2_error_deg * conv.deg_to_rad
    dec2 = dec2_deg * conv.deg_to_rad
    dec2_error = dec2_error_deg * conv.deg_to_rad
    pmra2 = pmra2_deg_per_s * conv.deg_to_rad
    pmra2_error = pmra2_error_deg_per_s * conv.deg_to_rad
    pmdec2 = pmdec2_deg_per_s * conv.deg_to_rad
    pmdec2_error = pmdec2_error_deg_per_s * conv.deg_to_rad

    # first: error propgate to vx, vy, vz (cartesian)

    dvx1_dra = -cos(dec1)*sin(ra1)*rv1 + r1*cos(dec1)*cos(ra1)*pmra1 - r1*sin(dec1)*sin(ra1)*pmdec1
    dvx1_ddec = -sin(dec1)*cos(ra1)*rv1 - r1*sin(dec1)*sin(ra1)*pmra1 + r1*cos(dec1)*cos(ra1)*pmdec1
    dvx1_dr = cos(dec1)*sin(ra1)*pmra1 + sin(dec1)*cos(ra1)*pmdec1
    dvx1_dpmra = r1*cos(dec1)*sin(ra1)
    dvx1_dpmdec = r1*sin(dec1)*cos(ra1)
    dvx1_drv = cos(dec1)*cos(ra1)

    dvy1_dra = cos(dec1)*cos(ra1)*rv1 + r1*cos(dec1)*sin(ra1)*pmra1 + r1*sin(dec1)*cos(ra1)*pmdec1
    dvy1_ddec = -sin(dec1)*sin(ra1)*rv1 + r1*sin(dec1)*cos(ra1)*pmra1 + r1*cos(dec1)*sin(ra1)*pmdec1
    dvy1_dr = -cos(dec1)*cos(ra1)*pmra1 + sin(dec1)*sin(ra1)*pmdec1
    dvy1_dpmra = -r1*cos(dec1)*cos(ra1)
    dvy1_dpmdec = r1*sin(dec1)*sin(ra1)
    dvy1_drv = cos(dec1)*sin(ra1)

    dvz1_dra = 0
    dvz1_ddec = cos(dec1)*rv1 + r1*sin(dec1)*pmdec1
    dvz1_dr = -cos(dec1)*pmdec1
    dvz1_dpmra = 0
    dvz1_dpmdec = -r1*cos(dec1)
    dvz1_drv = sin(dec1)

    error_vx1 = sqrt((dvx1_dra*ra1_error)**2 + (dvx1_ddec*dec1_error)**2 + (dvx1_dr*r1_error)**2 + (dvx1_dpmra*pmra1_error)**2 + (dvx1_dpmdec*pmdec1_error)**2 + (dvx1_drv*rv1_error)**2)
    error_vy1 = sqrt((dvy1_dra*ra1_error)**2 + (dvy1_ddec*dec1_error)**2 + (dvy1_dr*r1_error)**2 + (dvy1_dpmra*pmra1_error)**2 + (dvy1_dpmdec*pmdec1_error)**2 + (dvy1_drv*rv1_error)**2)
    error_vz1 = sqrt((dvz1_dra*ra1_error)**2 + (dvz1_ddec*dec1_error)**2 + (dvz1_dr*r1_error)**2 + (dvz1_dpmra*pmra1_error)**2 + (dvz1_dpmdec*pmdec1_error)**2 + (dvz1_drv*rv1_error)**2)

    dvx2_dra = -cos(dec2)*sin(ra2)*rv2 + r2*cos(dec2)*cos(ra2)*pmra2 - r2*sin(dec2)*sin(ra2)*pmdec2
    dvx2_ddec = -sin(dec2)*cos(ra2)*rv2 - r2*sin(dec2)*sin(ra2)*pmra2 + r2*cos(dec2)*cos(ra2)*pmdec2
    dvx2_dr = cos(dec2)*sin(ra2)*pmra2 + sin(dec2)*cos(ra2)*pmdec2
    dvx2_dpmra = r2*cos(dec2)*sin(ra2)
    dvx2_dpmdec = r2*sin(dec2)*cos(ra2)
    dvx2_drv = cos(dec2)*cos(ra2)

    dvy2_dra = cos(dec2)*cos(ra2)*rv2 + r2*cos(dec2)*sin(ra2)*pmra2 + r2*sin(dec2)*cos(ra2)*pmdec2
    dvy2_ddec = -sin(dec2)*sin(ra2)*rv2 + r2*sin(dec2)*cos(ra2)*pmra2 + r2*cos(dec2)*sin(ra2)*pmdec2
    dvy2_dr = -cos(dec2)*cos(ra2)*pmra2 + sin(dec2)*sin(ra2)*pmdec2
    dvy2_dpmra = -r2*cos(dec2)*cos(ra2)
    dvy2_dpmdec = r2*sin(dec2)*sin(ra2)
    dvy2_drv = cos(dec2)*sin(ra2)

    dvz2_dra = 0
    dvz2_ddec = cos(dec2)*rv2 + r2*sin(dec2)*pmdec2
    dvz2_dr = -cos(dec2)*pmdec2
    dvz2_dpmra = 0
    dvz2_dpmdec = -r2*cos(dec2)
    dvz2_drv = sin(dec2)
    
    error_vx2 = sqrt((dvx2_dra*ra2_error)**2 + (dvx2_ddec*dec2_error)**2 + (dvx2_dr*r2_error)**2 + (dvx2_dpmra*pmra2_error)**2 + (dvx2_dpmdec*pmdec2_error)**2 + (dvx2_drv*rv2_error)**2)
    error_vy2 = sqrt((dvy2_dra*ra2_error)**2 + (dvy2_ddec*dec2_error)**2 + (dvy2_dr*r2_error)**2 + (dvy2_dpmra*pmra2_error)**2 + (dvy2_dpmdec*pmdec2_error)**2 + (dvy2_drv*rv2_error)**2)
    error_vz2 = sqrt((dvz2_dra*ra2_error)**2 + (dvz2_ddec*dec2_error)**2 + (dvz2_dr*r2_error)**2 + (dvz2_dpmra*pmra2_error)**2 + (dvz2_dpmdec*pmdec2_error)**2 + (dvz2_drv*rv2_error)**2)

    v1 = cartesian_velocity_from_celestial(ra1, dec1, r1, pmra1, pmdec1, rv1)
    v2 = cartesian_velocity_from_celestial(ra2, dec2, r2, pmra2, pmdec2, rv2)
    
    # now: error propagate |v2 - v1| = vdiff using vx, vy, vz errors (v2, v1 are vectors here)

    s = len(sub(v2, v1))

    # note that abs(ds_vx1) == abs(ds_vx2), so we can use both for same 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*error_vx1)**2 + (ds_dvy*error_vy1)**2 + (ds_dvz*error_vz1)**2 + 
                (ds_dvx*error_vx2)**2 + (ds_dvy*error_vy2)**2 + (ds_dvz*error_vz2)**2)