In the question Is the Parker Solar Probe's semimajor axis being so close to that of Mercury's just coincidence? Does it help somehow? I used a Python script to plot the heliocentric orbit of the Parker Solar Probe in a frame rotating with Mercury's mean anomaly.

To do this I manually downloaded state vectors from Horizons as text files and then imported them, decoded the format, and converted to numpy arrays for simple processing and plotting.

Could someone show a specific example of how to get those state vectors into a numpy array using astroquery.jplhorizons? I've just learned in this answer of this package's existence.

This way I can just use np.save() or some other method to save the data for future, faster plotting or analysis.

Specifically, the data should be defined as follows:

start  2018-Aug-13 00:00:00
stop   2025-Aug-31 00:00:00
step   1 hour
Solar System Barycenter
Ecliptic and Mean Equinox of Reference
km and km/sec
csv format
x, y, z, vx, vy, vz (state vectors)
Objects: Sun, Mercury, Venus, Earth Geocenter, Parker Solar Probe

enter image description here

enter image description here

def makecubelimits(axis, centers=None, hw=None):
    lims = ax.get_xlim(), ax.get_ylim(), ax.get_zlim()
    if centers == None:
        centers = [0.5*sum(pair) for pair in lims] 

    if hw == None:
        widths  = [pair[1] - pair[0] for pair in lims]
        hw      = 0.5*max(widths)
        ax.set_xlim(centers[0]-hw, centers[0]+hw)
        ax.set_ylim(centers[1]-hw, centers[1]+hw)
        ax.set_zlim(centers[2]-hw, centers[2]+hw)
        print "hw was None so set to:", hw
            hwx, hwy, hwz = hw
            print "ok hw requested: ", hwx, hwy, hwz

            ax.set_xlim(centers[0]-hwx, centers[0]+hwx)
            ax.set_ylim(centers[1]-hwy, centers[1]+hwy)
            ax.set_zlim(centers[2]-hwz, centers[2]+hwz)
            print "nope hw requested: ", hw
            ax.set_xlim(centers[0]-hw, centers[0]+hw)
            ax.set_ylim(centers[1]-hw, centers[1]+hw)
            ax.set_zlim(centers[2]-hw, centers[2]+hw)

    return centers, hw

class Thing(object):
    def __init__(self, name):
        self.name = name

def rotate_it(vec, theta):
    cth, sth = [f(theta) for f in (np.cos, np.sin)]
    x, y, z = vec
    xr = cth*x - sth*y
    yr = cth*y + sth*x
    return np.vstack((xr, yr, z))

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fnames = ['Parker Solar Probe Sun 1h horizons_results.txt',
          'Parker Solar Probe Mercury 1h horizons_results.txt',
          'Parker Solar Probe Venus 1h horizons_results.txt',
          'Parker Solar Probe Earth 1h horizons_results.txt',
          'Parker Solar Probe Spacecraft 1h horizons_results.txt']

names  = ['Sun', 'Mercury', 'Venus', 'Earth', 'Parker']

halfpi, pi, twopi = [f*np.pi for f in [0.5, 1.0, 2.0]]
degs, rads        = 180./pi, pi/180.
AU                = 149597870.700  # kilometers

things = []
JDs, posns, vels, linez = [], [], [], []
for name, fname in zip(names, fnames):
    with open(fname, 'r') as infile:

        lines = infile.read().splitlines()

    iSOE = [i for i, line in enumerate(lines) if "$$SOE" in line][0]
    iEOE = [i for i, line in enumerate(lines) if "$$EOE" in line][0]

    print iSOE, iEOE, lines[iSOE], lines[iEOE]

    lines = [line.split(',') for line in lines[iSOE+1:iEOE]]
    JD  = np.array([float(line[0]) for line in lines])
    pos = np.array([[float(item) for item in line[2:5]] for line in lines])
    vel = np.array([[float(item) for item in line[5:8]] for line in lines])

    pos, vel = [thing.T for thing in pos, vel]

    thing = Thing(name)
    thing.JD    = JD
    thing.pos   = pos
    thing.vel   = vel
    thing.lines = lines

Sun, Mercury, Venus, Earth, Parker = things

theta_M = np.arctan2(Mercury.pos[1], Mercury.pos[0])
theta_P = np.arctan2(Parker.pos[1],  Parker.pos[0] )

T_mercury = 87.969 * 24. # hours
JDx  = 2459597.5
iJDx = 30096
theta_x0  = theta_M[iJDx]
hours_x   = np.arange(len(theta_M), dtype=float) - iJDx
theta_Mx  = theta_x0 + twopi * hours_x / T_mercury

for thing in things:
    thing.pos_rM = rotate_it(thing.pos, -theta_M)
    thing.vel_rM = rotate_it(thing.vel, -theta_M)

    thing.pos_rMx = rotate_it(thing.pos, -theta_Mx)
    thing.vel_rMx = rotate_it(thing.vel, -theta_Mx)

    thing.pos_rP = rotate_it(thing.pos, -theta_P)
    thing.vel_rP = rotate_it(thing.vel, -theta_P)

