Coordinates

Transformation.py

%matplotlib inline
import math
import matplotlib
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import astropy.time as time
import astropy.coordinates as coordinates
from astropy.coordinates.angles import Angle
import astropy.units as u

def planetRV2ICRS(body, epoch, r, v):

    j2000 = time.Time('J2000')
    T = (epoch - j2000).to('day').value/36525
    d = (epoch - j2000).to('day').value

    Ja = 99.360714 + 4850.4046 * T
    Jb = 175.895369 + 1191.9605 * T
    Jc = 300.323162 + 262.5475 * T
    Jd = 114.012305 + 6070.2476 * T
    Je = 49.511251 + 64.3000 * T

    N = 357.85 + 52.316 * T


    bodies_orientations = {
        "Mercury":{
            "ra": 281.01 - 0.033 * T,
            "dec": 61.45 - 0.005 * T,
            "W": 329.548 + 6.1385025 * d
            },
        "Venus":{
            "ra": 272.76,
            "dec": 67.16,
            "W": 160.20 - 1.4813688 * d
            },
        "Mars":{
            "ra": 317.68143 - 0.1061 * T,
            "dec": 52.88650 - 0.0609 * T,
            "W": 176.630 + 350.89198226 * d
            },
        "Jupiter":{
            "ra": 268.056595 - 0.006499 * T + 0.000117 * math.sin(Ja) + 0.000938 * math.sin(Jb)
            + 0.001432 * math.sin(Jc) + 0.000030 * math.sin(Jd) + 0.002150 * math.sin(Je),
            "dec": 64.495303 + 0.002413 * T + 0.000050 * math.cos(Ja) + 0.000404 * math.cos(Jb)
            + 0.000617 * math.cos(Jc) - 0.000013 * math.cos(Jd) + 0.000926 * math.cos(Je),
            "W": 284.95 + 870.5366420 * d
            },
        "Saturn":{
            "ra": 40.589 - 0.036 * T,
            "dec": 83.537 - 0.004 * T,
            "W": 38.90 + 810.7939024 * d
            },
        "Uranus":{
            "ra": 257.311,
            "dec": -15.175,
            "W": 203.81 - 501.1600928 * d
            },
        "Neptune":{
            "ra": 299.36 + 0.70 * math.sin(N),
            "dec": 43.46 - 0.51 * math.cos(N),
            "W": 253.18 + 536.3128492 * d - 0.48 * math.sin(N)
            },
        "Pluto":{
            "ra": 312.993,
            "dec": 6.163,
            "W": 237.305 - 56.3625225 * d
            },
        "Earth":{
            "ra":  0.00 - 0.641 * T,
            "dec": 90.00 - 0.557 * T,
            "W": 190.147 + 360.9856235 * d
            }
    }
    
    r_trans1 = transform(r, -(90 - bodies_orientations[body.name]['dec']), 'x', u.deg)
    r_trans2 = transform(r_trans1, -(90 + bodies_orientations[body.name]['ra']), 'z', u.deg)
    
    icrs_frame_pos_coord, icrs_frame_vel_coord = coordinates.get_body_barycentric_posvel(body.name, time=epoch, ephemeris=)
    
    r_f = icrs_frame_pos_coord.xyz + r_trans2

    return r_f

phobos_r_mars =[-3.036722254275202E-06, 6.315312728749861E-05, -8.608556793613191E-07] * u.au

planetRV2ICRS(bodies.Mars, time.Time("2017-08-17 00:00", scale='tdb'), phobos_r_mars, _)