I remember when I was a kid drawing with my Amstrad the orbits of the satellites of Jupiter and the planets. I was 12 when I started doing that. Later at 16 or so, with a PC, I did my first n-body simulation of the planetary orbits in QBasic based on the Newton's law of gravitation, improved later with some ideas from the codes by S. Aarseth to reduce the time step of the simulation when two particles are very close. After that I did not think anymore about these simple n-body simulations.

At that time the JPL DE200 integration was already available, and S. Moshier had just published a C code reproducing those results (de118i program). Understanding that program is hard, I sometimes think how interesting could be to port it to Java. But the question I asked myself this time is how necessary is to consider some of the effects implemented in that program, like the figure of the Earth and Moon, the tides, or non-linear relativistic corrections to the gravity, to reproduce the JPL DE4xx rectangular coordinates. In other words, how accurate can be a numerical integration of the planets when using only the Newton's law ? So my idea with this educational post/little summer project is to present a simple example program (included in the examples provided with JPARSEC) to propagate orbits and also to show/animate them in Java to obtain an updated implementation of my old programs.

Planetary ephemerides can be computed in a variety of ways, from numerical integration to analytical methods, or by using orbital elements. For an example of the two last methods you can take a look at my post about computing the position of the Sun, Moon, and planets. Here I will show the most simple case of numerical integration in which the target is to show the orbits visually (as in my old programs) but to also compute ephemerides up to the accuracy of the Newton's law of gravitation. As you will see this does not require necessarily an implementation with extended precision arithmetics. I will show some well-known but interesting effects like the movement of the Sun around the barycenter, or how to take advantage of this integrator to create a custom ephemerides file to obtain approximate results without messing with many GB of data from JPL files.

Although it is very simple, as an introduction I will revise how to implement the Newton's law programatically, with some notes about my particular implementation that uses as initial conditions the values obtained from a given JPL integration.

The Newton's law of Gravitation is a simple formula stating the gravity force a body of mass m suffers is proportional to the product of the masses of that body and the other atracting it (M), and inversely proportional to the square of the distance (d):

F (m) = G M m / d² = m a = m dv/dt

The value G is the Gravitational constant, around 6.674 10-¹¹ m³/(kg s²). When there are no other forces acting, the total force is due to gravity and this total force is also the product of the mass of the body (m) by the acceleration that the body experiments, which is defined as the change of the velocity with time (dv/dt).

Eliminating the common m factor we get dv/dt = G M / d², which means that the variation in the velocity of the object of mass m does not depend on its mass, but on the mass of the other body atracting it.

Another point to consider is the vectorial nature of the gravity force, acting in the line from one body to the other. Since the position of the bodies in the 3d space are represented respect axes not aligned with the bodies it is necessary to decompose the force and compute it in each of the axes. When decomposing it in the three ortogonal components or axes (x, y, z), the formula becomes:

dv_u/dt = -G M ū / d³. Notation is bad due to the limitations of this blog, but the idea is that to decompose in components we multiply by a vector ū which is the distance from m to the body M in each axis (dx, dy, dz), and divide by the norm of that vector (d). This will project the calculation of the force (dv/dt) in each axis x, y, z (dv_u/dt). I included the minus sign because it is an atractive force.

From the point of view of the numerical implementation that is the final formula to use. We start with a set of 'particles' with initial conditions for positions (px, py, pz) and velocities (vx, vy, vz), and compute the velocity change in a given time step (dt). Then we add that velocity change to the velocity in each of the three axes and update the position by moving them, for instance for x axis would be vx(after dt) = vx(before) + dt * dv_x/dt, and then the position with px(after dt) = px(before) + vx(after) * dt.

My intention here is to reproduce the DE4xx values from initial conditions obtained from the same integration, after a given integration time. The integration is done by repeating the previous computation in little steps dt, little enough to approximate accurately the effects of gravity. One problem is that when dt is very little the rounding errors of the computations can be important, so there is a lower limit for dt we should not cross unless when using extended precision arithmetics.

Another point is that each JPL integration can potentially use slightly different values for the masses of the planets and the Sun. The header of each JPL integration includes a set of constants called GMX (GMS for the Sun, GMB for the Earth-Moon barycenter, and GMx for the other planets). To reproduce the values as close as possible one should also read these values and not use fixed constants for the masses. The bodies are the Sun, all the planets, and Pluto, but instead of the Earth the integration includes the Earth-Moon barycenter and a separated set of polynomials to compute the geocentric position of the Moon. In this program I will give the option to use the Earth-Moon barycenter or the Earth and the Moon as separate bodies, but due to the complex interactions between the Earth and the Moon the second case cannot be extremely accurate.

As promised, here is the program listing, not specially long after all.

NBodyIntegrator.java

import java.awt.Color;
import java.awt.Graphics2D;
import java.awt.event.KeyEvent;
import java.awt.event.KeyListener;
import java.awt.geom.AffineTransform;
 
import jparsec.ephem.Ephem;
import jparsec.ephem.EphemerisElement;
import jparsec.ephem.EphemerisElement.ALGORITHM;
import jparsec.ephem.Functions;
import jparsec.ephem.Target.TARGET;
import jparsec.ephem.planets.JPLEphemeris;
import jparsec.ephem.planets.OrbitEphem;
import jparsec.ephem.planets.OrbitalElement;
import jparsec.ephem.probes.Spacecraft;
import jparsec.graph.chartRendering.AWTGraphics;
import jparsec.io.image.Picture;
import jparsec.math.Constant;
import jparsec.math.DoubleVector;
import jparsec.math.FastMath;
import jparsec.observer.LocationElement;
import jparsec.time.AstroDate;
import jparsec.util.JPARSECException;
 
/**
 * Simple n-body integrator program showing the Solar System using initial data from JPL ephemerides.
 * It allows to compare the orbits integrated after some time with the results from JPL integration.
 * It shows only the planets and optionally few asteroids or the Earth/Moon instead of the Earth-Moon 
 * barycenter. Use keys ESCAPE, Q, H, C, L, T, I, R, B, F, S, P, +, -, numbers, and cursor to control 
 * options. H will show a quick help. View direction can be changed to follow a given planet. When 
 * fixed on the Solar System Barycenter you can increase the scale and speed to see how the Sun moves 
 * around it, an effect used to detect exoplanets. Space probes can also be shown optionally, with the
 * entire history of space exploration. Simulation speed can be increased with F, but in that case 
 * accuracy will be reduced and after few decades Mercury or other inner planets will appear at wrong 
 * locations.
 * 
 * @author T. Alonso Albi - OAN (Spain)
 */
public class NBodyIntegrator {
 
	// Static content (initial variables to set or used in different methods)
	private static ALGORITHM jplVersion = ALGORITHM.JPL_DE405;
	private static AstroDate startTime = new AstroDate(2000, 1, 1, 0, 0, 0);
	//private static AstroDate startTime = new AstroDate(1969, 6, 28, 0, 0, 0);
	private static int w = 800, h = 800, r = 2, r2 = 2*r+1, fontSize = 20;
	private static double scale = w * 0.5 / 40, scaleFactor = 1.2, angFactor = 10 * Constant.DEG_TO_RAD,
			speed = 9.0, // in seconds per frame, use 20000 for a smooth animation, 10 for good accuracy
			lon = 0, lat = Constant.PI_OVER_TWO; // eye lon/lat in radians
	private static long wait = 3; // ms to wait after drawing, to reduce visual artifacts
	private static Color background = Color.BLACK, foreground = Color.WHITE;
	private static boolean showTrajectory = false, showLabels = true, loopForever = true, eclipticMode = false, 
			accurateInsteadOfFast = true, pause = false, showEarthMoon = false, showBodySize = true, 
			neglectWeakInteractions = false, showProbes = false;
	private static TARGET center = null;
 
	// More initial static content you do not need to modify
	private static double G, angularMomentum, jd, pos[][];
	private static Graphics2D g;
	private static final TARGET targetsWithoutEarthMoon[] = new TARGET[] {TARGET.SUN, TARGET.MERCURY, TARGET.VENUS,
			TARGET.Earth_Moon_Barycenter, TARGET.MARS, TARGET.JUPITER, TARGET.SATURN, TARGET.URANUS,
			TARGET.NEPTUNE};
	private static final TARGET targetsWithEarthMoon[] = new TARGET[] {TARGET.SUN, TARGET.MERCURY, TARGET.VENUS,
			TARGET.EARTH, TARGET.Moon, TARGET.MARS, TARGET.JUPITER, TARGET.SATURN, TARGET.URANUS,
			TARGET.NEPTUNE};
	private static TARGET targets[];
	private static final Color colorsWithoutEarthMoon[] = new Color[] {Color.YELLOW, Color.GRAY, Color.WHITE, 
			Color.CYAN, Color.RED, Color.MAGENTA, Color.ORANGE, Color.GREEN, Color.BLUE};
	private static final Color colorsWithEarthMoon[] = new Color[] {Color.YELLOW, Color.GRAY, Color.WHITE, 
			Color.CYAN, Color.LIGHT_GRAY, Color.RED, Color.MAGENTA, Color.ORANGE, Color.GREEN, Color.BLUE};
	private static Color colors[];
	private static EphemerisElement eph = new EphemerisElement(TARGET.SUN, EphemerisElement.COORDINATES_TYPE.GEOMETRIC,
			EphemerisElement.EQUINOX_J2000, EphemerisElement.GEOCENTRIC, EphemerisElement.REDUCTION_METHOD.IAU_2009,
			EphemerisElement.FRAME.DYNAMICAL_EQUINOX_J2000, jplVersion);
	private static JPLEphemeris jpl;
	private static double jd0;
 
	/**
	 * Main program.
	 * @param args Not used.
	 */
	public static void main(String args[]) {
		try {
			GUI();
		} catch (Exception exc) {
			exc.printStackTrace();
		}
	}
 
	/** Main GUI program. */
	public static void GUI() throws Exception {
		// Get initial ecliptic J2000 vector for the bodies
		colors = colorsWithoutEarthMoon;
		targets = targetsWithoutEarthMoon;
		if (showEarthMoon) {
			colors = colorsWithEarthMoon;
			targets = targetsWithEarthMoon;				
		}
		jpl = new JPLEphemeris(jplVersion);
		jd0 = startTime.jd();
		jd = jd0;
		pos = getJPLPositions(jpl, eph, jd);
		FastMath.EXACT_MODE = accurateInsteadOfFast;
 
		// Prepare the chart and key interaction
		final Picture pic = new Picture(w, h);
		g = pic.getImage().createGraphics();
		AWTGraphics.enableAntialiasing(g);
		g.setFont(g.getFont().deriveFont((float) fontSize));
		clearImage(false);
		pic.show(w, h, "Solar System Simulator", false, true, true);
		pic.getJFrame().addKeyListener(new KeyListener() {
			@Override
			public void keyTyped(KeyEvent e) {}
			@Override
			public void keyReleased(KeyEvent e) {}
			@Override
			public void keyPressed(KeyEvent e) {
				int code = e.getKeyCode();
				if (code == KeyEvent.VK_Q || code == KeyEvent.VK_ESCAPE) // quit
					loopForever = false;
				if (code == KeyEvent.VK_H) // help
					reportHelp();
				if (code == KeyEvent.VK_C) { // compare with JPL
					pause = true;
					reportComparison();
					System.out.println("Press P to continue");
					System.out.println();
					return;
				}
				if (code == KeyEvent.VK_L) // toogle labels
					showLabels = !showLabels;
				if (code == KeyEvent.VK_T) // toogle trajectories
					showTrajectory = !showTrajectory;
				if (code == KeyEvent.VK_I) // invert time direction
					speed *= -1;
				if (code == KeyEvent.VK_F) // faster
					speed *= scaleFactor;
				if (code == KeyEvent.VK_S) // slower
					speed /= scaleFactor;
				if (speed > 0) {
					speed = Math.min(speed, 1000000);
					speed = Math.max(speed, 1);
				} else {
					speed = Math.max(speed, -1000000);
					speed = Math.min(speed, -1);						
				}
				if (code == KeyEvent.VK_F || code == KeyEvent.VK_S) return;
				if (code == KeyEvent.VK_P) { // pause
					pause = !pause;
					System.out.println(pause ? "Paused, press P to continue" : "Running ...");
					return;
				}
				if (code == KeyEvent.VK_R) { // reset date
					try {
						boolean oldPause = pause;
						if (!pause) {
							pause = true;
							Thread.sleep(200);
						}
						jd = jd0;
						pos = getJPLPositions(jpl, eph, jd);
						pause = oldPause;
						System.out.println("Date reset to "+jd);
					} catch (Exception exc) {
						exc.printStackTrace();
					}
				}
				if (code == KeyEvent.VK_B) // Body size
					showBodySize = !showBodySize;
				if (code == KeyEvent.VK_PLUS) // scale up
					scale *= scaleFactor;
				if (code == KeyEvent.VK_MINUS) // scale down
					scale /= scaleFactor;
				scale = Math.min(scale, 100000);
				scale = Math.max(scale, 1);
				if (code == KeyEvent.VK_RIGHT) // lon +
					lon = Functions.normalizeRadians(lon + angFactor);
				if (code == KeyEvent.VK_LEFT) // lon -
					lon = Functions.normalizeRadians(lon - angFactor);
				if (code == KeyEvent.VK_UP) // lat +
					lat = Math.min(lat + angFactor, Constant.PI_OVER_TWO);
				if (code == KeyEvent.VK_DOWN) // lat -
					lat = Math.max(lat - angFactor, -Constant.PI_OVER_TWO);
				TARGET center0 = center;
				if (code == KeyEvent.VK_0) {
					try {
						boolean oldPause = pause;
						if (!pause) {
							pause = true;
							Thread.sleep(200);
						}
						center = null;
						g.setTransform(new AffineTransform());
						pause = oldPause;
						System.out.println("Setting center to SSB");
					} catch (Exception exc) {
						exc.printStackTrace();
					}
				}
				if (code == KeyEvent.VK_1 && targets.length > 0) center = targets[0];
				if (code == KeyEvent.VK_2 && targets.length > 1) center = targets[1];
				if (code == KeyEvent.VK_3 && targets.length > 2) center = targets[2];
				if (code == KeyEvent.VK_4 && targets.length > 3) center = targets[3];
				if (code == KeyEvent.VK_5 && targets.length > 4) center = targets[4];
				if (code == KeyEvent.VK_6 && targets.length > 5) center = targets[5];
				if (code == KeyEvent.VK_7 && targets.length > 6) center = targets[6];
				if (code == KeyEvent.VK_8 && targets.length > 7) center = targets[7];
				if (code == KeyEvent.VK_9 && targets.length > 8) center = targets[8];
				if (center != center0 && center != null) System.out.println("Setting center to "+center.getName());
				if (!pause) return;
				if (code == KeyEvent.VK_L || code == KeyEvent.VK_PLUS || code == KeyEvent.VK_MINUS || code == KeyEvent.VK_B ||
						code == KeyEvent.VK_RIGHT || code == KeyEvent.VK_LEFT || code == KeyEvent.VK_UP || 
						code == KeyEvent.VK_DOWN || code == KeyEvent.VK_I || code == KeyEvent.VK_R || center != center0) {
					try {
						drawAll();
						pic.update();
					} catch (Exception exc) {}
				}
			}
		});
 
		// Main loop
		while (loopForever) {
			if (pause) {
				Thread.sleep(200);
				continue;
			}
			try {
				// Draw all
				drawAll();
 
				// Update image and wait a little to reduce artifacts
				pic.update();
				if (wait > 0) Thread.sleep(wait);
 
				// Integrate for next draw
				Newton();
			} catch (Exception exc) {
				exc.printStackTrace();
			}
		}
		// For a clean close
		pic.getJFrame().dispose();
	}
 
