<h1>Recovering the motion of the Sun around the barycenter</h1>
<p>The accepted answer perfectly valid, since it contains figures that closely match the data provided in <a href="https://ssd.jpl.nasa.gov/horizons/app.html#/" rel="nofollow noreferrer">JPL Horizons</a>. In case anyone is interested, in my answer I will try to demonstrate how one can essentially reproduce this plot by <em>actually evolving</em> the positions and velocities of the planets in the solar system, only using:</p>
<ul>
<li>The vector data <span class="math-container">$(\vec{r}_i,\vec{v}_i)$</span> of all planets and the Sun at some initial time <span class="math-container">$t_0$</span>.</li>
<li>no specialized software - purely Python's <code>numpy</code> for all calculations.</li>
</ul>
<h1>Solar Motion 1945-1995</h1>
<p>Before going into the specifics, here are the final results that I get if I use the single-day data of January 1st 1945 provided by Horizons and evolve the system for 50 years:</p>
<p><a href="https://i.sstatic.net/P08dP.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/P08dP.png" alt="Motion of the Sun around the Solar-System Barycenter" /></a></p>
<p>On the upper plot, the scales of the x and y axes are the equal and the units are AU. Each black dot corresponds to January 1st of a year. For comparison, the Solar radius is about 0.00465 AU. The result closely matches the <a href="https://upload.wikimedia.org/wikipedia/commons/thumb/0/0c/Solar_system_barycenter.svg/1024px-Solar_system_barycenter.svg.png" rel="nofollow noreferrer">figure</a> found in Wikipedia.</p>
<p>The lower plot displays the evolution of the fractional energy change <span class="math-container">$\Delta E \;/\; E(t=0)$</span>, as a function of time (in days). Since this remains substantially lower than 1, this is an indication that the evolution of the system has not gone horribly wrong.</p>
<h1>Detailed instructions</h1>
<ol>
<li>Download the position and velocity data from JPL Horizons. If you prefer to use their <a href="https://ssd.jpl.nasa.gov/horizons/app.html#/" rel="nofollow noreferrer">online app</a>, select the following:</li>
</ol>
<ul>
<li><p>Ephemeris Type: Vector Table</p>
</li>
<li><p>Target Body: (Select Sun, Mercury, ..., Neptune iteratively)</p>
</li>
<li><p>Coordinate Center: Solar System Barycenter (SSB) [500@0]</p>
</li>
<li><p>Time Specification: Your prefered initial time, step=1 day</p>
</li>
<li><p>Table Settings: {&quot;Output Quantities&quot;: State Vector, &quot;Output Units&quot;: au and days, &quot;CSV format&quot;: True}</p>
<p>The output with contain <span class="math-container">$[x,y,z,v_x,v_y,v_z]$</span> (in AU and AU/days respectively), which is all we need, along with the masses and the gravitational constant <span class="math-container">$G$</span>.</p>
<p><strong>Edit</strong>: Following the suggestion in the comments, here is a Python script that uses directly the Horizons API to get this data, by simply specifying the date ('YYYY-MM-DD'):</p>
</li>
</ul>
<pre><code>    import requests
    import datetime
    def get_planet_vectors(start_time_str):
        &quot;&quot;&quot;Example usage: get_planet_vectors('2022-12-20')&quot;&quot;&quot;
        API_URL = 'https://ssd.jpl.nasa.gov/api/horizons.api'

        command_codes = ['10', '199', '299', '399', '499', '599', '699', '799', '899']
        options = {
            &quot;format&quot;: 'json',
            &quot;MAKE_EPHEM&quot;: 'YES',
            &quot;COMMAND&quot;: None,
            &quot;EPHEM_TYPE&quot;: 'VECTORS',
            &quot;CENTER&quot;: '500@0',
            &quot;START_TIME&quot;: None,
            &quot;STOP_TIME&quot;: None,
            &quot;STEP_SIZE&quot;: '2d',
            &quot;VEC_TABLE&quot;: '2',
            &quot;REF_SYSTEM&quot;: &quot;ICRF&quot;,
            &quot;REF_PLANE&quot;: &quot;ECLIPTIC&quot;,
            &quot;VEC_CORR&quot;: &quot;NONE&quot;,
            &quot;OUT_UNITS&quot;: 'au-d',
            &quot;VEC_LABELS&quot;: &quot;YES&quot;,
            &quot;VEC_DELTA_T&quot;: &quot;NO&quot;,
            &quot;CSV_FORMAT&quot;: &quot;YES&quot;,
            &quot;OBJ_DATA&quot;: &quot;YES&quot;,
        }

        start_time = datetime.datetime.strptime(start_time_str, '%Y-%m-%d')
        stop_time = start_time + datetime.timedelta(days=1)
        
        options['START_TIME'] = start_time.strftime('%Y-%m-%d')
        options['STOP_TIME'] = stop_time.strftime('%Y-%m-%d')
        planet_vectors = []
        for code in command_codes:
            options['COMMAND'] = code
            response = requests.get(API_URL, params=options)
            data = response.json()['result']
            # get the output csv data. It starts with <span class="math-container">$$SOE and ends with $$</span>EOE
            csv_data = data[data.find('<span class="math-container">$$SOE')+5:data.find('$$</span>EOE')-1]
            # strip any final commas and split the data into a list
            csv_data = csv_data.strip(',').split(',')
            # remove the first 2 elements, which are the time in two different formats
            csv_data = csv_data[2:]
            # convert the strings to floats
            csv_data = [float(x) for x in csv_data]
            planet_vectors.append(csv_data)
        return np.array(planet_vectors)
