/* A simulation of the n-body problem.  The bodies here are the planets of
   the solar system and the Sun, with their starting positions and velocities
   as of July 1, 1998, midnight.
   The units are Earth masses (= 5.9736e24 kg), astronomical units 
   (=1.49597870691e11 m) and days (=86400s).
   */

/* Some parameters of the interpreter are set here.  For example, we want 
   the minimum interval between successive samples to be 10 days (if we run
   the simulation over 10 years).  Also we want to use the Bulirsch-Stoer 
   integration method, as it works well for smooth equations */
#SAMPLE_STYLE 0
#MEMORY_CHUNK 102400
#TRAIL_CHUNK 20000
#SAMPLE_INTERVAL_MIN 10
#INTEGRATOR BS
#INTEGRATION_INIT 10.0
#define g    8.88769408e-10 
    
/* defines the dynamics of planets.  The accelerations on each axis are
   computed by Newton's laws.  The positions are sampled. */
always Planet = (m, initpx, initpy, initpz, initvx, initvy, initvz)[px, py, pz, mass]{
    realvar(px, py, pz, mass),
    px = initpx, py = initpy, pz = initpz,
    px' = initvx, py' = initvy, pz'=initvz,
    sample(py),sample(px),sample(pz),
    always {
	cont(px), cont(py), cont(pz),
	mass := m,
	px'' := sum(g * P.mass * (P.px - px)/((P.px - px)^2 + (P.py -py)^2 + (P.pz -pz)^2)^1.5, Planet(P), P != Self),
	py'' := sum(g * P.mass * (P.py - py)/((P.px - px)^2 + (P.py -py)^2 + (P.pz -pz)^2)^1.5, Planet(P), P != Self),
	pz'' := sum(g * P.mass * (P.pz - pz)/((P.px - px)^2 + (P.py -py)^2 + (P.pz -pz)^2)^1.5, Planet(P), P != Self)

    }
},

Planet(Sun, 332981.78652,
      -0.008348511782195148, 0.001967945668637627, 0.0002142251001467145,
      -0.000001148114436325131, -0.000008994958827348018, 0.00000006538635311283777),

Planet(Mercury, 0.0552765501,
      -0.4019000379850893, -0.04633361689674035, 0.03239207992750900,
      -0.002423875040887606, -0.02672168963230259, -0.001959654820981497),

Planet(Venus, 0.8150026784,
      0.6680247657009936, 0.2606201175567890, -0.03529355196193388,
      -0.007293563117650372, 0.01879420958390879, 0.0006778739390714113),

Planet(Earth, 1.01230246,
      0.1508758612501242, -1.002162526305211, 0.0002082851504420832,
      0.01671098890724774, 0.002627047365383169, -0.0000004771611907632339),

Planet(Mars, 0.10749598232,
      0.3095796131231540, 01.509103334193779, 0.02397249271102571,
      -0.01316091582511848, 0.004066946366453539, 0.0004089587638718586),

Planet(Jupiter, 317.89790076,
      04.824637433017612, -1.196812497464260, -0.1030599072585649,
      0.001716726975735163, 0.007676089019438129, -0.00007021734714029865),

Planet(Saturn, 95.18552281,
      08.348813878243238, 04.101757994088924, -0.4034245819779162,
      -0.002753597799584471, 0.004996451166988408, 0.00002272558599524301),

Planet(Uranus, 14.53748693,
      012.87118968572842, -15.12106483933579, -0.2229438226188330,
      0.002965946969806189, 0.002366610226063414, -0.00002962300113861107),

Planet(Neptune, 17.16975693,
      015.35793700864867, -25.92547713196480, 0.1799517749484475,
      0.002680784056046587, 0.001618204887705156, -0.00009497127021434280),

Planet(Pluto, 0.00244408732,
      -11.53274193833463, -27.06743920236740, 06.232330236968182,
      0.002977569298809042, -0.001704821621937649, -0.0006820673428000879)