SHARE
TWEET

Untitled

a guest Mar 19th, 2019 49 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # Linearly propagates the error of a velocity given in celestial coordinates to cartesian via the
  2. # formulas identical to the ones in function cartesian_velocity_from_celestial.
  3. def celestial_coords_to_cartesian_error(
  4.             ra, dec, r,
  5.             ra_error, dec_error, r_error,
  6.             pmra, pmdec, rv,
  7.             pmra_error, pmdec_error, rv_error):
  8.     dvx_dra = -cos(dec)*sin(ra)*rv + r*cos(dec)*cos(ra)*pmra - r*sin(dec)*sin(ra)*pmdec
  9.     dvx_ddec = -sin(dec)*cos(ra)*rv - r*sin(dec)*sin(ra)*pmra + r*cos(dec)*cos(ra)*pmdec
  10.     dvx_dr = cos(dec)*sin(ra)*pmra + sin(dec)*cos(ra)*pmdec
  11.     dvx_dpmra = r*cos(dec)*sin(ra)
  12.     dvx_dpmdec = r*sin(dec)*cos(ra)
  13.     dvx_drv = cos(dec)*cos(ra)
  14.  
  15.     dvy_dra = cos(dec)*cos(ra)*rv + r*cos(dec)*sin(ra)*pmra + r*sin(dec)*cos(ra)*pmdec
  16.     dvy_ddec = -sin(dec)*sin(ra)*rv + r*sin(dec)*cos(ra)*pmra + r*cos(dec)*sin(ra)*pmdec
  17.     dvy_dr = -cos(dec)*cos(ra)*pmra + sin(dec)*sin(ra)*pmdec
  18.     dvy_dpmra = -r*cos(dec)*cos(ra)
  19.     dvy_dpmdec = r*sin(dec)*sin(ra)
  20.     dvy_drv = cos(dec)*sin(ra)
  21.  
  22.     dvz_dra = 0
  23.     dvz_ddec = cos(dec)*rv + r*sin(dec)*pmdec
  24.     dvz_dr = -cos(dec)*pmdec
  25.     dvz_dpmra = 0
  26.     dvz_dpmdec = -r*cos(dec)
  27.     dvz_drv = sin(dec)
  28.  
  29.     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),
  30.             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),
  31.             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)]
  32.  
  33. # Linearly propagates the error of the expression |v2 - v1| where v2 and v1 are velocity vectors and
  34. # the errors are given in celestial coordinates. It first propagates the error
  35. # from celestial to cartesian velocity and then from cartesian velocity to the |v2 -v1| expression.
  36. def celestial_magnitude_of_velocity_difference_error(
  37.             ra1, dec1, r1,
  38.             ra1_error, dec1_error, r1_error,
  39.             pmra1, pmdec1, rv1,
  40.             pmra1_error, pmdec1_error, rv1_error,
  41.             ra2, dec2, r2,
  42.             ra2_error, dec2_error, r2_error,
  43.             pmra2, pmdec2, rv2,
  44.             pmra2_error, pmdec2_error, rv2_error):
  45.     v1_error = celestial_coords_to_cartesian_error(
  46.         ra1, dec1, r1,
  47.         ra1_error, dec1_error, r1_error,
  48.         pmra1, pmdec1, rv1,
  49.         pmra1_error, pmdec1_error, rv1_error)
  50.  
  51.     v2_error = celestial_coords_to_cartesian_error(
  52.         ra2, dec2, r2,
  53.         ra2_error, dec2_error, r2_error,
  54.         pmra2, pmdec2, rv2,
  55.         pmra2_error, pmdec2_error, rv2_error)
  56.  
  57.     v1 = cartesian_velocity_from_celestial(ra1, dec1, r1, pmra1, pmdec1, rv1)
  58.     v2 = cartesian_velocity_from_celestial(ra2, dec2, r2, pmra2, pmdec2, rv2)
  59.    
  60.     s = mag(sub(v2, v1))
  61.  
  62.     # note that abs(ds_vx1) == abs(ds_vx2), so we can use same expression
  63.     # for both since squared in error propagation also, 1/s is already
  64.     # taken outside (see return statement)
  65.     ds_dvx = (v2[0] - v1[0])
  66.     ds_dvy = (v2[1] - v1[1])
  67.     ds_dvz = (v2[2] - v1[2])
  68.  
  69.     return (1/s)*sqrt((ds_dvx*v1_error[0])**2 + (ds_dvy*v1_error[1])**2 + (ds_dvz*v1_error[2])**2 +
  70.                       (ds_dvx*v2_error[0])**2 + (ds_dvy*v2_error[1])**2 + (ds_dvz*v2_error[2])**2)
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top