	private static void drawAll() throws Exception {
		// Clean background and draw bodies
		if (!showTrajectory) clearImage(false);
		drawBodies(pos, g);
		if (showTrajectory) clearImage(true);
 
		// Draw labels
		AstroDate astro = new AstroDate(jd);
		g.drawString(astro.toString()+" TDB", fontSize, fontSize);
		String degSymbol = "\u00b0";
		String s = "Lon (bottom): "+Functions.formatAngleAsDegrees(lon, 1)+degSymbol+", Lat: "+Functions.formatAngleAsDegrees(lat, 1)+degSymbol;
		g.drawString(s, w-fontSize-g.getFontMetrics().stringWidth(s), fontSize);
		s = "Scale: "+Functions.formatValue(scale, 3)+" px/AU, speed: "+Functions.formatValue(speed, 1)+" s/frame";
		g.drawString(s, w/2-g.getFontMetrics().stringWidth(s)/2, (fontSize*5)/2);
		if (showLabels) g.drawString("Angular momentum (AU^2 Msun / s): "+Functions.formatValue(angularMomentum, 15), fontSize, h-fontSize);
	}
 
	private static void Newton() throws Exception {
		// 'Fit' the position of the Sun by computing the SSB and forcing it to be at 0
		// This helps to keep constant L long term but will degrade the position of the 
		// Sun respect JPL integrations
		/*
		double cx = 0, cy = 0, cz = 0;
		double vx = 0, vy = 0, vz = 0;
		int sunIndex = -1;
		for (int i=0;i<targets.length;i++) {			
			if (targets[i] == TARGET.SUN) {
				sunIndex = i;
				continue;
			}
			cx += pos[i][0] / targets[i].relativeMass;
			cy += pos[i][1] / targets[i].relativeMass;
			cz += pos[i][2] / targets[i].relativeMass;
			vx += pos[i][3] / targets[i].relativeMass;
			vy += pos[i][4] / targets[i].relativeMass;
			vz += pos[i][5] / targets[i].relativeMass;
		}
		pos[sunIndex][0] = -cx * targets[sunIndex].relativeMass;
		pos[sunIndex][1] = -cy * targets[sunIndex].relativeMass;
		pos[sunIndex][2] = -cz * targets[sunIndex].relativeMass;
		pos[sunIndex][3] = -vx * targets[sunIndex].relativeMass;
		pos[sunIndex][4] = -vy * targets[sunIndex].relativeMass;
		pos[sunIndex][5] = -vz * targets[sunIndex].relativeMass;
		*/
 
 
		// Integration for next drawing, credits for Newton and his apple
		double speedG = -speed * G;
		for (int i=0;i<targets.length-1;i++) {
			for (int j=i+1; j<targets.length;j++) {
				double dx = pos[i][0] - pos[j][0];
				double dy = pos[i][1] - pos[j][1];
				double dz = pos[i][2] - pos[j][2];
				double d2 = dx*dx + dy*dy + dz*dz;
				double d = FastMath.sqrt(d2);
 
				if (targets[j].relativeMass != 0) {
					double g = d2 * targets[j].relativeMass;
					if (!neglectWeakInteractions || g < 1000) { // Neglect interactions from far/little bodies
						double G2 = speedG / g;
						pos[i][3] += G2 * dx / d;
						pos[i][4] += G2 * dy / d;
						pos[i][5] += G2 * dz / d;
					}
				}
				if (targets[i].relativeMass != 0) {
					double g = d2 * targets[i].relativeMass;
					if (!neglectWeakInteractions || g < 1000) { // Neglect interactions from far/little bodies
						double G2 = speedG / g;
						pos[j][3] -= G2 * dx / d;
						pos[j][4] -= G2 * dy / d;
						pos[j][5] -= G2 * dz / d;
					}
				}
			} 
		}
 
		// Move the planets and compute L
		angularMomentum = 0;
		for (int i=0;i<targets.length;i++) {	
			pos[i][0] += speed * pos[i][3];
			pos[i][1] += speed * pos[i][4];
			pos[i][2] += speed * pos[i][5];
			if (showLabels && targets[i].relativeMass > 0) {
				DoubleVector vr = new DoubleVector(pos[i][0], pos[i][1], pos[i][2]);
				DoubleVector vv = new DoubleVector(pos[i][3], pos[i][4], pos[i][5]);
				angularMomentum += vr.crossProduct(vv).norm2() / targets[i].relativeMass;
			}
		}
 
		// Make integration effective on date
		jd += speed / Constant.SECONDS_PER_DAY;
	}
 
	private static double[][] getJPLPositions(JPLEphemeris jpl, EphemerisElement eph, double jd) 
			throws JPARSECException {
		double pos[][] = new double[targets.length][];
		for (int i=0; i<targets.length; i++) {
			double p[] = null;
			if (targets[i].isAsteroid()) {
				int index = OrbitEphem.getIndexOfAsteroid(targets[i].getName());
				p = new double[] {1E10, 1E10, 1E10, 0, 0, 0};
				if (index >= 0) {
					OrbitalElement orbit = OrbitEphem.getOrbitalElementsOfAsteroid(index);
					if (orbit.referenceEquinox != Constant.J2000) orbit.changeToEquinox(Constant.J2000);
					p = OrbitEphem.toEclipticPlane(orbit, OrbitEphem.elliptic(orbit, jd));
					if (!eclipticMode) p = Ephem.eclipticToEquatorial(p, orbit.referenceEquinox, eph);
				}
			} else {
				p = jpl.getPositionAndVelocity(jd, targets[i]);
				if (targets[i] == TARGET.Moon) {
					if (targets[i-1] != TARGET.EARTH) throw new JPARSECException("Earth not found before the Moon");
					if (eclipticMode) p = Ephem.equatorialToEcliptic(p, Constant.J2000, eph);
					p[3] /= Constant.SECONDS_PER_DAY;
					p[4] /= Constant.SECONDS_PER_DAY;
					p[5] /= Constant.SECONDS_PER_DAY;
					pos[i] = Functions.sumVectors(p, pos[i-1]);
					continue;
				}
				if (targets[i] == TARGET.Pluto)
					p = MoonEphem.fromPlutoCenterToPlutoBarycenter(p.clone(), jd, EphemerisElement.REDUCTION_METHOD.IAU_2009, true);
 
				if (eclipticMode) p = Ephem.equatorialToEcliptic(p, Constant.J2000, eph);
			}
			pos[i] = p;
			// Pass velocities to AU/s
			pos[i][3] /= Constant.SECONDS_PER_DAY;
			pos[i][4] /= Constant.SECONDS_PER_DAY;
			pos[i][5] /= Constant.SECONDS_PER_DAY;
		}
 
		// Calculate the Gravitational constant from the header constants, and set 
		// consistent planetary masses
        double gms = Double.parseDouble(jpl.getConstant("GMS"));
        G = gms / Math.pow(Constant.SECONDS_PER_DAY, 2); // AU^3/s^2, assuming Msun = 1
        TARGET.setPlanetaryMassFromJPLEphemeris(jpl);
		return pos;
	}
 
	private static void clearImage(boolean onlyBorders) {
		g.setColor(background);
		AWTGraphics.disableAntialiasing(g); // Better performance in fillRect
		g.fillRect(0, 0, w, onlyBorders ? fontSize*3 : h);
		if (onlyBorders) g.fillRect(0, h-fontSize*3, w, fontSize*3);
		AWTGraphics.enableAntialiasing(g);
		g.setColor(foreground);
	}
 
	private static synchronized void drawBodies(double pos[][], Graphics2D g) {
		// Compute center and translate to it
		int cx = 0, cy = 0, ci = -1;
		if (center != null) {
			for (int i=0; i<pos.length; i++) {
				if (center == targets[i]) {
					double p[] = new double[] {pos[i][0], -pos[i][1], pos[i][2]}; // Do not rotate velocities
					if (lon != -Constant.PI_OVER_TWO) p = Functions.rotateZ(p, lon + Constant.PI_OVER_TWO);
					if (lat != 0) p = Functions.rotateX(p, Constant.PI_OVER_TWO - lat);
					cx = -(int) (p[0] * scale + 0.5);
					cy = -(int) (p[1] * scale + 0.5);
					ci = i;
					g.translate(cx, cy);
				}
			}			
		}
 
		// Draw bodies
		for (int i=0; i<pos.length; i++) {
			double p[] = new double[] {pos[i][0], -pos[i][1], pos[i][2]}; // Do not rotate velocities
			if (lon != -Constant.PI_OVER_TWO) p = Functions.rotateZ(p, lon + Constant.PI_OVER_TWO);
			if (lat != 0) p = Functions.rotateX(p, Constant.PI_OVER_TWO - lat);
 
			g.setColor(foreground);
			if (i < colors.length) g.setColor(colors[i]);
			int px = (int) (p[0] * scale + w / 2 + 0.5);
			int py = (int) (p[1] * scale + h / 2 + 0.5);
 
			int pr = (int) (targets[i].equatorialRadius * scale / Constant.AU + 0.5);
			if (pr > w*2) pr = w*2;
			if (pr > r && showBodySize) {
				int pr2 = 2*pr+1;
				g.fillOval(px-pr, py-pr, pr2, pr2);
				if (showLabels && !showTrajectory)
					g.drawString(targets[i].getName(), px + pr + fontSize, py + fontSize / 2);
			} else {
				g.fillOval(px-r, py-r, r2, r2);
				if (showLabels && !showTrajectory)
					g.drawString(targets[i].getName(), px + r + fontSize, py + fontSize / 2);
			}
		}
		g.setColor(foreground);
		g.drawLine(w/2-r*2, h/2, w/2+r*2, h/2);
		g.drawLine(w/2, h/2-r*2, w/2, h/2+r*2);
 
		// Draw a mark on the 0 longitude reference position
		double ref1[] = LocationElement.parseLocationElementFast(new LocationElement(0, 0, 35 * scale));
		if (lon != -Constant.PI_OVER_TWO) ref1 = Functions.rotateZ(ref1, lon + Constant.PI_OVER_TWO);
		ref1 = Functions.rotateX(ref1, Constant.PI_OVER_TWO - lat);
		int x1 = w/2 + (int) ref1[0];
		int y1 = h/2 + (int) ref1[1];
		double ref2[] = LocationElement.parseLocationElementFast(new LocationElement(0, 0, 33 * scale));
		if (lon != -Constant.PI_OVER_TWO) ref2 = Functions.rotateZ(ref2, lon + Constant.PI_OVER_TWO);
		ref2 = Functions.rotateX(ref2, Constant.PI_OVER_TWO - lat);
		int x2 = w/2 + (int) ref2[0];
		int y2 = h/2 + (int) ref2[1];
		g.drawLine(x1, y1, x2, y2);
		//if (showLabels && !showTrajectory) g.drawString("0\u00b0 Lon", x2 + r + fontSize, y2 + fontSize / 2);
 
		if (showProbes) {
			try {
				OrbitalElement orbit[] = Spacecraft.getAllProbeElements();
				for (int i=0; i<orbit.length; i++) {
					if (orbit[i] == null) continue;
					if (!Spacecraft.isTimeApplicable(orbit[i], jd)) continue;
					if (orbit[i].referenceEquinox != Constant.J2000) orbit[i].changeToEquinox(Constant.J2000);
					double[] p = OrbitEphem.toEclipticPlane(orbit[i], OrbitEphem.elliptic(orbit[i], jd));
					if (!eclipticMode) p = Ephem.eclipticToEquatorial(p, orbit[i].referenceEquinox, eph);
 
					p = new double[] {p[0], -p[1], p[2]}; // Invert y, do not rotate velocities
					if (lon != -Constant.PI_OVER_TWO) p = Functions.rotateZ(p, lon + Constant.PI_OVER_TWO);
					if (lat != 0) p = Functions.rotateX(p, Constant.PI_OVER_TWO - lat);
					int px = (int) (p[0] * scale + w / 2 + 0.5);
					int py = (int) (p[1] * scale + h / 2 + 0.5);
					if (px <= 0 || py <= 0) continue;
					g.fillOval(px-r, py-r, r2, r2);
					if (showLabels && !showTrajectory) {
						String s = orbit[i].name;
						if (s.indexOf("-") > 0) s = s.substring(0, s.indexOf("-"));
						g.drawString(s, px + r + fontSize, py + fontSize / 2);
					}
				}
			} catch (Exception exc) {}
		}
 
		// Reset back translation in Graphics object
		if (ci >= 0) g.translate(-cx, -cy);
	}
 