</code></pre>
<ol start="2">
<li>Next we need a differential equation iterative solver to propagate the phase space. Note that Runge-Kutta-type solvers are not a great fit for this task, since they tend in the long-run to violate energy conservation quite severely. <a href="https://www.ias.edu/sites/default/files/sns/Numerical-integration(2).pdf" rel="nofollow noreferrer">This resource</a> is a contains a great introduction in accurately evolving celestial orbits. Below I have written below a 4th order symplectic integrator and a simple adaptive step-size routine, which appear to perform adequately well. Here they are:</li>
</ol>
<pre><code>import numpy as np
G = 0.000295912208 ## Gravitational Constant in AU^3 per Solar Masses per Days^2
class StarSystem:
    &quot;&quot;&quot;
    Class to simulate a solar system with stars and planets
    &quot;&quot;&quot;
    def __init__(self, phase_space=None, masses=None, n_dim=3):
        self.n_dim = n_dim
        self.step_size = 1E-3 # small initially, but adjusts automatically
        self.masses = masses
        self.phase_space = phase_space
        self.evolve = self.symplectic4

    def get_pairwise_separations(self):
        &quot;&quot;&quot;
        Return the array of the pairwise separations of all objects in the solar system
        in the same star system.
        &quot;&quot;&quot;
        num_objs = len(self.masses)
        positions = self.phase_space[:self.phase_space.size//2].reshape(num_objs, self.n_dim)
        # separations[i,j] is the separation vector pointing from the jth object to the ith object,
        # i.e. separations[i,j] = positions[i] - positions[j]
        separations = positions[:,np.newaxis,:] - positions[np.newaxis,:,:] # shape (num_objs, num_objs, n_dim)
        distances = np.linalg.norm(separations,axis=2) # shape (num_objs, num_objs)
        return separations, distances

    def get_accelerations(self):
        &quot;&quot;&quot;
        Return the array of the accelerations for all objects in the star system
        &quot;&quot;&quot;
        separations, distances = self.get_pairwise_separations()
        
        G_over_r3 = G/distances**3
        acceleration_components = - G_over_r3[:,:,np.newaxis] * separations * self.masses[np.newaxis,:,np.newaxis]
        
        accelerations = np.nansum(acceleration_components,axis=1) 
        accelerations = accelerations.flatten()
        return accelerations

    def symplectic4(self, inplace=True):
        &quot;&quot;&quot;
        Implement the symplectic integrator of 4th order to evolve
        Coefficients found in: 
        http://www.slac.stanford.edu/cgi-wrap/getdoc/slac-pub-5071.pdf
        &quot;&quot;&quot;
        h = self.step_size
        x = (2**(1/3) + 2**(-1/3) - 1)/6
        
        c1 = 0
        c3 = - 4*x - 1
        c2 = c4 = 2*x + 1
        d2 = d3 = -x
        d1 = d4 = x + 1/2

        c = np.array([c1,c2,c3,c4])
        d = np.array([d1,d2,d3,d4])

        original_configurations = self.phase_space.copy()

        positions, velocities = (
            self.phase_space.copy()[:self.phase_space.size//2],
            self.phase_space.copy()[self.phase_space.size//2:]
        )

        for i in range(4):
            accelerations = self.get_accelerations()
            velocities += c[i] * accelerations * h
            positions += d[i] * velocities * h
            new_phase_space = np.concatenate((positions, velocities))
            self.phase_space = new_phase_space

        if not inplace:
            # Restore the original configuration
            self.phase_space = original_configurations

        return new_phase_space

    def adapt_step(self, relError=1E-5, inplace=True):
        &quot;&quot;&quot;
        Adapt the step size to reduce local error in positions to some specified relative error.
        &quot;&quot;&quot;
        h = self.step_size
        ## Evolve with two steps of h
        original_configurations = self.phase_space.copy()
        self.evolve(inplace=True)
        self.evolve(inplace=True)
        positions = self.phase_space.copy()[:self.phase_space.size//2].reshape(len(self.masses), self.n_dim)

        self.phase_space = original_configurations
            
        ## Evolve with one step of 2*h
        self.step_size = 2*h
        self.evolve(inplace=True)
        positions2 = self.phase_space.copy()[:self.phase_space.size//2].reshape(len(self.masses), self.n_dim)
        self.phase_space = original_configurations
        self.step_size = h
        
        ## Calculate the maximum relative distance, defined (in 2d) as the maximum of sqrt((x1-x2)^2 + (y1-y2)^2) / sqrt(x1^2 + y1^2)
        # where (x1,y1) and (x2,y2) are the positions achieved with two steps of h and one step of 2*h
        position_differences = positions - positions2
        norm_positions = np.linalg.norm(positions, axis=1)
        norm_position_differences = np.linalg.norm(position_differences, axis=1)

        maxRelDist = np.max(norm_position_differences / norm_positions)
        
        # Using error estimates. Multiply by 0.85 for extra safety.
        h = 0.85 * h * (relError / maxRelDist )**(1/5)

        if inplace:
            self.step_size = h
        
        return h
</code></pre>
<ol start="3">
<li>Evolve the system until you have reached the desired duration, e.g.</li>
</ol>
<pre><code>n_dim = 3
start_time_str = '1945-01-01'
duration = 50 # duration in years

masses_p = np.array([3.285E23, 4.867E24, 5.972E24, 6.39E23, 1.898E27, 5.683E26, 8.681E25, 1.024E26])
masses_p = masses_p / 1.989E30
masses = np.concatenate(([1], masses_p))

planet_vectors = get_planet_vectors(start_time_str)

# create the initial positions as an array [x1, y1, z1, x2, y2, z2, ...]
initial_positions = planet_vectors[:, :n_dim].flatten()
# create the initial velocities as an array [vx1, vy1, vz1, vx2, vy2, vz2, ...]
initial_velocities = planet_vectors[:, n_dim:].flatten()
phase_space = np.concatenate(
    (initial_positions,
    initial_velocities)
)
solar_system = StarSystem(phase_space, masses, n_dim=n_dim)
t_0 = time = 0
t_end = duration * 365.25
sun_position = []
while time&lt;t_end:
    print(f&quot;\r{time:.2f} days. Sun at {solar_system.phase_space[:3]} AU&quot;, end=&quot;&quot;)
    solar_system.evolve(inplace=True)
    # update your plot data using solar_system.phase_space
    time += solar_system.step_size
    sun_position.append([time, *solar_system.phase_space[:3]])
    solar_system.adapt_step(inplace=True)

data = np.array(sun_position)
time = data[:, 0]
x = data[:, 1]
y = data[:, 2]
import matplotlib.pyplot as plt
plt.plot(x, y)
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
</code></pre>

CC BY-SA 4.0
2022-12-20T13:51:41.240Z
Nice work! You may enjoy my version https://astronomy.stackexchange.com/a/28036/16685 which gets Horizons to do the orbit calculation. I have 4th order Verlet / Leapfrog symplectic integration Python code here:  https://astronomy.stackexchange.com/a/48477/16685 which y9u might find interesting to compare against the algorithm you're using.
PM 2Ring
CC BY-SA 4.0
2022-12-20T14:00:01.870Z
You don't seem to have mass values in your code... JPL has a GM table [here](https://ssd.jpl.nasa.gov/astro_par.html); the corresponding value for the Sun (derived from DE440) is 132712440041.279419 km^3/s^2. Using that data, G in km^3/s^2 per solar mass is 0.000295912208284119497, but it's easier to just use combined GM values in code, as I mention here: https://astronomy.stackexchange.com/a/48616/16685
PM 2Ring
CC BY-SA 4.0
2022-12-20T14:03:42.707Z
BTW, it's easy to fetch that vector data directly using the [Horizons API](https://ssd-api.jpl.nasa.gov/doc/horizons.html). Here's an example query for [Jupiter](https://ssd.jpl.nasa.gov/api/horizons.api?format=text&EPHEM_TYPE=V&OUT_UNITS=AU-D&CSV_FORMAT=YES&VEC_TABLE=2&CENTER=%400&TIME_TYPE=UT&TLIST=%271945-Jan-1%27&COMMAND=5)
PM 2Ring
CC BY-SA 4.0
2022-12-20T21:26:25.237Z
@PM2Ring, your link has the Earth as the center, it needs to be the SSB.  This one is correct.  You can change the "**command='0'**" to 1,2,3....10 for each object barycenter.


[link][1]


  [1]: https://ssd.jpl.nasa.gov/api/horizons.api?COMMAND=%270%27&EPHEM_TYPE=%27VECTORS%27&CENTER=%27500@0%27&START_TIME=%272022-01-01%27&STOP_TIME=%272022-01-02%27&STEP_SIZE=%272d%27&format=text&CAL_FORMAT=%27BOTH%27&TIME_DIGITS=%27FRACSEC%27&EXTRA_PREC=%27YES%27&CSV_FORMAT=%27YES%27&VEC_DELTA_T=%27YES%27&VEC_CORR=%27NONE%27
Greg Miller
CC BY-SA 4.0
2022-12-20T21:32:15.847Z
thanks @PM2Ring, these are useful resources! The value of G is now fixed and I added the masses (but I agree that using directly GM is better). I added a routine that gets the data from horizons via its API, so that now one can easily use these scripts for any initial date (and any duration). Also added a simple plotting snippet to plot the 2D projection of the data. Your 3D plot of the Horizons data shows that the 2D projection is enough (there is very little motion perpendicularly to the ecliptic).
Ioannis Michaloliakos
CC BY-SA 4.0
2022-12-20T21:44:47.130Z
@GregMiller No, my query uses the SSB as the center. Your query uses the SSB as both the center and the command, which results in null vectors.
PM 2Ring
CC BY-SA 4.0
2022-12-20T21:51:00.300Z
@Ioannis No worries. :) It's probably best to use the planet barycentres (1, 2, 3, etc) for this application, rather than the planets themselves (199, 299, 399, etc). That also requires adding the GM of the Earth & Moon. The other GMs in that table are already correct.
PM 2Ring
CC BY-SA 4.0
2022-12-20T22:41:48.203Z
Ok! I was able to run it in Sage (in Python mode). The plot looks good, but it reports a divide by zero error in `G_over_r3 = G/distances**3`. I guess that happens when you find a planet's distance from itself.
PM 2Ring
CC BY-SA 4.0
2022-12-20T07:34:29.297Z
2022-12-20T21:22:12.053Z
2022-12-20T21:22:12.053Z
Ioannis Michaloliakos
<p>I quickly found two pictures that illustrate the motion of the Sun with respect to the known Solar System barycentre: from: <a href="https://commons.wikimedia.org/wiki/File:Solar_system_barycenter.svg" rel="nofollow noreferrer">here</a> and <a href="https://web.archive.org/web/20150524071045/https://qph.is.quoracdn.net/main-qimg-118a1175c48f6ff2afd0183dd0040dd4?convert_to_webp=true" rel="nofollow noreferrer">here</a>.</p>
<p><a href="https://upload.wikimedia.org/wikipedia/commons/thumb/0/0c/Solar_system_barycenter.svg/463px-Solar_system_barycenter.svg.png" rel="nofollow noreferrer"><img src="https://upload.wikimedia.org/wikipedia/commons/thumb/0/0c/Solar_system_barycenter.svg/463px-Solar_system_barycenter.svg.png" alt="Motion of the barycenter of the Solar System relative to the Sun" /></a></p>
<p><a href="https://i.sstatic.net/re77M.gif" rel="nofollow noreferrer"><img src="https://i.sstatic.net/re77M.gif" alt="track of the solar centre" /></a></p>
<p>The first claims to show the track of the solar system barycentre in the heliocentric reference frame. The outer yellow circle marks the photosphere of the Sun. The second plot claims to show the track of the centre of the Sun in the barycentric reference frame. The yellow circle shows the photosphere of the Sun to scale. As you can see, the plots are actually (almost) the same! Given that to go from one frame to the other is just a translation, I suppose they can both be right providing the x and y axes are defined appropriately.</p>
<p>To answer the questions posed: &quot;What does it look like&quot; - it looks like these two pictures. &quot;How large is it?&quot; As you can see, the maximum separation between the barycentre and the solar centre appears to be about 2 solar radii over the timescale covered by these plots, but is as small as a tenth of a solar radius (e.g. in 1950). &quot;How elliptic?&quot; Not at all really, it is a complicated superposition caused by the orbits mainly of Jupiter and Saturn, but all the planets contribute to a greater or lesser extent.</p>
<p>The barycentre is calculated from the instantaneous positions of all the discrete masses in the Solar System, by finding their common centre of mass.
<span class="math-container">$$\sum_i m_i({\bf r}_i -{\bf r}_{\rm bc}) = 0\ .$$</span>
I do not know for sure, but I assume the summation for these pictures includes the Sun and all of the planets, and that everything else is negligible at the scale of the thickness of the line. If we define a coordinate system with <span class="math-container">${\bf r}_{\rm bc}=0$</span>, then the position of the Sun is
<span class="math-container">$${\bf r}_\odot = -\frac{\sum_j m_j {\bf r}_j}{M_\odot}\ , $$</span>
where the sum now is over just the planets.</p>

CC BY-SA 3.0
2016-04-27T14:39:16.517Z
Could you include the links and explain the sources of these images in your answer? It would be a lot more interesting and informative if you could explain how these are generated rather than just saying "here they are!".
FJC
CC BY-SA 3.0
2016-04-27T18:28:23.423Z
You can use JPL [HORIZONS](http://ssd.jpl.nasa.gov/horizons.cgi) output to make your own plot; choose Vectors, Sun, and Solar System Barycenter. If each planet's product of mass and orbital radius is proportional to its contribution to the Sun's offset from the barycenter, I figure Jupiter accounts for 49%, Saturn 27%, Uranus 8%, Neptune 15%.
Mike G
CC BY-SA 3.0
2017-06-21T19:44:44.430Z
(Rob, I know this is old but,) sorry for a newbie question, can you point towards some starting point/some details of _how exactly is an orbit sketched if one knows the instantaneous positions of the an object in orbit_?
299792458
CC BY-SA 4.0
2021-01-17T05:47:12.867Z
On studying a [clearer version](http://www.blackholeformulas.com/files/orb2020.gif) of the second image, I think the "average position" of the Sun is about 0.1 $R_\odot$ to the upper right of the barycenter.
Mike G
CC BY-SA 4.0
2021-01-17T08:44:52.630Z
@MikeG an average over what period of time? I would presume that the average position of the Sun (on long timescales) is at the barycentre. A more meaningful statistic would be the average scalar separation from the barycentre.
ProfRob
CC BY-SA 4.0
2021-01-17T18:04:10.980Z
In DE430, the mean barycentric (x, y) of the Sun in years 1800-2300 is (3.04e-4, 3.38e-4) au. I think it's related to the giant planets' orbital eccentricities.
Mike G
CC BY-SA 4.0
2016-04-27T10:05:17.620Z
2022-12-20T09:56:54.163Z
2022-12-20T09:56:54.163Z
ProfRob
<p>Everything in the Solar System revolves around the "barycenter": the overall center of mass. This barycenter is not in the center of the Sun. Some articles and essays I've read go so far as to suggest that the position of the barycenter does not have a set of fixed coordinates within the System: it fluctuates. </p>

<p>Well. Since everything, including the Sun, revolves around this barycenter, the Sun must have its own orbit around it. <strong>What does it look like? How large is it? How elliptic?</strong></p>

<p>(In my research, I have tried and failed to establish whether the barycenter is within the Sun or outside it. Either way, an orbit is an orbit).</p>

CC BY-SA 3.0
2016-04-26T23:29:48.917Z
@RobJeffries: As did yours, obviously. So, were you able to establish, based on those two links, the size and shape of the Sun's orbit? Just curious. Cause I wasn't.
Ricky
CC BY-SA 3.0
2016-04-27T01:44:02.287Z
@RobJeffries: The picture is just as beautiful as the other ones. The fact that it does not inform is just a minor flaw, I suppose.
Ricky
<the-sun><solar-system>
What does the Sun's orbit within the Solar System look like?
Ricky
# What does the Sun's orbit within the Solar System look like?

(Asked by: username_0 on 2016-04-26)

Everything in the Solar System revolves around the "barycenter": the overall center of mass. This barycenter is not in the center of the Sun. Some articles and essays I've read go so far as to suggest that the position of the barycenter does not have a set of fixed coordinates within the System: it fluctuates.

Well. Since everything, including the Sun, revolves around this barycenter, the Sun must have its own orbit around it. **What does it look like? How large is it? How elliptic?**

(In my research, I have tried and failed to establish whether the barycenter is within the Sun or outside it. Either way, an orbit is an orbit).

**Tags:** the-sun, solar-system

> **Comments:**
> **username_0**: @RobJeffries: As did yours, obviously. So, were you able to establish, based on those two links, the size and shape of the Sun's orbit? Just curious. Cause I wasn't.  
> **username_0**: @RobJeffries: The picture is just as beautiful as the other ones. The fact that it does not inform is just a minor flaw, I suppose.  

---

## 2 answers

### Answer by username_4 on 2016-04-27. Score: 18 | Accepted answer

I quickly found two pictures that illustrate the motion of the Sun with respect to the known Solar System barycentre: from: [here](https://commons.wikimedia.org/wiki/File:Solar_system_barycenter.svg) and [here](https://web.archive.org/web/20150524071045/https://qph.is.quoracdn.net/main-qimg-118a1175c48f6ff2afd0183dd0040dd4?convert_to_webp=true).

[![Motion of the barycenter of the Solar System relative to the Sun](https://upload.wikimedia.org/wikipedia/commons/thumb/0/0c/Solar_system_barycenter.svg/463px-Solar_system_barycenter.svg.png)](https://upload.wikimedia.org/wikipedia/commons/thumb/0/0c/Solar_system_barycenter.svg/463px-Solar_system_barycenter.svg.png)

[![track of the solar centre](https://i.sstatic.net/re77M.gif)](https://i.sstatic.net/re77M.gif)

The first claims to show the track of the solar system barycentre in the heliocentric reference frame. The outer yellow circle marks the photosphere of the Sun. The second plot claims to show the track of the centre of the Sun in the barycentric reference frame. The yellow circle shows the photosphere of the Sun to scale. As you can see, the plots are actually (almost) the same! Given that to go from one frame to the other is just a translation, I suppose they can both be right providing the x and y axes are defined appropriately.

To answer the questions posed: "What does it look like" - it looks like these two pictures. "How large is it?" As you can see, the maximum separation between the barycentre and the solar centre appears to be about 2 solar radii over the timescale covered by these plots, but is as small as a tenth of a solar radius (e.g. in 1950). "How elliptic?" Not at all really, it is a complicated superposition caused by the orbits mainly of Jupiter and Saturn, but all the planets contribute to a greater or lesser extent.

The barycentre is calculated from the instantaneous positions of all the discrete masses in the Solar System, by finding their common centre of mass.
$$\sum\_i m\_i({\bf r}\_i -{\bf r}\_{\rm bc}) = 0\ .$$
I do not know for sure, but I assume the summation for these pictures includes the Sun and all of the planets, and that everything else is negligible at the scale of the thickness of the line. If we define a coordinate system with ${\bf r}\_{\rm bc}=0$, then the position of the Sun is
$${\bf r}\_\odot = -\frac{\sum\_j m\_j {\bf r}\_j}{M\_\odot}\ , $$
where the sum now is over just the planets.

> **Comments:**
> **username_5**: Could you include the links and explain the sources of these images in your answer? It would be a lot more interesting and informative if you could explain how these are generated rather than just saying "here they are!".  
> **username_6**: You can use JPL [HORIZONS](http://ssd.jpl.nasa.gov/horizons.cgi) output to make your own plot; choose Vectors, Sun, and Solar System Barycenter. If each planet's product of mass and orbital radius is proportional to its contribution to the Sun's offset from the barycenter, I figure Jupiter accounts for 49%, Saturn 27%, Uranus 8%, Neptune 15%.  
> **username_7**: (Rob, I know this is old but,) sorry for a newbie question, can you point towards some starting point/some details of _how exactly is an orbit sketched if one knows the instantaneous positions of the an object in orbit_?  
> **username_6**: On studying a [clearer version](http://www.blackholeformulas.com/files/orb2020.gif) of the second image, I think the "average position" of the Sun is about 0.1 $R_\odot$ to the upper right of the barycenter.  
> **username_4**: @username_6 an average over what period of time? I would presume that the average position of the Sun (on long timescales) is at the barycentre. A more meaningful statistic would be the average scalar separation from the barycentre.  
> **username_6**: In DE430, the mean barycentric (x, y) of the Sun in years 1800-2300 is (3.04e-4, 3.38e-4) au. I think it's related to the giant planets' orbital eccentricities.  

---

### Answer by username_1 on 2022-12-20. Score: 6

# Recovering the motion of the Sun around the barycenter

The accepted answer perfectly valid, since it contains figures that closely match the data provided in [JPL Horizons](https://ssd.jpl.nasa.gov/horizons/app.html#/). In case anyone is interested, in my answer I will try to demonstrate how one can essentially reproduce this plot by *actually evolving* the positions and velocities of the planets in the solar system, only using:

* The vector data $(\vec{r}\_i,\vec{v}\_i)$ of all planets and the Sun at some initial time $t\_0$.
* no specialized software - purely Python's `numpy` for all calculations.

# Solar Motion 1945-1995

Before going into the specifics, here are the final results that I get if I use the single-day data of January 1st 1945 provided by Horizons and evolve the system for 50 years:

[![Motion of the Sun around the Solar-System Barycenter](https://i.sstatic.net/P08dP.png)](https://i.sstatic.net/P08dP.png)

On the upper plot, the scales of the x and y axes are the equal and the units are AU. Each black dot corresponds to January 1st of a year. For comparison, the Solar radius is about 0.00465 AU. The result closely matches the [figure](https://upload.wikimedia.org/wikipedia/commons/thumb/0/0c/Solar_system_barycenter.svg/1024px-Solar_system_barycenter.svg.png) found in Wikipedia.

The lower plot displays the evolution of the fractional energy change $\Delta E \;/\; E(t=0)$, as a function of time (in days). Since this remains substantially lower than 1, this is an indication that the evolution of the system has not gone horribly wrong.

# Detailed instructions

1. Download the position and velocity data from JPL Horizons. If you prefer to use their [online app](https://ssd.jpl.nasa.gov/horizons/app.html#/), select the following:

* Ephemeris Type: Vector Table
* Target Body: (Select Sun, Mercury, ..., Neptune iteratively)
* Coordinate Center: Solar System Barycenter (SSB) [500@0]
* Time Specification: Your prefered initial time, step=1 day
* Table Settings: {"Output Quantities": State Vector, "Output Units": au and days, "CSV format": True}

  The output with contain $[x,y,z,v\_x,v\_y,v\_z]$ (in AU and AU/days respectively), which is all we need, along with the masses and the gravitational constant $G$.

  **Edit**: Following the suggestion in the comments, here is a Python script that uses directly the Horizons API to get this data, by simply specifying the date ('YYYY-MM-DD'):

```
    import requests
    import datetime
    def get_planet_vectors(start_time_str):
        """Example usage: get_planet_vectors('2022-12-20')"""
        API_URL = 'https://ssd.jpl.nasa.gov/api/horizons.api'

        command_codes = ['10', '199', '299', '399', '499', '599', '699', '799', '899']
        options = {
            "format": 'json',
            "MAKE_EPHEM": 'YES',
            "COMMAND": None,
            "EPHEM_TYPE": 'VECTORS',
            "CENTER": '500@0',
            "START_TIME": None,
            "STOP_TIME": None,
            "STEP_SIZE": '2d',
            "VEC_TABLE": '2',
            "REF_SYSTEM": "ICRF",
            "REF_PLANE": "ECLIPTIC",
            "VEC_CORR": "NONE",
            "OUT_UNITS": 'au-d',
            "VEC_LABELS": "YES",
            "VEC_DELTA_T": "NO",
            "CSV_FORMAT": "YES",
            "OBJ_DATA": "YES",
        }

        start_time = datetime.datetime.strptime(start_time_str, '%Y-%m-%d')
        stop_time = start_time + datetime.timedelta(days=1)
        
        options['START_TIME'] = start_time.strftime('%Y-%m-%d')
        options['STOP_TIME'] = stop_time.strftime('%Y-%m-%d')
        planet_vectors = []
        for code in command_codes:
            options['COMMAND'] = code
            response = requests.get(API_URL, params=options)
            data = response.json()['result']
            # get the output csv data. It starts with $$SOE and ends with $$EOE
            csv_data = data[data.find('$$SOE')+5:data.find('$$EOE')-1]
            # strip any final commas and split the data into a list
            csv_data = csv_data.strip(',').split(',')
            # remove the first 2 elements, which are the time in two different formats
            csv_data = csv_data[2:]
            # convert the strings to floats
            csv_data = [float(x) for x in csv_data]
            planet_vectors.append(csv_data)
        return np.array(planet_vectors)

```

2. Next we need a differential equation iterative solver to propagate the phase space. Note that Runge-Kutta-type solvers are not a great fit for this task, since they tend in the long-run to violate energy conservation quite severely. [This resource](https://www.ias.edu/sites/default/files/sns/Numerical-integration(2).pdf) is a contains a great introduction in accurately evolving celestial orbits. Below I have written below a 4th order symplectic integrator and a simple adaptive step-size routine, which appear to perform adequately well. Here they are:

```
import numpy as np
G = 0.000295912208 ## Gravitational Constant in AU^3 per Solar Masses per Days^2
class StarSystem:
    """
    Class to simulate a solar system with stars and planets
    """
    def __init__(self, phase_space=None, masses=None, n_dim=3):
        self.n_dim = n_dim
        self.step_size = 1E-3 # small initially, but adjusts automatically
        self.masses = masses
        self.phase_space = phase_space
        self.evolve = self.symplectic4

    def get_pairwise_separations(self):
        """
        Return the array of the pairwise separations of all objects in the solar system
        in the same star system.
        """
        num_objs = len(self.masses)
        positions = self.phase_space[:self.phase_space.size//2].reshape(num_objs, self.n_dim)
        # separations[i,j] is the separation vector pointing from the jth object to the ith object,
        # i.e. separations[i,j] = positions[i] - positions[j]
        separations = positions[:,np.newaxis,:] - positions[np.newaxis,:,:] # shape (num_objs, num_objs, n_dim)
        distances = np.linalg.norm(separations,axis=2) # shape (num_objs, num_objs)
        return separations, distances

    def get_accelerations(self):
        """
        Return the array of the accelerations for all objects in the star system
        """
        separations, distances = self.get_pairwise_separations()
        
        G_over_r3 = G/distances**3
        acceleration_components = - G_over_r3[:,:,np.newaxis] * separations * self.masses[np.newaxis,:,np.newaxis]
        
        accelerations = np.nansum(acceleration_components,axis=1) 
        accelerations = accelerations.flatten()
        return accelerations

    def symplectic4(self, inplace=True):
        """
        Implement the symplectic integrator of 4th order to evolve
        Coefficients found in: 
        http://www.slac.stanford.edu/cgi-wrap/getdoc/slac-pub-5071.pdf
        """
        h = self.step_size
        x = (2**(1/3) + 2**(-1/3) - 1)/6
        
        c1 = 0
        c3 = - 4*x - 1
        c2 = c4 = 2*x + 1
        d2 = d3 = -x
        d1 = d4 = x + 1/2

        c = np.array([c1,c2,c3,c4])
        d = np.array([d1,d2,d3,d4])

        original_configurations = self.phase_space.copy()

        positions, velocities = (
            self.phase_space.copy()[:self.phase_space.size//2],
            self.phase_space.copy()[self.phase_space.size//2:]
        )

        for i in range(4):
            accelerations = self.get_accelerations()
            velocities += c[i] * accelerations * h
            positions += d[i] * velocities * h
            new_phase_space = np.concatenate((positions, velocities))
            self.phase_space = new_phase_space

        if not inplace:
            # Restore the original configuration
            self.phase_space = original_configurations

        return new_phase_space

    def adapt_step(self, relError=1E-5, inplace=True):
        """
        Adapt the step size to reduce local error in positions to some specified relative error.
        """
        h = self.step_size
        ## Evolve with two steps of h
        original_configurations = self.phase_space.copy()
        self.evolve(inplace=True)
        self.evolve(inplace=True)
        positions = self.phase_space.copy()[:self.phase_space.size//2].reshape(len(self.masses), self.n_dim)

        self.phase_space = original_configurations
            
        ## Evolve with one step of 2*h
        self.step_size = 2*h
        self.evolve(inplace=True)
        positions2 = self.phase_space.copy()[:self.phase_space.size//2].reshape(len(self.masses), self.n_dim)
        self.phase_space = original_configurations
        self.step_size = h
        
        ## Calculate the maximum relative distance, defined (in 2d) as the maximum of sqrt((x1-x2)^2 + (y1-y2)^2) / sqrt(x1^2 + y1^2)
        # where (x1,y1) and (x2,y2) are the positions achieved with two steps of h and one step of 2*h
        position_differences = positions - positions2
        norm_positions = np.linalg.norm(positions, axis=1)
        norm_position_differences = np.linalg.norm(position_differences, axis=1)

        maxRelDist = np.max(norm_position_differences / norm_positions)
        
        # Using error estimates. Multiply by 0.85 for extra safety.
        h = 0.85 * h * (relError / maxRelDist )**(1/5)

        if inplace:
            self.step_size = h
        
        return h

```

3. Evolve the system until you have reached the desired duration, e.g.

```
n_dim = 3
start_time_str = '1945-01-01'
duration = 50 # duration in years

masses_p = np.array([3.285E23, 4.867E24, 5.972E24, 6.39E23, 1.898E27, 5.683E26, 8.681E25, 1.024E26])
masses_p = masses_p / 1.989E30
masses = np.concatenate(([1], masses_p))

planet_vectors = get_planet_vectors(start_time_str)

# create the initial positions as an array [x1, y1, z1, x2, y2, z2, ...]
initial_positions = planet_vectors[:, :n_dim].flatten()
# create the initial velocities as an array [vx1, vy1, vz1, vx2, vy2, vz2, ...]
initial_velocities = planet_vectors[:, n_dim:].flatten()
phase_space = np.concatenate(
    (initial_positions,
    initial_velocities)
)
solar_system = StarSystem(phase_space, masses, n_dim=n_dim)
t_0 = time = 0
t_end = duration * 365.25
sun_position = []
while time<t_end:
    print(f"\r{time:.2f} days. Sun at {solar_system.phase_space[:3]} AU", end="")
    solar_system.evolve(inplace=True)
    # update your plot data using solar_system.phase_space
    time += solar_system.step_size
    sun_position.append([time, *solar_system.phase_space[:3]])
    solar_system.adapt_step(inplace=True)

data = np.array(sun_position)
time = data[:, 0]
x = data[:, 1]
y = data[:, 2]
import matplotlib.pyplot as plt
plt.plot(x, y)
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

```

> **Comments:**
> **username_2**: Nice work! You may enjoy my version https://astronomy.stackexchange.com/a/28036/16685 which gets Horizons to do the orbit calculation. I have 4th order Verlet / Leapfrog symplectic integration Python code here:  https://astronomy.stackexchange.com/a/48477/16685 which y9u might find interesting to compare against the algorithm you're using.  
> **username_2**: You don't seem to have mass values in your code... JPL has a GM table [here](https://ssd.jpl.nasa.gov/astro_par.html); the corresponding value for the Sun (derived from DE440) is 132712440041.279419 km^3/s^2. Using that data, G in km^3/s^2 per solar mass is 0.000295912208284119497, but it's easier to just use combined GM values in code, as I mention here: https://astronomy.stackexchange.com/a/48616/16685  
> **username_2**: BTW, it's easy to fetch that vector data directly using the [Horizons API](https://ssd-api.jpl.nasa.gov/doc/horizons.html). Here's an example query for [Jupiter](https://ssd.jpl.nasa.gov/api/horizons.api?format=text&EPHEM_TYPE=V&OUT_UNITS=AU-D&CSV_FORMAT=YES&VEC_TABLE=2&CENTER=%400&TIME_TYPE=UT&TLIST=%271945-Jan-1%27&COMMAND=5)  
> **username_3**: @username_2, your link has the Earth as the center, it needs to be the SSB.  This one is correct.  You can change the "**command='0'**" to 1,2,3....10 for each object barycenter.   [link][1]     [1]: https://ssd.jpl.nasa.gov/api/horizons.api?COMMAND=%270%27&EPHEM_TYPE=%27VECTORS%27&CENTER=%27500@0%27&START_TIME=%272022-01-01%27&STOP_TIME=%272022-01-02%27&STEP_SIZE=%272d%27&format=text&CAL_FORMAT=%27BOTH%27&TIME_DIGITS=%27FRACSEC%27&EXTRA_PREC=%27YES%27&CSV_FORMAT=%27YES%27&VEC_DELTA_T=%27YES%27&VEC_CORR=%27NONE%27  
> **username_1**: thanks @username_2, these are useful resources! The value of G is now fixed and I added the masses (but I agree that using directly GM is better). I added a routine that gets the data from horizons via its API, so that now one can easily use these scripts for any initial date (and any duration). Also added a simple plotting snippet to plot the 2D projection of the data. Your 3D plot of the Horizons data shows that the 2D projection is enough (there is very little motion perpendicularly to the ecliptic).  
> **username_2**: @username_3 No, my query uses the SSB as the center. Your query uses the SSB as both the center and the command, which results in null vectors.  
> **username_2**: @username_1 No worries. :) It's probably best to use the planet barycentres (1, 2, 3, etc) for this application, rather than the planets themselves (199, 299, 399, etc). That also requires adding the GM of the Earth & Moon. The other GMs in that table are already correct.  
> **username_2**: Ok! I was able to run it in Sage (in Python mode). The plot looks good, but it reports a divide by zero error in `G_over_r3 = G/distances**3`. I guess that happens when you find a planet's distance from itself.  

---