JD_today      = 2458360.5 # 2018-Aug-30
i_today       = np.argmax(Earth.JD >= JD_today)
i_08_Nov_2024 = np.argmax(Earth.JD >= 2460622.5)

colors = '-y', '-r', '-g', '-b', '-m'

if True:
    fig = plt.figure()  # [12, 10]
    ax1 = fig.add_subplot(2, 1, 1)
    for thing, co in zip(things, colors):
        x, y, z = thing.pos_rMx
        ax1.plot(x, y, co)
        if True:
                     y[i_today:i_today+1], 'ok', markersize=8)
    ax1.set_xlim(-1.6E+08, 1.6E+08)
    ax1.set_ylim(-1.6E+08, 1.6E+08)

    ax2  = fig.add_subplot(2, 1, 2)
    i_start = i_08_Nov_2024
    for thing, co in zip(things, colors):
        x, y, z = thing.pos_rMx[:, i_start:]
        ax2.plot(x, y, co)
        if False:
                     y[i_today:i_today+1], 'ok', markersize=8)
    ax2.set_xlim(-1.6E+08, 1.6E+08)
    ax2.set_ylim(-1.6E+08, 1.6E+08)

if True:    
    fig = plt.figure(figsize=[10, 8])  # [12, 10]
    ax  = fig.add_subplot(1, 1, 1, projection='3d')
    for thing, co in zip(things, colors):
        x, y, z = thing.pos_rM
        ax.plot(x, y, z, co)
                z[i_today:i_today+1], 'ok')

    c, h = makecubelimits(ax, centers=(0, 0, 0), hw=None)
    print c, h

  • 1
    $\begingroup$ You have an example here of how to combine it with poliastro docs.poliastro.space/en/latest/examples/… $\endgroup$
    – user10716
    Commented Sep 14, 2018 at 4:40
  • $\begingroup$ @astrojuanlu I noticed that I'd forgotten to add the link to the other answer (have edited) and just a few days ago I've reinstalled anaconda, meaning I can get back to your answer and give it another go. In both cases thanks for the reminders! $\endgroup$
    – uhoh
    Commented Sep 14, 2018 at 4:49

1 Answer 1


Not sure exactly what you want as the origin and what corrections for aberrations are wanted and whether you want Mercury itself or its barycenter but there didn't appear to be much/any difference in this case (no Mercurian moons...) but the following code defines a mapping for both and also for Parker Solar Probe (PSP) to the correct HORIZONS id. You can find these HORIZONS SPK ids out by typing the names into the web version of HORIZONS. This example produces a daily position for 2020; change the start, end and ephem_step_size to taste but don't abuse the JPL service by asking for giant amounts of output calculations.

from datetime import datetime
from astropy.time import Time
from astroquery.jplhorizons import Horizons

object_id_map = { 'MercuryBody' : 199,
                  'MercuryBary' : 1,
                  'PSP' : -96
start = datetime(2020,1,1)
end = datetime(2021,1,1)
origin = '@ssb'
obj_name = 'MercuryBody'

# Setup HORIZONS query for the desired body by SPK id
eph = Horizons(id=object_id_map.get(obj_name, 1), id_type='id', epochs={'start' : start.strftime("%Y-%m-%d %H:%M"),
            'stop' : end.strftime("%Y-%m-%d %H:%M"), 'step' : ephem_step_size}, location=origin)

# Query HORIZONS for the state vectors, set the reference plane to
# default 'ecliptic' (ecliptic and mean equinox of reference epoch),
# no aberrations (light travel time and stellar aberration not included),
# include time difference of TDB-UTC in output for possible further use (although
# astropy.Time can do this just as well)
table = eph.vectors(refplane='ecliptic', aberrations='geometric', delta_T=True)

# Add a `datetime` column to the table for easier filtering. Needs to be done by making a
# Time array not a Column due to numpy issue (https://github.com/astropy/astropy/issues/9374) May not apply in Astropy 4.0+
dates = Time([datetime.strptime(d, "A.D. %Y-%b-%d %H:%M:%S.%f") for d in table['datetime_str']])
table.add_column(dates, name='datetime')

['targetname', 'datetime_jd', 'datetime_str', 'delta_T', 'x', 'y', 'z', 'vx', 'vy', 'vz', 'lighttime', 'range', 'range_rate']
x = table['x']

The resulting table contains the X, Y, Z components of the vector (in AU) in the x, y, z columns and the velocities (in AU/day) in vz, vy, vz. AstroPy Columns are built on top of numpy arrays and provide the same functionality with additional functionality to handle Quantitys (Astropy Table documentation reference) as illustrated in the last line where we extract the X component.

  • 1
    $\begingroup$ Fixed the missing comma vs period; update datetime parser (format is quite a bit different between the vectors option and the ephemerides option which I'm more used to using) $\endgroup$ Commented Aug 1, 2020 at 0:59
  • $\begingroup$ fixed one more typo, works great thank you very much! $\endgroup$
    – uhoh
    Commented Aug 1, 2020 at 1:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.