	private static void reportComparison() {
		try {
			Thread.sleep(200); // Ensure the loop is paused and not computing positions
			double fromJ2000 = (jd - Constant.J2000) / Constant.JULIAN_DAYS_PER_CENTURY;
			if (Math.abs(fromJ2000) > 1) {
				System.out.println("- Comparison with JPL not available so far from J2000 -");
				return;
			}
			double jplPos[][] = getJPLPositions(jpl, eph, jd);
			int dec = 12;
			String sep = "  ", mode = eclipticMode ? "ecliptic" : "equatorial";
			System.out.println("- Comparison with JPL DE"+jpl.getJPLVersion()+" ("+mode+" J2000 coordinates) -");
			System.out.println("JD (TDB): "+jd+" ("+Functions.formatValue(jd-jd0, 3)+" days from starting time)");
			double maxDif = -1, maxDifFromEarth = -1;
			TARGET maxT = null, maxTEarth = null;
			int baryc = -1;
			for (int i=0;i<targets.length;i++) {
				if (targets[i] == TARGET.Earth_Moon_Barycenter) {
					baryc = i;
					break;
				}
			}
			for (int i=0;i<targets.length;i++) {
				String ps = "x-y-z (JPL):     "+sep+Functions.formatValue(jplPos[i][0], dec)+sep+
						Functions.formatValue(jplPos[i][1], dec)+sep+
						Functions.formatValue(jplPos[i][2], dec);
				String vs = "vx-vy-vz (JPL):  "+sep+Functions.formatValue(jplPos[i][3], dec)+sep+
						Functions.formatValue(jplPos[i][4], dec)+sep+
						Functions.formatValue(jplPos[i][5], dec);
				System.out.println(targets[i].getName());
				System.out.println(ps);
				System.out.println(vs);
				ps = "x-y-z (NBody):   "+sep+Functions.formatValue(pos[i][0], dec)+sep+
						Functions.formatValue(pos[i][1], dec)+sep+
						Functions.formatValue(pos[i][2], dec);
				vs = "vx-vy-vz (NBody):"+sep+Functions.formatValue(pos[i][3], dec)+sep+
						Functions.formatValue(pos[i][4], dec)+sep+
						Functions.formatValue(pos[i][5], dec);
				System.out.println(ps);
				System.out.println(vs);
				if (targets[i].isAsteroid()) continue;
 
				double dif = Math.max(Math.abs(jplPos[i][0]-pos[i][0]), Math.abs(jplPos[i][1]-pos[i][1]));
				dif = Math.max(dif, Math.abs(jplPos[i][2]-pos[i][2]));
				if (dif > maxDif || maxDif < 0) {
					maxDif = dif;
					maxT = targets[i];
				}
				if (baryc >= 0 && baryc != i) {
					double JPLdx = jplPos[i][0]-jplPos[baryc][0];
					double JPLdy = jplPos[i][1]-jplPos[baryc][1];
					double JPLdz = jplPos[i][2]-jplPos[baryc][2];
					double JPLdr = Math.sqrt(JPLdx*JPLdx + JPLdy*JPLdy + JPLdz*JPLdz);
 
					double dx = pos[i][0]-pos[baryc][0] - JPLdx;
					double dy = pos[i][1]-pos[baryc][1] - JPLdy;
					double dz = pos[i][2]-pos[baryc][2] - JPLdz;
 
					double r = Math.sqrt(dx*dx + dy*dy + dz*dz);
					double dr = Math.atan(r / JPLdr) * Constant.RAD_TO_ARCSEC;
					if (dr > maxDifFromEarth || maxDifFromEarth < 0) {
						maxDifFromEarth = dr;
						maxTEarth = targets[i];
					}
				}
			}
			System.out.println("Maximum difference (AU): "+maxDif+" (for "+maxT.getName()+")");
			if (maxDifFromEarth >= 0)
				System.out.println("Maximum difference from Earth (\"): "+maxDifFromEarth+" (for "+maxTEarth.getName()+")");
		} catch (Exception e1) {
			e1.printStackTrace();
		}
	}
 
	private static void reportHelp() {
		System.out.println("- LIST OF KEYS -");
		System.out.println("Q/Esc to quit");
		System.out.println("P to pause");
		System.out.println("C to compare with JPL DE"+jpl.getJPLVersion());
		System.out.println("L to toggle labels");
		System.out.println("T to toggle trajectories");
		System.out.println("I to invert time direction");
		System.out.println("R to reset to start time");
		System.out.println("B to toggle drawing bodies with their sizes");
		System.out.println("F/S for faster/slower speed");
		System.out.println("+/- to increase/decrease the scale");
		System.out.println("Cursor keys to change view angle (lon/lat)");
		System.out.println("[0-9] to choose object to follow. 0 to reset to SSB");
		System.out.println("H for this list");
		System.out.println();
	}
}

The implementation is very simple with some static variables at the beginning that can be modified to choose a given option. There are some keys to control the animation, and the time step can be even inverted to go back in time. The default speed value is 9s, which is the minimum value that can be used for maximum accuracy in double precision. You can increase it to really see the planets moving in the screen in real time (or press f key for faster animation), but the accuracy will be lower. In case you use very high values of speed the position of Mercury will be completely wrong even after few decades.

The array pos[][] contains the positions and velocities of all the bodies. In method getJPLPositions the array is initialized with data from a given integration. This method can include some asteroids (considering their masses) like Ceres, Pallas, and Vesta, allowing to propagate their orbits with the planetary perturbations. Then the method Newton computes the gravity force as described above, taking advantage of its simmetric nature to speed up the calculations. It is important to remark here the selection of units I made: AU for distances and AU/s for velocities. The JPL integration provides velocities in AU/day, but here I decided to use seconds as time unit for the Newton method. The mass unit is the solar mass.

Fig. 1: Two screenshots of the program, the second one showing the orbits of the inner planets.

In method drawBodies everything is put on the screen. The only remarkable option here is to show the space probes with orbital elements included in JPARSEC (without propagating perturbations). This allows to show the entire history of space exploration, although it seems I should update it since many probes are being launched again to Mars. When you zoom in the bodies will appear with their real sizes to scale, but this can be disabled to show them as points. This feature combined with the option to show trajectories can be used to represent the movement of the Sun around the Solar System barycenter (static point representing the location of the center of mass of the entire Solar System). Since most of the mass is in the Sun the barycenter is very close to the solar surface and playing a little you can easily identify how the different planets affect with their orbital periods and masses to the little wobble of the Sun. The velocity variations will also produce a wobble effect in the radial velocity and frequency of the spectral lines of the stars, and this has been used to detect the first exoplanets. Many of them has been detected also because of the little brightness drop when a planet passes in front of the star, but this requires a perfect alignment of the orbital plane of the exoplanet with us and many exoplanets cannot be detected with that method (especially far away ones, that may transit once per many years). Now the great positional accuracy of the Gaia satellite is used to direcly detect the wobble and infer the presence of exoplanets around any star, independently of the orientation of the orbital planes of their planets respect us.

Fig. 2: Movement of the Sun around the Solar System barycenter (crux) between 2020 and 2120, as visible from the north ecliptic pole, obtained from the NBodyIntegrator program. The Sun to scale is drawn on the bottom-right corner.

The method reportComparison is triggered with the c key to compare the positions and velocities computed by the program with those computed with the JPL integration. When using between 9s and 20s for the speed you will notice excellent agreement for the outer planets. The accuracy from Earth for the planets close to us will be around 0.02” after one day and 0.1” after ten days. For outer planets about ten times better or more. We will see that with numbers below.

The previous program is started with a main method. Instead of executing the graphical display you can call something similar to the following numerical integration test program:

/** Numerical integration test without a GUI */
public static void numericalIntegrationTest() throws Exception {
	targets = showEarthMoon ? targetsWithEarthMoon : targetsWithoutEarthMoon;
	jpl = new JPLEphemeris(jplVersion);
	jd0 = startTime.jd();
	//jd0 -= 5000; // Modify start date for other quick tests
	speed = 20;
 
	jd = jd0;
	pos = getJPLPositions(jpl, eph, jd);
	FastMath.EXACT_MODE = accurateInsteadOfFast;
 
	double endDate = jd + 180;
	if (endDate < jd) speed = -Math.abs(speed);
	long t0 = System.currentTimeMillis();
	int iter = 0;
	while (Math.abs(endDate-jd) > Math.abs(speed/Constant.SECONDS_PER_DAY)) {
		Newton();
		iter ++;
		jd = jd0 + iter * speed / Constant.SECONDS_PER_DAY;
	}
	speed = (endDate - jd) * Constant.SECONDS_PER_DAY;
	if (speed != 0) Newton();
	long t1 = System.currentTimeMillis();
	if (endDate != jd) System.out.println("*** ERROR *** Target end date "+endDate+
			" different from the integration final date "+jd);
	reportComparison();
	System.out.println("Integration took "+(float) ((t1-t0)*0.001)+" s");
}

The previous method computes the initial conditions and then use the Newton method to integrate with a given speed forward or backward in time until the integration end time. Then a comparison is shown against the JPL integration. The integration with 20s of speed is performed over 180 days (a lot!), and returns the following results in my laptop:

- Comparison with JPL DE405 (equatorial J2000 coordinates) -
JD (TDB): 2451724.5 (180.000 days from starting time)
Sun
x-y-z (JPL):       -0.006032658552  -0.003741542907  -0.001420728930
vx-vy-vz (JPL):    0.000000000080  -0.000000000062  -0.000000000029
x-y-z (NBody):     -0.006032658640  -0.003741542065  -0.001420728592
vx-vy-vz (NBody):  0.000000000080  -0.000000000062  -0.000000000029
Mercury
x-y-z (JPL):       -0.058103868552  -0.414395272621  -0.215375832955
vx-vy-vz (JPL):    0.000000258350  -0.000000008197  -0.000000031160
x-y-z (NBody):     -0.058104802862  -0.414394668634  -0.215375413557
vx-vy-vz (NBody):  0.000000258350  -0.000000008198  -0.000000031161
Venus
x-y-z (JPL):       -0.242561212403  0.609559831691  0.289464459217
vx-vy-vz (JPL):    -0.000000221752  -0.000000076642  -0.000000020442
x-y-z (NBody):     -0.242563460323  0.609561187029  0.289465211195
vx-vy-vz (NBody):  -0.000000221751  -0.000000076641  -0.000000020442
Earth-Moon barycenter
x-y-z (JPL):       0.127769131873  -0.928376971857  -0.402297058027
vx-vy-vz (JPL):    0.000000194230  0.000000023296  0.000000010098
x-y-z (NBody):     0.127775059703  -0.928376235696  -0.402296754735
vx-vy-vz (NBody):  0.000000194230  0.000000023296  0.000000010098
Mars
x-y-z (JPL):       -0.250313769347  1.422640449575  0.659416143757
vx-vy-vz (JPL):    -0.000000153806  -0.000000011682  -0.000000001198
x-y-z (NBody):     -0.250316226093  1.422638328034  0.659415237094
vx-vy-vz (NBody):  -0.000000153806  -0.000000011682  -0.000000001198
Jupiter
x-y-z (JPL):       3.031599995432  3.669672591247  1.499124814774
vx-vy-vz (JPL):    -0.000000070447  0.000000051951  0.000000023984
x-y-z (NBody):     3.031600458188  3.669671954894  1.499124530731
vx-vy-vz (NBody):  -0.000000070447  0.000000051951  0.000000023984
Saturn
x-y-z (JPL):       5.591660709347  6.764515860547  2.553313687645
vx-vy-vz (JPL):    -0.000000054434  0.000000035560  0.000000017030
x-y-z (NBody):     5.591661175088  6.764515473002  2.553313507543
vx-vy-vz (NBody):  -0.000000054434  0.000000035560  0.000000017030
Uranus
x-y-z (JPL):       14.897614837489  -12.060738395909  -5.493007971301
vx-vy-vz (JPL):    0.000000029919  0.000000029372  0.000000012441
x-y-z (NBody):     14.897614541978  -12.060738665190  -5.493008085067
vx-vy-vz (NBody):  0.000000029919  0.000000029372  0.000000012441
Neptune
x-y-z (JPL):       17.265901881373  -22.680489235848  -9.713104446249
vx-vy-vz (JPL):    0.000000029531  0.000000019747  0.000000007347
x-y-z (NBody):     17.265901597032  -22.680489418269  -9.713104513845
vx-vy-vz (NBody):  0.000000029531  0.000000019747  0.000000007347
Maximum difference (AU): 5.927829813873187E-6 (for Earth-Moon barycenter)
Maximum difference from Earth ("): 2.4511604012101915 (for Mercury)
Integration took 0.76 s

This test show accuracies around 10^-9 AU for the Sun (only 150 m!) and 10^-6 AU (150 km) for the outer planets using a relatively high value of the speed (20s), a great integration time of 180 days, and calculations in double precision mode only. Of course the accuracy can be increased by reducing the speed and the integration time. Let's see what happens with 9s and 10 days:

- Comparison with JPL DE405 (equatorial J2000 coordinates) -
JD (TDB): 2451554.5 (10.000 days from starting time)
Sun
x-y-z (JPL):       -0.007084986624  -0.002710890366  -0.000951626978
vx-vy-vz (JPL):    0.000000000063  -0.000000000077  -0.000000000035
x-y-z (NBody):     -0.007084986566  -0.002710890425  -0.000951627005
vx-vy-vz (NBody):  0.000000000063  -0.000000000077  -0.000000000035
Mercury
x-y-z (JPL):       0.073742795504  -0.396185425220  -0.219513620709
vx-vy-vz (JPL):    0.000000255220  0.000000075210  0.000000013718
x-y-z (NBody):     0.073743216686  -0.396184693214  -0.219513273372
vx-vy-vz (NBody):  0.000000255220  0.000000075210  0.000000013718
Venus
x-y-z (JPL):       -0.692409955824  -0.220934911643  -0.055755857812
vx-vy-vz (JPL):    0.000000070655  -0.000000202548  -0.000000095591
x-y-z (NBody):     -0.692409592344  -0.220935061208  -0.055755948099
vx-vy-vz (NBody):  0.000000070655  -0.000000202548  -0.000000095591
Earth-Moon barycenter
x-y-z (JPL):       -0.344426805301  0.844813611199  0.366494228132
vx-vy-vz (JPL):    -0.000000190223  -0.000000063430  -0.000000027501
x-y-z (NBody):     -0.344426935789  0.844813414644  0.366494143116
vx-vy-vz (NBody):  -0.000000190223  -0.000000063430  -0.000000027501
Mars
x-y-z (JPL):       1.383134535178  0.129701457011  0.022192225611
vx-vy-vz (JPL):    -0.000000008884  0.000000159005  0.000000073173
x-y-z (NBody):     1.383134447791  0.129701588317  0.022192288202
vx-vy-vz (NBody):  -0.000000008884  0.000000159005  0.000000073173
Jupiter
x-y-z (JPL):       3.950257913739  2.789441305836  1.099449361810
vx-vy-vz (JPL):    -0.000000053869  0.000000067260  0.000000030142
x-y-z (NBody):     3.950257864311  2.789441357815  1.099449385294
vx-vy-vz (NBody):  -0.000000053869  0.000000067260  0.000000030142
Saturn
x-y-z (JPL):       6.358438049285  6.205358229132  2.289379644655
vx-vy-vz (JPL):    -0.000000049887  0.000000040498  0.000000018872
x-y-z (NBody):     6.358438006900  6.205358261266  2.289379659751
vx-vy-vz (NBody):  -0.000000049887  0.000000040498  0.000000018872
Uranus
x-y-z (JPL):       14.450191832477  -12.485560799374  -5.672740512421
vx-vy-vz (JPL):    0.000000030999  0.000000028469  0.000000012030
x-y-z (NBody):     14.450191857738  -12.485560775673  -5.672740502397
vx-vy-vz (NBody):  0.000000030999  0.000000028469  0.000000012030
Neptune
x-y-z (JPL):       16.829465510886  -22.966959644683  -9.819493012110
vx-vy-vz (JPL):    0.000000029895  0.000000019260  0.000000007139
x-y-z (NBody):     16.829465535424  -22.966959628684  -9.819493006173
vx-vy-vz (NBody):  0.000000029895  0.000000019260  0.000000007139
Maximum difference (AU): 7.320058659776585E-7 (for Mercury)
Maximum difference from Earth ("): 0.16726068313039802 (for Mercury)
Integration took 0.127 s

As you see the accuracy is around 10x better. These results show that a custom ephemerides method could be created from these results. One just need to write to a file the JPL initial conditions once or twice per year, or more if greater accuracy is needed, and then read the record for the closest initial conditions for a given calculation time and integrate from the record's time to the calculation time.

Results are not significantly different when using other start times. You may be interested in the results when separating the Earth-Moon barycenter into the two bodies. Results are more than decent, but I leave that test to you. Another question that remains is if extended precision arithmetics can improve these results, and if the performace could be good enough for ephemerides calculations.

Java offers BigDecimal for extended precision maths. If you write the NBodyIntegrator using BigDecimal you inmediately notice the extremely poor performance of the code. Too slow to be useful. Fortunately there is a library called Apfloat which is easier to use than BigDecimal, and a lot faster.

I present below the main methods of the previous program implemented using Apfloat. In case you are interested in the complete program you can ask me to post or send it to you.

private static Apfloat speedG;
private static Apfloat G2, d, d2, massRatio, gravityX, gravityY, gravityZ, dx, dy, dz, dx2, dy2, dz2;
private static Apfloat timePlus;
private static Apfloat mass[];
private static void Newton() throws Exception {
	// Prepare some variables
	if (speedG == null) speedG = speed.multiply(G);
	if (timePlus == null) timePlus = speed.divide(new Apfloat(Constant.SECONDS_PER_DAY, precisionDigits));
	if (mass == null) {
		mass = new Apfloat[targets.length];
		for (int i=0;i<targets.length;i++) {
			double rm = targets[i].relativeMass;
			if (rm == 0) rm = 1.0 / Constant.SUN_MASS; // minimum mass of 1 kg to avoid zeros
			mass[i] = new Apfloat(rm, precisionDigits);
		}
	}
 
	// Integration for next drawing, credits for Newton and his apple
	for (int i=0;i<targets.length-1;i++) {
		for (int j=i+1; j<targets.length;j++) {
			if (targets[j].relativeMass == 0) continue;
			dx = pos[i][0].subtract(pos[j][0]);
			dy = pos[i][1].subtract(pos[j][1]);
			dz = pos[i][2].subtract(pos[j][2]);
			dx2 = dx.multiply(dx);
			dy2 = dy.multiply(dy);
			dz2 = dz.multiply(dz);
			d2 = dx2.add(dy2).add(dz2);
			d = ApfloatMath.sqrt(d2);
 
			G2 = speedG.divide(d2.multiply(mass[j]));
 
			gravityX = G2.multiply(dx).divide(d);
			gravityY = G2.multiply(dy).divide(d);
			gravityZ = G2.multiply(dz).divide(d);
			pos[i][3] = pos[i][3].subtract(gravityX);
			pos[i][4] = pos[i][4].subtract(gravityY);
			pos[i][5] = pos[i][5].subtract(gravityZ);
 
			massRatio = mass[j].divide(mass[i]);
 
			pos[j][3] = pos[j][3].add(gravityX.multiply(massRatio));
			pos[j][4] = pos[j][4].add(gravityY.multiply(massRatio));
			pos[j][5] = pos[j][5].add(gravityZ.multiply(massRatio));
		}
	}
 
	// Move the planets and compute L
	if (showLabels) angularMomentum = new Apfloat(0, precisionDigits);
	for (int i=0;i<targets.length;i++) {			
		pos[i][0] = pos[i][0].add(speed.multiply(pos[i][3]));
		pos[i][1] = pos[i][1].add(speed.multiply(pos[i][4]));
		pos[i][2] = pos[i][2].add(speed.multiply(pos[i][5]));
		if (showLabels) {
			DoubleVector vr = new DoubleVector(pos[i][0].doubleValue(), pos[i][1].doubleValue(), pos[i][2].doubleValue());
			DoubleVector vv = new DoubleVector(pos[i][3].doubleValue(), pos[i][4].doubleValue(), pos[i][5].doubleValue());
			angularMomentum = angularMomentum.add(new Apfloat(vr.crossProduct(vv).norm2() / targets[i].relativeMass, precisionDigits));
		}
	}
 
	// Make integration effective on date
	jd = jd.add(timePlus);
}
 
private static Apfloat[][] getJPLPositions(JPLEphemeris jpl, EphemerisElement eph, Apfloat jd) 
		throws JPARSECException {
	Apfloat pos[][] = new Apfloat[targets.length][6];
	for (int i=0; i<targets.length; i++) {
		double p[] = null;
		if (targets[i].isAsteroid()) {
			int index = OrbitEphem.getIndexOfAsteroid(targets[i].getName());
			p = new double[] {1E10, 1E10, 1E10, 0, 0, 0};
			if (index >= 0) {
				OrbitalElement orbit = OrbitEphem.getOrbitalElementsOfAsteroid(index);
				if (orbit.referenceEquinox != Constant.J2000) orbit.changeToEquinox(Constant.J2000);
				p = OrbitEphem.toEclipticPlane(orbit, OrbitEphem.elliptic(orbit, jd.doubleValue()));
				if (!eclipticMode) p = Ephem.eclipticToEquatorial(p, orbit.referenceEquinox, eph);
			}
		} else {
			p = jpl.getPositionAndVelocity(BigDecimal.valueOf(jd.doubleValue()), targets[i]);
			if (targets[i] == TARGET.Moon) {
				if (targets[i-1] != TARGET.EARTH) throw new JPARSECException("Earth not found before the Moon");
				if (eclipticMode) p = Ephem.equatorialToEcliptic(p, Constant.J2000, eph);
				for (int j=0;j<p.length;j++) {
					if (j > 2) p[j] /= Constant.SECONDS_PER_DAY;
					pos[i][j] = new Apfloat(p[j] + pos[i-1][j].doubleValue(), precisionDigits);						
				}
				continue;
			}
			if (targets[i] == TARGET.Pluto)
				p = MoonEphem.fromPlutoCenterToPlutoBarycenter(p.clone(), jd.doubleValue(), EphemerisElement.REDUCTION_METHOD.IAU_2009, true);
 
			if (eclipticMode) p = Ephem.equatorialToEcliptic(p, Constant.J2000, eph);
		}
		double f[] = p;
		for (int j=0;j<f.length;j++) {
			// Pass velocities to AU/s
			if (j > 2) f[j] /= Constant.SECONDS_PER_DAY;
			pos[i][j] = new Apfloat(f[j], precisionDigits);
		}
	}
 
	// Calculate the Gravitational constant from the header constants, and set 
	// consistent planetary masses
    G = new Apfloat(jpl.getConstant("GMS"), precisionDigits).divide(
					ApfloatMath.pow(new Apfloat(Constant.SECONDS_PER_DAY, precisionDigits), 2)); // AU^3/s^2, assuming Msun = 1
    TARGET.setPlanetaryMassFromJPLEphemeris(jpl);
 
    // Dirty tweaking factors fitted by hand for even greater accuracy. For speed <= 0.1 it is not necessary
    if (speed.doubleValue() == 10.0)
    	G = G.divide(new Apfloat(1.00035868, precisionDigits).divide(new Apfloat(1.00034708, precisionDigits))); // for speed = 10
    //G = G.divide(new Apfloat(1.00034824, precisionDigits).divide(new Apfloat(1.00034708, precisionDigits))); // for speed = 1
 
 
	// Applying DE405 initial conditions as documented in the DE405 explanatory supplement 
	// increases the residuals of the integration respect using previous values from JPL integration
	if (targets == targetsWithoutEarthMoon && jpl.getJPLVersion() == 405 && 
			jd.doubleValue() == 2440400.5) { // June 28 1969
		// Work in equatorial mode in this case to avoid losing precision when 
		// converting to ecliptic and comparing with JPL after some time
		eclipticMode = false;
		System.out.println("Applying DE405 initial conditions");
		String initialConditionsDE405[] = new String[] {
			"Mercury",
			"0.35726020644727541518 -0.09154904243051842990 -0.08598103998694037053",
			"0.00336784566219378527 0.02488934284224928990 0.01294407158679596663",
			"Venus",
			"0.60824943318560406033 -0.34913244319590053792 -0.19554434578540693592",
			"0.01095242010990883868 0.01561250673986247038 0.00632887645174666542",
			"Earth_Moon_Barycenter",
			"0.11601490913916648627 -0.92660555364038517604 -0.40180627760698804496",
			"0.01681162005220228976 0.00174313168798203152 0.00075597376713614610",
			"Mars",
			"-0.11468858243909270380 -1.32836653083348816476 -0.60615518941938081574",
			"0.01448200480794474793 0.00023728549236071136 -0.00028374983610239698",
			"Jupiter",
			"-5.38420940699214597622 -0.83124765616108382433 -0.22509475703354987777",
			"0.00109236329121849772 -0.00652329419119226767 -0.00282301226721943903",
			"Saturn",
			"7.88988993382281817537 4.59571072692601301962 1.55843151672508969735",
			"-0.00321720349109366378 0.00433063223355569175 0.00192641746379945286",
			"Uranus",
			"-18.26990081497826660524 -1.16271158021904696130 -0.25036954074255487461",
			"0.00022154016562741063 -0.00376765355824616179 -0.00165324380492239354",
			"Neptune",
			"-16.05954509192446441763 -23.94294829087950141524 -9.40042278035400838599",
			"0.00264312279157656145 -0.00150349208075879462 -0.00068127100487234772",
			"Pluto",
			"-30.48782211215555045830 -0.87324543019672926542 8.91129698418475509659",
			"0.00032256219593323324 -0.00314874792755160542 -0.00108017793159369583",
			"Sun",
			"0.00450250884530125842 0.00076707348146464055 0.00026605632781203556",
			"-0.00000035174423541454 0.00000517762777222281 0.00000222910220557907",
			"Moon",
			"-0.00080817732791148419 -0.00199463000162039941 -0.00108726266083810178",
			"0.00060108481665912983 -0.00016744546061515148 -0.00008556214497398616"
		};
		int sunIndex = -1;
		for (int i=0; i<initialConditionsDE405.length; i++) {
			String t = initialConditionsDE405[i].trim();
			TARGET target = null;
			try {
				target = TARGET.valueOf(t);
			} catch (Exception exc) {
				target = TARGET.valueOf(t.toUpperCase());
			}
			Apfloat spd = new Apfloat(Constant.SECONDS_PER_DAY, precisionDigits);
			for (int j=0;j<targets.length;j++) {
				if (targets[j] == target) {
					if (target == TARGET.SUN) sunIndex = j;
					System.out.println("   Found "+target.getName());
					String l1[] = DataSet.toStringArray(initialConditionsDE405[i+1], " ", true);
					String l2[] = DataSet.toStringArray(initialConditionsDE405[i+2], " ", true);
					pos[j][0] = new Apfloat(l1[0].trim(), precisionDigits);
					pos[j][1] = new Apfloat(l1[1].trim(), precisionDigits);
					pos[j][2] = new Apfloat(l1[2].trim(), precisionDigits);
					pos[j][3] = new Apfloat(l2[0].trim(), precisionDigits).divide(spd);
					pos[j][4] = new Apfloat(l2[1].trim(), precisionDigits).divide(spd);
					pos[j][5] = new Apfloat(l2[2].trim(), precisionDigits).divide(spd);
					break;
				}
			}
			i += 2;
		}
		for (int j=0;j<targets.length;j++) {
			if (targets[j] != TARGET.SUN && targets[j] != TARGET.Moon) { // Planets are heliocentric in the data above
				pos[j][0] = pos[j][0].add(pos[sunIndex][0]);
				pos[j][1] = pos[j][1].add(pos[sunIndex][1]);
				pos[j][2] = pos[j][2].add(pos[sunIndex][2]);
				pos[j][3] = pos[j][3].add(pos[sunIndex][3]);
				pos[j][4] = pos[j][4].add(pos[sunIndex][4]);
				pos[j][5] = pos[j][5].add(pos[sunIndex][5]);
			}
		}
	}
 
	return pos;
}

The previous code shows some differences besides the Apfloat implementation. I tested as initial conditions the values shown in the DE405 chapter 8 supplement. For some reason they show little differencies respect the results of the DE405 integration for the same date, so they are not adequate as initial conditions. Another difference are the 'dirty tweaking factors' I included to improve the results when using 10s as speed. They are minor changes to the G constant I fitted by hand to compensate for the numerical effects of using an excesive value for the speed, at least for the outer planets.

Here are the results for speed = 1s, no tweaking factors, and 10 days of integration time:

- Comparison with JPL DE405 (equatorial J2000 coordinates) -
JD (TDB): 2.4515545e6 (10.000 days from starting time)
Sun
x-y-z (JPL):       -0.007084986624  -0.002710890366  -0.000951626978
vx-vy-vz (JPL):    0.000000000063  -0.000000000077  -0.000000000035
x-y-z (NBody):     -0.007084986622  -0.002710890365  -0.000951626978
vx-vy-vz (NBody):  0.000000000063  -0.000000000077  -0.000000000035
Mercury
x-y-z (JPL):       0.073742795504  -0.396185425220  -0.219513620709
vx-vy-vz (JPL):    0.000000255220  0.000000075210  0.000000013718
x-y-z (NBody):     0.073742819595  -0.396185346988  -0.219513581420
vx-vy-vz (NBody):  0.000000255220  0.000000075210  0.000000013718
Venus
x-y-z (JPL):       -0.692409955824  -0.220934911643  -0.055755857812
vx-vy-vz (JPL):    0.000000070655  -0.000000202548  -0.000000095591
x-y-z (NBody):     -0.692409920864  -0.220934909566  -0.055755859090
vx-vy-vz (NBody):  0.000000070655  -0.000000202548  -0.000000095591
Earth-Moon barycenter
x-y-z (JPL):       -0.344426805301  0.844813611199  0.366494228132
vx-vy-vz (JPL):    -0.000000190223  -0.000000063430  -0.000000027501
x-y-z (NBody):     -0.344426803783  0.844813596258  0.366494221857
vx-vy-vz (NBody):  -0.000000190223  -0.000000063430  -0.000000027501
Mars
x-y-z (JPL):       1.383134535178  0.129701457011  0.022192225611
vx-vy-vz (JPL):    -0.000000008884  0.000000159005  0.000000073173
x-y-z (NBody):     1.383134526142  0.129701457033  0.022192225865
vx-vy-vz (NBody):  -0.000000008884  0.000000159005  0.000000073173
Jupiter
x-y-z (JPL):       3.950257913739  2.789441305836  1.099449361810
vx-vy-vz (JPL):    -0.000000053869  0.000000067260  0.000000030142
x-y-z (NBody):     3.950257913177  2.789441305452  1.099449361659
vx-vy-vz (NBody):  -0.000000053869  0.000000067260  0.000000030142
Saturn
x-y-z (JPL):       6.358438049285  6.205358229132  2.289379644655
vx-vy-vz (JPL):    -0.000000049887  0.000000040498  0.000000018872
x-y-z (NBody):     6.358438049143  6.205358228994  2.289379644605
vx-vy-vz (NBody):  -0.000000049887  0.000000040498  0.000000018872
Uranus
x-y-z (JPL):       14.450191832477  -12.485560799374  -5.672740512421
vx-vy-vz (JPL):    0.000000030999  0.000000028469  0.000000012030
x-y-z (NBody):     14.450191832445  -12.485560799347  -5.672740512408
vx-vy-vz (NBody):  0.000000030999  0.000000028469  0.000000012030
Neptune
x-y-z (JPL):       16.829465510886  -22.966959644683  -9.819493012110
vx-vy-vz (JPL):    0.000000029895  0.000000019260  0.000000007139
x-y-z (NBody):     16.829465510875  -22.966959644669  -9.819493012104
vx-vy-vz (NBody):  0.000000029895  0.000000019260  0.000000007139
Maximum difference (AU): 7.823123476802252E-8 (for Mercury)
Maximum difference from Earth ("): 0.015260370588755635 (for Mercury)
Integration took 1497.053 s

And here are the results for speed = 10s, applying the tweak to G, and again 10 days of integration:

- Comparison with JPL DE405 (equatorial J2000 coordinates) -
JD (TDB): 2.4515545e6 (10.000 days from starting time)
Sun
x-y-z (JPL):       -0.007084986624  -0.002710890366  -0.000951626978
vx-vy-vz (JPL):    0.000000000063  -0.000000000077  -0.000000000035
x-y-z (NBody):     -0.007084986623  -0.002710890366  -0.000951626978
vx-vy-vz (NBody):  0.000000000063  -0.000000000077  -0.000000000035
Mercury
x-y-z (JPL):       0.073742795504  -0.396185425220  -0.219513620709
vx-vy-vz (JPL):    0.000000255220  0.000000075210  0.000000013718
x-y-z (NBody):     0.073742918414  -0.396185406264  -0.219513623332
vx-vy-vz (NBody):  0.000000255220  0.000000075208  0.000000013717
Venus
x-y-z (JPL):       -0.692409955824  -0.220934911643  -0.055755857812
vx-vy-vz (JPL):    0.000000070655  -0.000000202548  -0.000000095591
x-y-z (NBody):     -0.692409947626  -0.220934938543  -0.055755870433
vx-vy-vz (NBody):  0.000000070654  -0.000000202548  -0.000000095591
Earth-Moon barycenter
x-y-z (JPL):       -0.344426805301  0.844813611199  0.366494228132
vx-vy-vz (JPL):    -0.000000190223  -0.000000063430  -0.000000027501
x-y-z (NBody):     -0.344426816935  0.844813609856  0.366494227752
vx-vy-vz (NBody):  -0.000000190223  -0.000000063429  -0.000000027501
Mars
x-y-z (JPL):       1.383134535178  0.129701457011  0.022192225611
vx-vy-vz (JPL):    -0.000000008884  0.000000159005  0.000000073173
x-y-z (NBody):     1.383134534759  0.129701459906  0.022192226950
vx-vy-vz (NBody):  -0.000000008884  0.000000159005  0.000000073173
Jupiter
x-y-z (JPL):       3.950257913739  2.789441305836  1.099449361810
vx-vy-vz (JPL):    -0.000000053869  0.000000067260  0.000000030142
x-y-z (NBody):     3.950257913724  2.789441305867  1.099449361823
vx-vy-vz (NBody):  -0.000000053869  0.000000067260  0.000000030142
Saturn
x-y-z (JPL):       6.358438049285  6.205358229132  2.289379644655
vx-vy-vz (JPL):    -0.000000049887  0.000000040498  0.000000018872
x-y-z (NBody):     6.358438049285  6.205358229137  2.289379644658
vx-vy-vz (NBody):  -0.000000049887  0.000000040498  0.000000018872
Uranus
x-y-z (JPL):       14.450191832477  -12.485560799374  -5.672740512421
vx-vy-vz (JPL):    0.000000030999  0.000000028469  0.000000012030
x-y-z (NBody):     14.450191832477  -12.485560799375  -5.672740512421
vx-vy-vz (NBody):  0.000000030999  0.000000028469  0.000000012030
Neptune
x-y-z (JPL):       16.829465510886  -22.966959644683  -9.819493012110
vx-vy-vz (JPL):    0.000000029895  0.000000019260  0.000000007139
x-y-z (NBody):     16.829465510886  -22.966959644683  -9.819493012110
vx-vy-vz (NBody):  0.000000029895  0.000000019260  0.000000007139
Maximum difference (AU): 1.2290968719186335E-7 (for Mercury)
Maximum difference from Earth ("): 0.019564751377336518 (for Mercury)
Integration took 164.227 s

Thanks to the tweak the accuracy with 10s speed is quite similar with only 10^-7 AU (15 km) of error in the worst case for Mercury (without it would be almost 10x more). Notice also the perfect match for the outer planets Uranus and Neptune, not achieved with 1s speed unless you also apply the tweak to G. As an exercise I suggest to patiently try a speed of 0.1s. You will see again a perfect match for the outer planets (but with no tweak) and the worst result again for Mercury, with a position error of 10^-8 AU (2 km only!).

As you see these results are great, but the program requires a lot of time to perform the integration. It is interesting to see the great results of the Newton's law applied directly without any other consideration, but for ephemerides calculations it is too slow, so it is mandatory to use double precision only. Another point to discuss would be the mathematical approach. The point here is to use a direct approach, which is simple to implement but slow when looking for high accuracies. For that purpose there are more practical options like Runge-Kutta integrators as we will see soon.

The direct approach of the Newton method presented before is inadequate for several reasons. The velocities are updated first and then the positions using a very rude method which is not continuous at all. This produces inaccuracies that are summed over the time. Of course in the limit speed = 0 the method works quite well, but in that limit we have also rounding errors or slow execution in extended precision.

There are many techniques to solve the differential equation (for dv/dt) in a more consistent way for all the bodies in a nbody problem. After testing several of them like 4th order Runge-Kutta or even more complex ones, I found that the modified 2nd order Euler method (also known as Heun's method) is probably the best for the problem presented in this post, since the accuracy increases a lot and the performance is also very good. This method estimates the forces in the next iteration to go back to the current one to correct the velocities and positions using an averaged value with the next iteration. For an adequate explanation with clear formulas you can go to the Wikipedia with the previous link. I prefer to present directly the code with some comments so you will get the idea.

private static void Newton() throws Exception {
	// Integration for next drawing, credits for Newton and his apple
	double gravity[][] = new double[targets.length][3];
	for (int i=0;i<targets.length-1;i++) {
		for (int j=i+1; j<targets.length;j++) {
			double dx = pos[i][0] - pos[j][0];
			double dy = pos[i][1] - pos[j][1];
			double dz = pos[i][2] - pos[j][2];
			double d2 = dx*dx + dy*dy + dz*dz;
			double d = FastMath.sqrt(d2);
 
			if (targets[j].relativeMass != 0) {
				double g = -G / (d2 * targets[j].relativeMass);
				if (!neglectWeakInteractions || g < 1000) { // Neglect interactions from far/little bodies
					gravity[i][0] += g * dx / d;
					gravity[i][1] += g * dy / d;
					gravity[i][2] += g * dz / d;
				}
			}
			if (targets[i].relativeMass != 0) {
				double g = -G / (d2 * targets[i].relativeMass);
				if (!neglectWeakInteractions || g < 1000) { // Neglect interactions from far/little bodies
					gravity[j][0] -= g * dx / d;
					gravity[j][1] -= g * dy / d;
					gravity[j][2] -= g * dz / d;
				}
			}
		} 
	}
 
	/* Compute poor estimates for future positions and velocities */
	double poorNewPosVel[][] = new double[targets.length][6];
	for (int i=0;i<targets.length;i++) {	
		poorNewPosVel[i][3] = pos[i][3] + gravity[i][0] * speed;
		poorNewPosVel[i][4] = pos[i][4] + gravity[i][1] * speed;
		poorNewPosVel[i][5] = pos[i][5] + gravity[i][2] * speed;
 
		poorNewPosVel[i][0] = pos[i][0] + poorNewPosVel[i][3] * speed;
		poorNewPosVel[i][1] = pos[i][1] + poorNewPosVel[i][4] * speed;
		poorNewPosVel[i][2] = pos[i][2] + poorNewPosVel[i][5] * speed;
	}
    /* Use previous poor future positions to predict poor future accelerations for all objects */
	double futurePosVel[][] = new double[targets.length][6];
	for (int i = 0; i < targets.length; i++) {
		double poorNewGravity[] = new double[3];
		for (int j = 0; j < targets.length; j++) {
			if (i == j) continue;
 
			double dxyz[] = new double[] {
				poorNewPosVel[i][0] - poorNewPosVel[j][0],
				poorNewPosVel[i][1] - poorNewPosVel[j][1],
				poorNewPosVel[i][2] - poorNewPosVel[j][2]
			};
 
			double d2 = Functions.scalarProduct(dxyz, dxyz);
			double d = Math.sqrt(d2);
			double newGravity = -G / (d2 * targets[j].relativeMass);
			poorNewGravity[0] += newGravity * dxyz[0] / d;
			poorNewGravity[1] += newGravity * dxyz[1] / d;
			poorNewGravity[2] += newGravity * dxyz[2] / d;
		}
 
		/* Compute average of current and future accelerations */
		double averageGravity[] = new double[] {
				(gravity[i][0] + poorNewGravity[0]) * 0.5,
				(gravity[i][1] + poorNewGravity[1]) * 0.5,
				(gravity[i][2] + poorNewGravity[2]) * 0.5
		};
 
		/* Use average accelerations to compute better future velocities */
		double betterNewVelocity[] = new double[] {
				pos[i][3] + averageGravity[0] * speed,
				pos[i][4] + averageGravity[1] * speed,
				pos[i][5] + averageGravity[2] * speed,
		};
 
		/* Compute average of current and future velocities */
		double averageVelocity[] = new double[] {
				(pos[i][3] + betterNewVelocity[0]) * 0.5,
				(pos[i][4] + betterNewVelocity[1]) * 0.5,
				(pos[i][5] + betterNewVelocity[2]) * 0.5,
		};
 
		/* Use average velocity to compute better future position */
		double betterNewPosition[] = new double[] {
				pos[i][0] + averageVelocity[0] * speed,
				pos[i][1] + averageVelocity[1] * speed,
				pos[i][2] + averageVelocity[2] * speed,
		};
		futurePosVel[i] = DataSet.addDoubleArray(betterNewPosition, betterNewVelocity);
	}
	pos = futurePosVel;
 
	// Compute L
	angularMomentum = 0;
	for (int i=0;i<targets.length;i++) {	
		if (showLabels && targets[i].relativeMass > 0) {
			DoubleVector vr = new DoubleVector(pos[i][0], pos[i][1], pos[i][2]);
			DoubleVector vv = new DoubleVector(pos[i][3], pos[i][4], pos[i][5]);
			angularMomentum += vr.crossProduct(vv).norm2() / targets[i].relativeMass;
		}
	}
 
	// Make integration effective on date
	jd += speed / Constant.SECONDS_PER_DAY;
}

This new method is far more accurate and stable when using high values of speed (unless you use huge values and make the inner planets crash). With speed = 200s the results of the previous numerical integration test now becomes:

- Comparison with JPL DE405 (equatorial J2000 coordinates) -
JD (TDB): 2451724.5 (180.000 days from starting time)
Sun
x-y-z (JPL):       -0.006032658552  -0.003741542907  -0.001420728930
vx-vy-vz (JPL):    0.000000000080  -0.000000000062  -0.000000000029
x-y-z (NBody):     -0.006032658079  -0.003741542799  -0.001420728920
vx-vy-vz (NBody):  0.000000000080  -0.000000000062  -0.000000000029
Mercury
x-y-z (JPL):       -0.058103868552  -0.414395272621  -0.215375832955
vx-vy-vz (JPL):    0.000000258350  -0.000000008197  -0.000000031160
x-y-z (NBody):     -0.058103543150  -0.414395275938  -0.215375868482
vx-vy-vz (NBody):  0.000000258350  -0.000000008197  -0.000000031160
Venus
x-y-z (JPL):       -0.242561212403  0.609559831691  0.289464459217
vx-vy-vz (JPL):    -0.000000221752  -0.000000076642  -0.000000020442
x-y-z (NBody):     -0.242561563368  0.609559689232  0.289464417339
vx-vy-vz (NBody):  -0.000000221752  -0.000000076642  -0.000000020442
Earth-Moon barycenter
x-y-z (JPL):       0.127769131873  -0.928376971857  -0.402297058027
vx-vy-vz (JPL):    0.000000194230  0.000000023296  0.000000010098
x-y-z (NBody):     0.127768908839  -0.928377054611  -0.402297109777
vx-vy-vz (NBody):  0.000000194230  0.000000023296  0.000000010098
Mars
x-y-z (JPL):       -0.250313769347  1.422640449575  0.659416143757
vx-vy-vz (JPL):    -0.000000153806  -0.000000011682  -0.000000001198
x-y-z (NBody):     -0.250313807781  1.422640406519  0.659416125053
vx-vy-vz (NBody):  -0.000000153806  -0.000000011682  -0.000000001198
Jupiter
x-y-z (JPL):       3.031599995432  3.669672591247  1.499124814774
vx-vy-vz (JPL):    -0.000000070447  0.000000051951  0.000000023984
x-y-z (NBody):     3.031599994671  3.669672590569  1.499124814496
vx-vy-vz (NBody):  -0.000000070447  0.000000051951  0.000000023984
Saturn
x-y-z (JPL):       5.591660709347  6.764515860547  2.553313687645
vx-vy-vz (JPL):    -0.000000054434  0.000000035560  0.000000017030
x-y-z (NBody):     5.591660709250  6.764515860462  2.553313687611
vx-vy-vz (NBody):  -0.000000054434  0.000000035560  0.000000017030
Uranus
x-y-z (JPL):       14.897614837489  -12.060738395909  -5.493007971301
vx-vy-vz (JPL):    0.000000029919  0.000000029372  0.000000012441
x-y-z (NBody):     14.897614837516  -12.060738395880  -5.493007971298
vx-vy-vz (NBody):  0.000000029919  0.000000029372  0.000000012441
Neptune
x-y-z (JPL):       17.265901881373  -22.680489235848  -9.713104446249
vx-vy-vz (JPL):    0.000000029531  0.000000019747  0.000000007347
x-y-z (NBody):     17.265901881418  -22.680489235838  -9.713104446255
vx-vy-vz (NBody):  0.000000029531  0.000000019747  0.000000007347
Maximum difference (AU): 3.5096551720403824E-7 (for Venus)
Maximum difference from Earth ("): 0.19796559374329475 (for Mercury)
Integration took 0.363 s

Which is around 10x more accurate (similar to the Apfloat implementation) and 3x faster than the direct approach, thanks to the increase in the speed value. Obviously this is the best method for all purposes, both for graphics and numerical integration, and can be used with the previous suggested custom ephemerides files to obtain ephemerides with great accuracy (with 200s of speed my tests show a peak error and mean error of 0.25” and 0.01” for Mercury, lower for the other bodies) and very good performace. In case you enable the Earth and Moon as separate bodies the peak error after 180 days will be around 40 km only. And by including Pluto the outer planets will reproduce DE405 almost perfectly. The 4th order Runge-Kutta is almost 2x slower and, overall, only slightly more accurate, although it is particularly better for Mercury and maybe the Moon. With discrepancies of 3 10^-7 AU for Mercury after one year, this method goes close to the limits without applying other complex corrections, since the precession of the perihelion of Mercury due to relativistic effects is 43” per century, equivalent to a deviation of almost 10^-7 AU per year.

As a final test, I leave here the results obtained with 200s speed and 180 days of integration when including Pluto and the three main asteroids. The test was executed with a slightly revised version of the Heun's method, in which I basically repeat the estimates considering the array futurePosVel again as a new (less poor) poorNewPosVel. This modification seems to return better results than the 4th order Runge-Kutta while maintaining a similar performance and a more simple implementation. Notice the improvement in the position of both the Sun and the outer planets (with discrepancies in the cm range!), despite the asteroid initial positions come from approximate elliptic elements:

- Comparison with JPL DE405 (equatorial J2000 coordinates) -
JD (TDB): 2451724.5 (180.000 days from starting time)
Sun
x-y-z (JPL):       -0.006032658552  -0.003741542907  -0.001420728930
vx-vy-vz (JPL):    0.000000000080  -0.000000000062  -0.000000000029
x-y-z (NBody):     -0.006032658535  -0.003741542885  -0.001420728913
vx-vy-vz (NBody):  0.000000000080  -0.000000000062  -0.000000000029
Mercury
x-y-z (JPL):       -0.058103868552  -0.414395272621  -0.215375832955
vx-vy-vz (JPL):    0.000000258350  -0.000000008197  -0.000000031160
x-y-z (NBody):     -0.058103685495  -0.414395268228  -0.215375849594
vx-vy-vz (NBody):  0.000000258350  -0.000000008197  -0.000000031160
Venus
x-y-z (JPL):       -0.242561212403  0.609559831691  0.289464459217
vx-vy-vz (JPL):    -0.000000221752  -0.000000076642  -0.000000020442
x-y-z (NBody):     -0.242561547063  0.609559695939  0.289464419335
vx-vy-vz (NBody):  -0.000000221752  -0.000000076642  -0.000000020442
Earth
x-y-z (JPL):       0.127753535978  -0.928400923101  -0.402304963432
vx-vy-vz (JPL):    0.000000194306  0.000000023259  0.000000010077
x-y-z (NBody):     0.127753710193  -0.928400842769  -0.402304928419
vx-vy-vz (NBody):  0.000000194306  0.000000023259  0.000000010077
Moon
x-y-z (JPL):       0.129037086833  -0.926429722287  -0.401654344187
vx-vy-vz (JPL):    0.000000188025  0.000000026306  0.000000011825
x-y-z (NBody):     0.129037319867  -0.926429666202  -0.401654334317
vx-vy-vz (NBody):  0.000000188026  0.000000026306  0.000000011825
Mars
x-y-z (JPL):       -0.250313769347  1.422640449575  0.659416143757
vx-vy-vz (JPL):    -0.000000153806  -0.000000011682  -0.000000001198
x-y-z (NBody):     -0.250313807503  1.422640406982  0.659416125273
vx-vy-vz (NBody):  -0.000000153806  -0.000000011682  -0.000000001198
Jupiter
x-y-z (JPL):       3.031599995432  3.669672591247  1.499124814774
vx-vy-vz (JPL):    -0.000000070447  0.000000051951  0.000000023984
x-y-z (NBody):     3.031599994598  3.669672590513  1.499124814480
vx-vy-vz (NBody):  -0.000000070447  0.000000051951  0.000000023984
Saturn
x-y-z (JPL):       5.591660709347  6.764515860547  2.553313687645
vx-vy-vz (JPL):    -0.000000054434  0.000000035560  0.000000017030
x-y-z (NBody):     5.591660709218  6.764515860424  2.553313687600
vx-vy-vz (NBody):  -0.000000054434  0.000000035560  0.000000017030
Uranus
x-y-z (JPL):       14.897614837489  -12.060738395909  -5.493007971301
vx-vy-vz (JPL):    0.000000029919  0.000000029372  0.000000012441
x-y-z (NBody):     14.897614837475  -12.060738395899  -5.493007971296
vx-vy-vz (NBody):  0.000000029919  0.000000029372  0.000000012441
Neptune
x-y-z (JPL):       17.265901881373  -22.680489235848  -9.713104446249
vx-vy-vz (JPL):    0.000000029531  0.000000019747  0.000000007347
x-y-z (NBody):     17.265901881370  -22.680489235844  -9.713104446247
vx-vy-vz (NBody):  0.000000029531  0.000000019747  0.000000007347
Pluto
x-y-z (JPL):       -9.336180264872  -28.180274953790  -5.981237412137
vx-vy-vz (JPL):    0.000000035331  -0.000000012505  -0.000000014547
x-y-z (NBody):     -9.336180264871  -28.180274953785  -5.981237412136
vx-vy-vz (NBody):  0.000000035331  -0.000000012505  -0.000000014547
Maximum difference (AU): 3.346603522913494E-7 (for Venus)
Integration took 0.876 s

Can these results be improved even more, without having to implement non-linear relativistic corrections to gravity or the effects of the figure of the Earth, Moon, and other bodies ? Mmmmmmm… Implementing the latest Newton method with Apfloat and using lower values of speed (again around 20 instead of 200) could certainly help, especially in long-term integrations, but the performance drop prevents any practical use, would be only for academic interest. However, the major improvement probably will come from a higher order integrator. Another problem with the previous implementation is that I consider at most three asteroids. In long-term integrations the inner planets will always be affected more than the outer ones by their mutual perturbations. As an exercise I leave to you the task to repeat the tests including the computation of 300 asteroides using the initial conditions detailed at the end of this DE430 report. You will have to first pass the coordinates to barycentric.

Well, it seems this post is going beyond my initial intentions, we will cross the line only a little… The DE405 explanatory supplement chapter 8 presents the relativistic corrections in equation 8-1. Most of these terms are difficult to quantify, while others are easy. I present below a little code you can copy to the top of the Newton method to calculate a partial relativistic correction relCor[] for each planet, a value which is slightly below 1. Later, in the gravity computation, where the uppercase constant 'G' is used, you can multiply it by relCor[i] to compute the corrected value of gravity (both in the first loop and also after in the Euler method. In the first loop use relCor[j] for the second symmetric case). The code for the partial correction is:

double relCor[] = new double[targets.length];
double c2 = Math.pow(Constant.SPEED_OF_LIGHT * 0.001 / Constant.AU, 2);
for (int i=0;i<targets.length;i++) {
	relCor[i] = 0;
	for (int j=0; j<targets.length;j++) {
		if (j == i) continue;
		double dx = pos[i][0] - pos[j][0];
		double dy = pos[i][1] - pos[j][1];
		double dz = pos[i][2] - pos[j][2];
		double d2 = dx*dx + dy*dy + dz*dz;
		double d = FastMath.sqrt(d2);
		relCor[i] += G / (d * targets[j].relativeMass);
	}
	relCor[i] = -relCor[i] * 4.0 / c2;
	double dx = pos[i][3];
	double dy = pos[i][4];
	double dz = pos[i][5];
	double d2 = dx*dx + dy*dy + dz*dz;
	relCor[i] += (d2 / c2) + 1.0;
}

With this correction included and changing the speed to 450s, which according to my tests is the best value when including this correction, the discrepancy with the DE4xx integration goes down dramatically from 4E-7 AU, as you can see in the new results for the previous example:

- Comparison with JPL DE405 (equatorial J2000 coordinates) -
JD (TDB): 2451724.5 (180.000 days from starting time)
Sun
x-y-z (JPL):       -0.006032658552  -0.003741542907  -0.001420728930
vx-vy-vz (JPL):    0.799550639E-10  -0.619466189E-10  -0.286648994E-10
x-y-z (NBody):     -0.006032658536  -0.003741542885  -0.001420728913
vx-vy-vz (NBody):  0.799550655E-10  -0.619466158E-10  -0.286648972E-10
Mercury
x-y-z (JPL):       -0.058103868552  -0.414395272621  -0.215375832955
vx-vy-vz (JPL):    2583.499231E-10  -81.971220E-10  -311.602330E-10
x-y-z (NBody):     -0.058103958263  -0.414395265729  -0.215375819970
vx-vy-vz (NBody):  2583.499211E-10  -81.971581E-10  -311.602521E-10
Venus
x-y-z (JPL):       -0.242561212403  0.609559831691  0.289464459217
vx-vy-vz (JPL):    -2217.515527E-10  -766.416577E-10  -204.418962E-10
x-y-z (NBody):     -0.242561210712  0.609559831861  0.289464459198
vx-vy-vz (NBody):  -2217.515532E-10  -766.416569E-10  -204.418959E-10
Earth
x-y-z (JPL):       0.127753535978  -0.928400923101  -0.402304963432
vx-vy-vz (JPL):    1943.064760E-10  232.586664E-10  100.766664E-10
x-y-z (NBody):     0.127753534173  -0.928400920896  -0.402304962283
vx-vy-vz (NBody):  1943.064746E-10  232.586648E-10  100.766655E-10
Moon
x-y-z (JPL):       0.129037086833  -0.926429722287  -0.401654344187
vx-vy-vz (JPL):    1880.2545301E-10  263.0551844E-10  118.2473451E-10
x-y-z (NBody):     0.129037152492  -0.926429746777  -0.401654369917
vx-vy-vz (NBody):  1880.2557446E-10  263.0566905E-10  118.2482177E-10
Mars
x-y-z (JPL):       -0.250313769347  1.422640449575  0.659416143757
vx-vy-vz (JPL):    -1538.060447E-10  -116.818862E-10  -11.979152E-10
x-y-z (NBody):     -0.250313767401  1.422640446134  0.659416142147
vx-vy-vz (NBody):  -1538.060445E-10  -116.818867E-10  -11.979154E-10
Jupiter
x-y-z (JPL):       3.031599995432  3.669672591247  1.499124814774
vx-vy-vz (JPL):    -704.467905E-10  519.514300E-10  239.840055E-10
x-y-z (NBody):     3.031599995459  3.669672591221  1.499124814762
vx-vy-vz (NBody):  -704.467905E-10  519.514300E-10  239.840055E-10
Saturn
x-y-z (JPL):       5.591660709347  6.764515860547  2.553313687645
vx-vy-vz (JPL):    -544.340531E-10  355.602819E-10  170.295232E-10
x-y-z (NBody):     5.591660709341  6.764515860551  2.553313687648
vx-vy-vz (NBody):  -544.340531E-10  355.602819E-10  170.295232E-10
Uranus
x-y-z (JPL):       14.897614837489  -12.060738395909  -5.493007971301
vx-vy-vz (JPL):    299.1921531E-10  293.7190011E-10  124.4093746E-10
x-y-z (NBody):     14.897614837488  -12.060738395910  -5.493007971301
vx-vy-vz (NBody):  299.1921531E-10  293.7190011E-10  124.4093746E-10
Neptune
x-y-z (JPL):       17.265901881373  -22.680489235848  -9.713104446249
vx-vy-vz (JPL):    295.3110154E-10  197.4679538E-10  73.4727591E-10
x-y-z (NBody):     17.265901881373  -22.680489235848  -9.713104446249
vx-vy-vz (NBody):  295.3110154E-10  197.4679538E-10  73.4727591E-10
Pluto
x-y-z (JPL):       -9.336180264872  -28.180274953790  -5.981237412137
vx-vy-vz (JPL):    353.3082307E-10  -125.0479803E-10  -145.4736878E-10
x-y-z (NBody):     -9.336180264872  -28.180274953789  -5.981237412137
vx-vy-vz (NBody):  353.3082307E-10  -125.0479803E-10  -145.4736878E-10
Maximum difference (AU): 8.971095644821903E-8 (for Mercury)
Integration took 0.599 s

These last results were obtained as the previous ones, including the three main asteroids and Pluto, and shows again the Earth and Moon as separated bodies. It is notorious that all objects improve a lot except the Sun. As you see after 180 days of integration the discrepancy in the inner planets except Mercury/Moon is below the km level, few meters for the massive bodies Jupiter and the Sun, and in the cm level beyond Saturn. From Earth this means to leave the lunar position discrepancy in few arcseconds and all the planets except Mercury inside the milliarcsecond range. The position of Mercury improves with both the speed change and the relativistic correction, and the discrepancy seems to be just below and close to 1E-7 AU most of the times (remember this is the expected accuracy limit due to the precession of the perihelion of Mercury). Another nice result are the discrepancies after five years: they increase only to few km and they keep relatively under control in Mercury. In longer integrations like 25 years the position of Mercury and the Moon degrades about 400 km (for Mercury it is the expected value when accumulating the effects of the precession of the perihelion), while Mars about 150 km, and Earth, Venus, and Jupiter less than 30 km. Mars is more affected and Venus less, suggesting more asteroids are needed in the model, although my tests seems to indicate more asteroids will only improve the position of the Sun.

Now it is evident we reached the limits of a reasonably simple implementation, more accuracy means a more robust integrator, extended precision, more asteroids, and especially a complete and rigurous implementation of relativistic corrections and the figure of the Earth, Sun, and Moon tides.

So the conclusion of this post is that the accuracy limit when following directly and strictly the Newton's law ranges from around the meter level in outer planets to around or below 100 km for inner ones for integration periods of 180 days, although the accuracy depends on using extended precision and reducing the integration step. Far better and faster results are obtained with an integrator, and among the possible options the implementation chosen as more adequate (simple, educative, with good accuracy and performance) is the modified 2nd order Euler integrator (Heun's method) in double precision, but slightly modified by repeating it in 2-pass mode. Discrepancies with JPL ephemerides goes down to the cm range for external planets when including Pluto in the integration, and the inner planets and the Sun also improve to few km when including the three main asteroids. Even in double precision arithmetics the accuracy can be increased about 5x more with a partial relativistic correction and a speed value of 450s, leaving the discrepancy below 1 km after 180 days in almost all cases. However, Mercury and the Moon are special cases and, without an specific rigorous model for them (considering the figures of their perturbing bodies Sun and Earth), the discrepancy in these bodies cannot be consistently reduced below the 10-15 km level (peak errors of 0.04” and 10” from Earth, accumulating with time). Considering more asteroids does not help significantly for the short integration periods analyzed. Although other minor improvements are still feasible, the program de118i-2 by S. Moshier mentioned in the introduction provides results several orders of magnitude better, so it is the implementation to follow when more accuracy is required.

As a final practical application I will take one of the latest JPL integrations to reduce it to a little binary file including the GMX constants and the initial conditions once per year during several centuries. For this purpose I wrote the following piece of code:

public static void planetaryEphemerisWriter() throws Exception {
	TARGET jplTarget[] = new TARGET[] {
			TARGET.NOT_A_PLANET, TARGET.MERCURY, TARGET.VENUS, TARGET.EARTH, TARGET.MARS, 
			TARGET.JUPITER, TARGET.SATURN, TARGET.URANUS, TARGET.NEPTUNE, TARGET.Pluto, 
			TARGET.Moon, TARGET.SUN, TARGET.Solar_System_Barycenter, TARGET.Earth_Moon_Barycenter, 
			TARGET.Nutation, TARGET.Libration
	};
	jplVersion = ALGORITHM.JPL_DE430;
	targets = new TARGET[] {TARGET.SUN, TARGET.MERCURY, TARGET.VENUS,
			TARGET.EARTH, TARGET.Moon, TARGET.MARS, TARGET.JUPITER, TARGET.SATURN, TARGET.URANUS,
			TARGET.NEPTUNE, TARGET.Pluto};
	eclipticMode = true;
	startTime = new AstroDate(1850, 1, 1, 0, 0, 0);
	AstroDate endTime = new AstroDate(2150, 1, 1, 0, 0, 0);
	double step = 360;
 
	jd0 = startTime.jd();
	double jdf = endTime.jd();
	int n = (int) ((jdf - jd0) / step);
	jdf = jd0 + step * n;
	jpl = new JPLEphemeris(jplVersion);
	jd = jd0;
 
	String path = "/home/alonso/";
	File file = new File(path + "de"+jpl.getJPLVersion()+".jparsec");
	RandomAccessFile f = new RandomAccessFile(file, "rw");
	f.writeInt(jpl.getJPLVersion());
	f.writeDouble(jd0);
	f.writeDouble(jdf);
	f.writeDouble(step);
	f.writeInt(targets.length);
	for (int i=0;i<targets.length;i++) {
		int index = -1;
		for (int j=0;j<jplTarget.length;j++) {
			if (jplTarget[j] == targets[i]) {
				index = j;
				break;
			}
		}
		f.writeInt(index);
	}
	for (int i=0;i<targets.length;i++) {
		if (targets[i] == TARGET.SUN) {
	        double gms = Double.parseDouble(jpl.getConstant("GMS"));
	        f.writeDouble(gms);
		} else {
			if (targets[i] == TARGET.Earth_Moon_Barycenter) {
		        double gmb = Double.parseDouble(jpl.getConstant("GMB"));
		        f.writeDouble(gmb);					
			} else {
				if (targets[i] == TARGET.EARTH || targets[i] == TARGET.Moon) {
			        double emrat = Double.parseDouble(jpl.getConstant("EMRAT"));
			        double gmb = Double.parseDouble(jpl.getConstant("GMB"));
			        double gme = gmb / (1.0 + 1.0 / emrat);
			        double gmm = gmb / (1.0 + emrat);
			        double gm = targets[i] == TARGET.EARTH ? gme : gmm;
			        f.writeDouble(gm);
				} else {
			        double gm = Double.parseDouble(jpl.getConstant("GM"+targets[i].ordinal()));
			        f.writeDouble(gm);
				}
			}
		}
	}
	f.writeDouble(Double.parseDouble(jpl.getConstant("AU")));
	for (int i=0;i<=n;i++) {
		pos = getJPLPositions(jpl, eph, jd);
		for (int j=0;j<pos.length;j++) {
			for (int k=0;k<pos[j].length;k++) {
				f.writeDouble(pos[j][k]);
			}
		}
		jd += step;
	}
	f.close();
}

This method writes to a binary file the DE version, the start and end times for the ephemeris file, and the time step in days between the records. Then it writes the number of targets considered, and for each of them it writes the integer identifier used by the JPL for that body. Then, again for each body, it writes the GM constant as provided in the header file of the JPL integration. The AU constant is the only additional value written. Finally, the set of positions and velocities for each date is written in ecliptic coordinates.

The program works between 1850 and 2150, producing a binary file of only 150 kb. The corresponding reader and a comparison with the JPL integration would be something like this:

public static void planetaryEphemerisReader() throws Exception {
	double jdCalc = Constant.J2000 - 170 + 360*5;
	jplVersion = ALGORITHM.JPL_DE430;
	jpl = new JPLEphemeris(jplVersion);
	TARGET jplTarget[] = new TARGET[] {
			TARGET.NOT_A_PLANET, TARGET.MERCURY, TARGET.VENUS, TARGET.EARTH, TARGET.MARS, 
			TARGET.JUPITER, TARGET.SATURN, TARGET.URANUS, TARGET.NEPTUNE, TARGET.Pluto, 
			TARGET.Moon, TARGET.SUN, TARGET.Solar_System_Barycenter, TARGET.Earth_Moon_Barycenter, 
			TARGET.Nutation, TARGET.Libration
	};
	eclipticMode = true;
 
	long t0 = System.currentTimeMillis();
	String path = "/home/alonso/";
	File file = new File(path + "de"+jpl.getJPLVersion()+".jparsec");
	RandomAccessFile f = new RandomAccessFile(file, "rw");
	int de = f.readInt();
	jd0 = f.readDouble();
	double jdf = f.readDouble();
	double step = f.readDouble();
	int nTarget = f.readInt();
	targets = new TARGET[nTarget];
	int sunIndex = -1;
	for (int i=0;i<nTarget;i++) {
		targets[i] = jplTarget[f.readInt()];
		if (targets[i] == TARGET.SUN) sunIndex = i;
		if (targets[i] == TARGET.NOT_A_PLANET) {
			f.close();
			throw new Exception("Incorrect target identifier");				
		}
	}
	double gms[] = new double[nTarget];
	for (int i=0;i<nTarget;i++) {
        gms[i] = f.readDouble();
	}
	double au = f.readDouble();
	for (int i=0;i<nTarget;i++) {
		targets[i].relativeMass = gms[sunIndex] / gms[i];
	}
	if (de != jpl.getJPLVersion()) {
		f.close();
		throw new Exception("Incorrect file or integration version");
	}
	int recordSkip = (int) ((jdCalc - jd0) / step + 0.5);
	int recordSkipMax = (int) ((jdf - jd0) / step - 0.5);
	if (recordSkip < 0 || recordSkip > recordSkipMax) {
		f.close();
		throw new Exception("Incorrect date, out of range"); 			
	}
	int bytesPerRecord = 8 * nTarget * 6;
	int n = f.skipBytes(bytesPerRecord * recordSkip);
	if (n != bytesPerRecord * recordSkip) {
		f.close();
		throw new Exception("Could not find the record, corrupt file ?");
	}
	pos = new double[nTarget][6];
	for (int i=0;i<nTarget;i++) {
		for (int j=0;j<6;j++) {
			pos[i][j] = f.readDouble();
		}
	}
	f.close();
	long t1 = System.currentTimeMillis();
	System.out.println("File reading took "+(float) ((t1-t0)*0.001)+" s");
 
	t0 = System.currentTimeMillis();
	// Calculate the Gravitational constant from the header constants
    G = gms[sunIndex] / Math.pow(Constant.SECONDS_PER_DAY, 2); // AU^3/s^2, assuming Msun = 1
 
    // Perform integration from closest date to final date
	speed = 450;
	double endDate = jdCalc;
	jd = jd0 + recordSkip * step;
	if (endDate < jd) speed = -Math.abs(speed);
 
	jd0 = jd;
	int iter = 0;
	while (Math.abs(endDate-jd) > Math.abs(speed/Constant.SECONDS_PER_DAY)) {
		Newton();
		iter ++;
		jd = jd0 + iter * speed / Constant.SECONDS_PER_DAY;
	}
	speed = (endDate - jd) * Constant.SECONDS_PER_DAY;
	if (speed != 0) Newton();
	t1 = System.currentTimeMillis();
	if (endDate != jd) System.out.println("*** ERROR *** Target end date "+endDate+
			" different from the integration final date "+jd);
	System.out.println("Integration took "+(float) ((t1-t0)*0.001)+" s");
	System.out.println();
	reportComparison();
}

And the output to the console:

File reading took 0.001 s
Integration took 0.29 s

- Comparison with JPL DE430 (ecliptic J2000 coordinates) -
JD (TDB): 2453175.0 (-103.500 days from starting time)
Sun
x-y-z (JPL):       0.003990643104  -0.001693103029  -0.000083882386
vx-vy-vz (JPL):    0.280646846E-10  0.808016691E-10  -0.013543813E-10
x-y-z (NBody):     0.003990643196  -0.001693103066  -0.000083882383
vx-vy-vz (NBody):  0.280646651E-10  0.808016788E-10  -0.013543821E-10
Mercury
x-y-z (JPL):       0.028259330560  0.304540976661  0.022706909935
vx-vy-vz (JPL):    -3898.933351E-10  377.518954E-10  388.607564E-10
x-y-z (NBody):     0.028259347951  0.304540970309  0.022706907832
vx-vy-vz (NBody):  -3898.933375E-10  377.519022E-10  388.607572E-10
Venus
x-y-z (JPL):       0.053438644149  -0.727094210348  -0.012860766505
vx-vy-vz (JPL):    2320.175833E-10  151.558747E-10  -131.853425E-10
x-y-z (NBody):     0.053438646912  -0.727094209264  -0.012860766639
vx-vy-vz (NBody):  2320.175835E-10  151.558755E-10  -131.853425E-10
Earth
x-y-z (JPL):       -0.039959748872  -1.016878246015  -0.000075683973
vx-vy-vz (JPL):    1958.000629E-10  -92.667774E-10  -0.052915E-10
x-y-z (NBody):     -0.039959747319  -1.016878245237  -0.000075683995
vx-vy-vz (NBody):  1958.000628E-10  -92.667761E-10  -0.052910E-10
Moon
x-y-z (JPL):       -0.040175145938  -1.014177761687  0.000121336287
vx-vy-vz (JPL):    1893.4082504E-10  -98.5029951E-10  3.2901527E-10
x-y-z (NBody):     -0.040175186712  -1.014177763797  0.000121339441
vx-vy-vz (NBody):  1893.4083363E-10  -98.5040310E-10  3.2897869E-10
Mars
x-y-z (JPL):       -1.142725939406  1.188491115884  0.053025037416
vx-vy-vz (JPL):    -1104.769477E-10  -984.942747E-10  6.479997E-10
x-y-z (NBody):     -1.142725938718  1.188491115345  0.053025037115
vx-vy-vz (NBody):  -1104.769479E-10  -984.942747E-10  6.479998E-10
Jupiter
x-y-z (JPL):       -5.376855052145  0.752951376040  0.117191692073
vx-vy-vz (JPL):    -131.936555E-10  -823.729550E-10  6.367927E-10
x-y-z (NBody):     -5.376855052275  0.752951376035  0.117191692069
vx-vy-vz (NBody):  -131.936554E-10  -823.729550E-10  6.367927E-10
Saturn
x-y-z (JPL):       -2.492288216405  8.688240427242  -0.051982283070
vx-vy-vz (JPL):    -655.136728E-10  -179.521257E-10  29.194241E-10
x-y-z (NBody):     -2.492288216410  8.688240427269  -0.051982283070
vx-vy-vz (NBody):  -655.136728E-10  -179.521257E-10  29.194241E-10
Uranus
x-y-z (JPL):       18.020641218188  -8.790799865008  -0.266144455503
vx-vy-vz (JPL):    196.2492898E-10  387.9291982E-10  -1.1006700E-10
x-y-z (NBody):     18.020641218192  -8.790799865011  -0.266144455503
vx-vy-vz (NBody):  196.2492898E-10  387.9291982E-10  -1.1006700E-10
Neptune
x-y-z (JPL):       20.755785980924  -21.764370998469  -0.030145641469
vx-vy-vz (JPL):    260.5476930E-10  252.8587288E-10  -11.2116914E-10
x-y-z (NBody):     20.755785980926  -21.764370998471  -0.030145641469
vx-vy-vz (NBody):  260.5476930E-10  252.8587288E-10  -11.2116914E-10
Pluto
x-y-z (JPL):       -4.821628395886  -30.077988684924  4.613220276365
vx-vy-vz (JPL):    365.7376005E-10  -121.5658970E-10  -92.7842624E-10
x-y-z (NBody):     -4.821628395886  -30.077988684926  4.613220276365
vx-vy-vz (NBody):  365.7376005E-10  -121.5658970E-10  -92.7842624E-10
Maximum difference (AU): 4.077392814355596E-8 (for Moon)
Maximum difference from Earth ("): 3.230750768894152 (for Moon)

In the previous example the calculation date is -103.5 days far from the closest record (close to a mean value, can be between 0 and 180), so the integration time (0.3s) is not that bad. In addition, as compensation all the planets are computed. But the main asteroids are not included and this has some minor effects in accuracy. To complete this method I would allow to use as initial conditions (when calling the method) the final position computed and the final date so minor effects like light-time corrections could be computed faster (without searching the file for a record and integration from a far date) with an ephemerides reduction program similar to the ones present in JPARSEC. However, with only the cartesian coordinates the main object of this section is achieved and you can use the algorithms described in my post about computing the position of the Sun, Moon, and planets to reduce these ecliptic J2000 coordinates and get the same output data from Earth returned in that post. As described in the section for planetary ephemerides a precession correction from J2000 to the equinox of date is necessary, after that the doCalc method can be used. As mentioned before the light-time correction would be a little integration from a previous initial conditions close in time. The accuracy would be highly superior to what is possible with the methods described in that post, but computation would be also slower. With the reduction method described in that post the accuracy is limited to 1”, if you are interested I can post the code.

Of course, integration times can be reduced and accuracy increased by writting two or more records per year, but the file would get bigger and more accuracy is hardly necessary. In case of computing only the position of an external planet the speed value can be increased to reduce the computing time. Respect the size, a 1.5 MB file could be used to compute planetary ephemerides with a typical error always well below 0.1” during 3 millenia, and would be easy to create for any JPL integration version. By contrast, the JPL files would require 3 GB of disk. In many cases the maximum accuracy or the libration angles are not needed so something like this can be useful.

Another practical use of this program would be to propagate orbits of other minor bodies or a mass-less particle like a meteor observed from different positions on Earth. With several photographs and by triangulation it is possible to derive the coordinates of the meteor when it burns (start and end times of the trail), and that information can be used to compute the position and velocity of the meteor, derive the orbit, and to identify the original comet where it came from.