Numerical model of the Solar System: Difference between revisions

Revision as of 03:17, 12 April 2005

Latest revision as of 07:06, 14 August 2025

(108 intermediate revisions by 75 users not shown)

Line 1:

{{Short description|Equations to predict the position of planets}}

==A Solar System Numerical Model==

{{No footnotes|article|date=April 2009}}A '''numerical model of the Solar System''' is a set of mathematical equations, which, when solved, give the approximate positions of the planets as a function of time. Attempts to create such a model established the more general field of [[celestial mechanics]]. The results of this simulation can be compared with past measurements to check for accuracy and then be used to predict future positions. Its main use therefore is in preparation of almanacs.

==Older efforts==

===A Law of Motion for Entities in a Solar System Numerical Model===

The simulations can be done in either [[Cartesian coordinate system|Cartesian]] or in [[Spherical coordinate system|spherical]] coordinates. The former are easier, but extremely calculation intensive, and only practical on an electronic computer. As such only the latter was used in former times. Strictly speaking, the latter was not much less calculation intensive, but it was possible to start with some simple approximations and then to add [[Perturbation (astronomy)|perturbations]], as much as needed to reach the wanted accuracy.

In essence this mathematical simulation of the [[Solar System]] is a form of the ''[[N-body problem]]''. The symbol '''''N''''' represents the number of bodies, which can grow quite large if one includes the Sun, 8 planets, dozens of moons, and countless planetoids, comets and so forth. However the influence of the Sun on any other body is so large, and the influence of all the other bodies on each other so small, that the problem can be reduced to the analytically solvable 2-body problem. The result for each planet is an orbit, a simple description of its position as function of time. Once this is solved the influences moons and planets have on each other are added as small corrections. These are small compared to a full planetary orbit. Some corrections might be still several degrees large, while measurements can be made to an accuracy of better than 1″.

The solar system may be modeled numerically using both [[Newton's Law of Gravitation]] and [[Newton's Second Law of Motion]].

Although this method is no longer used for simulations, it is still useful to find an approximate [[ephemeris]] as one can take the relatively simple main solution, perhaps add a few of the largest perturbations, and arrive without too much effort at the wanted planetary position. The disadvantage is that perturbation theory is very advanced mathematics.

Newton's Universal Law of Gravitation states that if a test [[mass]] m resides at a [[coordinate]] [[vector]] <math>\vec{r}</math>, then the [[gravitational]] [[force]] acting on the test mass, due to another mass that is present is:

==Modern method==

<math>\vec{F} = G \frac{m M_i}{|\vec{r}_i - \vec{r}|^3} (\vec{r}_i - \vec{r})</math>

The modern method consists of numerical integration in 3-dimensional space. One starts with a high accuracy value for the position (''x'', ''y'', ''z'') and the velocity (''v_x'', ''v_y'', ''v_z'') for each of the bodies involved. When also the mass of each body is known, the acceleration (''a_x'', ''a_y'', ''a_z'') can be calculated from [[Newton's law of gravitation]]. Each body attracts each other body, the total acceleration being the sum of all these attractions. Next one chooses a small time-step Δ''t'' and applies [[Newton's second law of motion]]. The acceleration multiplied with Δ''t'' gives a correction to the velocity. The velocity multiplied with Δ''t'' gives a correction to the position. This procedure is repeated for all other bodies.

The result is a new value for position and velocity for all bodies. Then, using these new values one starts over the whole calculation for the next time-step Δ''t''. Repeating this procedure often enough, and one ends up with a description of the positions of all bodies over time.

where:

The advantage of this method is that for a computer it is a very easy job to do, and it yields highly accurate results for all bodies at the same time, doing away with the complex and difficult procedures for determining perturbations. The disadvantage is that one must start with highly accurate figures in the first place, or the results will drift away from the reality in time; that one gets ''x'', ''y'', ''z'' positions which are often first to be transformed into more practical ecliptical or equatorial coordinates before they can be used; and that it is an all or nothing approach. If one wants to know the position of one planet on one particular time, then all other planets and all intermediate time-steps are to be calculated too.

:<math>\vec{F}</math> is the [[gravitational force]] acting on mass m.

==Integration==

:<math>G</math> is the [[gravitation constant]].

In the previous section it was assumed that acceleration remains constant over a small timestep Δt so that the calculation reduces to simply the addition of V × Δt to R and so forth. In reality this is not the case, except when one takes Δt so small that the number of steps to be taken would be prohibitively high. Because while at any time the position is changed by the acceleration, the value of the acceleration is determined by the instantaneous position. Evidently a full integration is needed.

Several methods are available. First notice the needed equations:

:<math>M_i</math> is the mass of an entity i. For purposes of modeling a [[solar system]], it is assumed to be a mass that is large enough to make a difference in the motion of test mass m. It could be a [[star]], a [[planet]], a [[moon]], a large [[asteroid]], or maybe even a large [[comet]]. A small [[rock (geology)|rock]] or a [[spacecraft]] would not be massive enough.

:<math>\vec{r}_i</math> is the coordinate vector of entity i.

In a solar system, there will be many entities i. Therefore, the total gravitational force acting on test mass m will be the sum of the forces for each i, over all n forces that are present:

<math>\vec{F} = \sum_i^n G \frac{m M_i}{|\vec{r}_i - \vec{r}|^3} (\vec{r}_i - \vec{r})</math>

Here we have changed the definition of <math>\vec{F}</math> to mean the total gravitational force due to all entities i.

<math>\vec{F}</math> is the total gravitational force acting on mass m.

[[Newton's Second Law of Motion]] states that if a net force of <math>\vec{F}</math> acts on a test mass m, that the test mass will undergo an [[acceleration]] <math>\vec{a}</math> such that all three quantities are related by:

<math>\vec{F} = m \vec{a}</math>

These two laws have the common left-hand-side <math>\vec{F}</math>. They may therefore be combined. When this is done and when the resulting equation is divided through by the common factor m, an equation of motion for test mass m results:

<math>\vec{a} = \sum_i^n G \frac{M_i}{|\vec{r}_i - \vec{r}|^3} (\vec{r}_i - \vec{r})</math>

At this point the concept of a test mass m is dispensed with. In order to model a solar system, each of the entities j of mass <math>M_j</math>, in turn, take the place of the test mass. The coordinate vector <math>\vec{r}</math> of the test mass therefore becomes the coordinate vector <math>\vec{r}_j</math> of entity j, and the summation is done over all other entities i except for the case where i equals j itself.

<math>\vec{a}_j = \sum_{i \neq j}^n G \frac{M_i}{|\vec{r}_i - \vec{r}_j|^3} (\vec{r}_i - \vec{r}_j)</math>

This equation describes the acceleration all bodies '''i''' running from 1 to N exercise on a particular body '''j'''. It is a vector equation, so it is to be split in 3 equations for each of the X, Y, Z components, yielding:

where obviously both i and j run from 1 to n.

The above equation is written in vector form. In order to have a set of formulas that is useful, the vector equation must be split into its three [[vector component]] equations.

<math>(a_j)_x = \sum_{i \neq j}^n G \frac{M_i}{d^3} (x_i - x_j)</math>

<math>(a_j)_y = \sum_{i \neq j}^n G \frac{M_i}{d^3} (y_i - y_j)</math>

<math>(a_j)_z = \sum_{i \neq j}^n G \frac{M_i}{d^3} (z_i - z_j)</math>

where

:<math>d = ( (x_i - x_j)^2 + (y_i - y_j)^2 + (z_i - z_j)^2 )^{1/2} </math>

At this point we have an [[equation of motion]] for all entities j in the solar system to be modeled.

===Reduced Units===

The [[mks]] system of [[units]] consists of three basic units for the three basic mechanical scalars. Length is in [[meters]] [m], [[mass]] is in [[kilograms]] [kg], and [[time]] is in [[seconds]] [s].

This would be a terribly awkward system of units for modeling a solar system. A far more convenient system of units would be: length in astronomical units [ [[AU]] ], mass in [[Terran]] masses [mTerran], and time in [[days]] [d].

An astronomical unit [AU] is the mean distance from the center of the [[Earth]] to the center of the [[Sun]], over the course of any given recent year. In this system of units, the Earth could be given an initial coordinate of, say, (x = 1.0, y = 0.0). A planet that is toward the outer edge of the solar system, such as [[Pluto]], could be given an initial coordinate of, say, (x = 40.0, y = 0.0). The magnitude of the distances that occur in a solar system model is quite manageable in this system of units.

A Terran mass [mTerran] is simply the mass of the Earth. In this system of units, the Earth would have a mass of (m = 1.0). A large planet, (the gas giant [[Jupiter]] being an example of such), could have a mass of, say (m = 300.0). A star, (our Sun being an example of such), could have a mass of, say, (m = 300,000.0) The magnitude of masses that occur in a solar system model is quite intuitive in this system of units. It is fairly easy to picture that another planet could have, say 10 times, or 1000 times the mass of the Earth, and that a star could have, say 1 million times the mass of the Earth.

A [[mean solar day]] of Earth [d] is a familiar unit that is easy to picture. It takes our Earth a [[year]] made up of 364.24 means solar days to [[orbit]] 360 [[degree (angle)|degree]]s around our Sun. If a planet at the outer edge of the solar system took, say, 200.0 years to orbit around the Sun, then this would be an interval of time that could be easily pictured, and easily converted to days. The magnitude of the time intervals that occur in a solar system model is both manageable and intuitive in this system of units.

For reference, the conversion between the MKS system, and our AU-mTerran-d system, is given below:

1.0 AU = 149.5985 x 10^9 m

1.0 mTerran = 5.977 x 10^24 kg

1.0 d = 60. * 60. * 24. s = 86,400.0 s

The Universal Constant of Gravitation, G is, in MKS units:

G = 6.673 x 10^-11 m^3 kg^-1 s^-2

If we are going to go ahead and use our AU-mTerran-d system of units, then we must convert this [[constant]] to that system of units and use the converted constant in our [[numerical model]] of the solar system. Let us call this converted constant 'g'.

g = 6.673 x 10^-11 (m^3 kg^-1 s^-2) (149.5985 x 10^9 m/AU)^-3 (5.977 x 10^24 kg/mTerran)^+1 (86,400.0 s/d)^+2

= (6.673 x 10^-11) (149.5985 x 10^9)^-3 (5.977 x 10^24)^+1 (86,400.0)^-2 AU^3 mTerran^-1 d^-2

= 1.595 861 x 10^-29 AU^3 mTerran^-1 d^-2

It is this converted constant g that will be used in the code. The converted constant g may appear to be a less convenient number than the original constant G. However, the more convenient value for G appears in one of the input files, and the awkward conversion to g is done automatically in the code.

The initial velocities of the solar system entities are given in units of [km/s]. However, users who have become used to our AU-mTerran-d system may wish to assign velocities in units of [AU/d]. The conversion factors are:

1.0 AU / d = 1731.464 120 km / s

1.0 km / s = 0.000 577 459 AU / d

The [[code]] expects the initial velocities to be in units of [km/s]. The code was written this way for the reason that users may wish to generalize it so that the [[trajectories]] of spacecraft through the solar system may be calculated. Spacecraft using conventional [[chemical]] [[rocket fuels]] function by applying an [[engine]] [[burn]] in order to attain a change in [[velocity]] (a [[delta-v]]) in any given direction. The distance over which the spacecraft travels while a given engine burn is applied, is negligible compared to the distance scale of the solar system through which it will travel. Therefore a given delta-v may be considered to occur [[instantaneously]]. A delta-v is typically in units of [km/s]. Therefore while distances are often given in [AU], and time intervals in [d], velocities are often given in [km/s].

Now that we have both a set of equations of motion and a convenient system of units, we next need an [[algorithm]] for computing the [[coordinates]] of the solar system entities, as time passes.

===A Time-Stepping Algorithm===

The next step is to come up with an algorithm in which the coordinates <math>\vec{r}_j</math> may be calculated as time passes.

We start with the definitions for velocity and acceleration:

<math>\vec{a} = \frac{d \vec{v}}{dt}</math>

<math>\vec{v} = \frac{d \vec{r}}{dt}</math>

If instead of an [[infinitesimal]] [[time step]] <math>dt</math> we consider a [[finite]] time step <math>\delta t</math>, then we obtain definitions for the [[average]] velocity and acceleration over the time interval <math>\delta t</math>

<math>\vec{a}_{ave} = \frac{\Delta \vec{v}}{\Delta t}</math>

<math>\vec{v}_{ave} = \frac{\Delta \vec{r}}{\Delta t}</math>

or

<math>\Delta \vec{v} = \vec{a}_{ave} \Delta t </math>

<math>\Delta \vec{x} = \vec{v}_{ave} \Delta t </math>

Splitting these vector equations into the vector component equations:

<math>(a_{ave})_x = \frac{\Delta v_x}{\Delta t}</math>

<math>(a_{ave})_y = \frac{\Delta v_y}{\Delta t}</math>

<math>(a_{ave})_z = \frac{\Delta v_z}{\Delta t}</math>

<math>(v_{ave})_x = \frac{\Delta x}{\Delta t}</math>

<math>(v_{ave})_y = \frac{\Delta y}{\Delta t}</math>

<math>(v_{ave})_z = \frac{\Delta z}{\Delta t}</math>

or:

<math>\Delta v_x = (a_{ave})_x \Delta t </math>

<math>\Delta v_y = (a_{ave})_y \Delta t </math>

<math>\Delta v_z = (a_{ave})_z \Delta t </math>

<math>\Delta x = (v_{ave})_x \Delta t </math>

<math>\Delta y = (v_{ave})_y \Delta t </math>

<math>\Delta z = (v_{ave})_z \Delta t </math>

Next we need an algorithm for stepping the coordinates (x, y, z)_j through time. Let us consider just coordinate x for just one entity, and let us consider the following algorithm:

input t_f, (delta)t (We need input values for the time step size (delta)t

and the final time t_f)

t = 0. (We start the time at zero)

x = x_0 (We set the coordinate x to its initial value)

v = v_0 (We set the velocity v to its initial value)

a = f(x) (The acceleration is a function of x and we calculate it here) (3.)

output t, x

loop until (t = t_f)

| t = t + (delta)t (We advance the time)

| x = x + v * (delta)t (We advance the coordinate) (1.)

| v = v + a * (delta)t (We advance the velocity) (2.)

| a = f(x) (We calculate the new acceleration as a function of x) (3.)

| output t, x

This algorithm is called the [[Euler]] Algorithm. It is the most basic [[time stepping algorithm]] that can be used to solve a [[point particle]] [[classical mechanics]] problem. The Euler Algorithm, however, has problems. The formulas (1.) and (2.) are based on the definitions given earlier of [[average velocity]] and [[average acceleration]] over a finite time interval, and yet it uses velocitites and accelerations that occur at the beginning of the time interval.

A better algorithm is to use velocities and accelerations that occur in the middle of the time intervals. This is quite simply accomplished. We start by advancing the velocity by only a half step rather than a full step. From that point on, we advance the coordinate and the velocity by a full step, except that in this algorithm, the velocity that is used to advance the coordinate, occurs in the middle of the coordinate time step, and the acceleration that is used to advance the velocity, occurs in the middle of the velocity time step.

input t_f, (delta)t (We need input values for the time step size (delta)t

and the final time t_f)

t = 0. (We start the time at zero)

x = x_0 (We set the coordinate x to its initial value)

a = f(x) (The acceleration is a function of x and we calculate it here) (3.)

v = v_0 (We set the velocity v to its initial value)

output t, x

v = v + a * (delta)t/2. (We advance the velocity by half a time step)

loop until (t = t_f)

| t = t + (delta)t (We advance the time)

| x = x + v * (delta)t (We advance the coordinate by a full time step) (1.)

| a = f(x) (We calculate the new acceleration as a function of x) (3.)

| v = v + a * (delta)t (We advance the velocity by a full time step) (2.)

| output t, x

This algorithm is called the [[Feynman]] Algorithm, or the [[Verlet]] Algorithm.

We now have all the tools we need to construct a [[numerical model]] of the solar system. We have an equation of motion for all the solar system entities j, which provides us with the acceleration of all the solar system entities j at any given instant. We also have a time stepping algorithm which uses this acceleration information for a given entity to step its coordinates and velocity through time.

===A Working Code===

The following [[code]] constitutes such a [[numerical model]]. The code is written in [[fortran 77]]. It reads an [[input]] file]] containing start-up information for a solar system, and writes an [[output]] file that contains the coordinates of each solar system entity, at each time step in the simulation, from an initial time of zero, up until the final time.

program solarV2f

implicit none

c

call readConsts()

call readSystem()

call moveSystem()

c

stop

end

c

subroutine readConsts()

implicit none

c

double precision rAu, mTerran, g, deltaT, tF

double precision tDay

common /fundConsts/ rAu, mTerran, g, deltaT, tF, tDay

c

integer nin, nout, nlog

parameter(nin=1, nout=2, nlog=3)

c

open(nin, file='const.dat', status='old')

read(nin, *) rAu

read(nin, *) mTerran

read(nin, *) g

read(nin, *) deltaT

read(nin, *) tF

close(nin)

c

tDay = 60. * 60. * 24.

g = g * mTerran * tDay * tDay / (rAu * rAu * rAu)

c

return

end

c

subroutine readSystem()

implicit none

c

double precision PI

parameter(PI=3.141592)

c

double precision rAu, mTerran, g, deltaT, tF

double precision tDay

common /fundConsts/ rAu, mTerran, g, deltaT, tF, tDay

c

integer nMax, n

parameter(nMax = 20)

character*20 name(nMax)

double precision mass(nMax)

double precision rx(nMax), ry(nMax), rz(nMax)

double precision vOldx(nMax), vOldy(nMax), vOldz(nMax)

double precision vAvex(nMax), vAvey(nMax), vAvez(nMax)

common /size/ n

common /system/ name, mass, rx, ry, rz,

& vOldx, vOldy, vOldz, vAvex, vAvey, vAvez

c

double precision r, theta, phi

double precision vr, vtheta, vphi

double precision vx(nMax), vy(nMax), vz(nMax)

integer i

c

integer nin, nout, nlog

parameter(nin=1, nout=2, nlog=3)

c

open(nin, file='solar.dat', status='old')

read(nin, *) n

do i = 1, n

read(nin, 10) name(i)

read(nin, *) mass(i)

read(nin, *) r, theta, phi

theta = theta * (PI/180.)

phi = phi * (PI/180.)

rx(i) = r * sin(theta) * cos(phi)

ry(i) = r * sin(theta) * sin(phi)

rz(i) = r * cos(theta)

c

read(nin, *) vr, vtheta, vphi

vx(i) = vr * sin(theta) * cos(phi)

& + vtheta * cos(theta) * cos(phi)

& - vphi * sin(phi)

vy(i) = vr * sin(theta) * sin(phi)

& + vtheta * cos(theta) * sin(phi)

& + vphi * cos(phi)

vz(i) = vr * cos(theta)

& - vtheta * sin(theta)

c

vx(i) = (tDay * 1000. / rAu) * vx(i)

vy(i) = (tDay * 1000. / rAu) * vy(i)

vz(i) = (tDay * 1000. / rAu) * vz(i)

c

vOldx(i) = vx(i)

vOldy(i) = vy(i)

vOldz(i) = vz(i)

c

vAvex(i) = vx(i)

vAvey(i) = vy(i)

vAvez(i) = vz(i)

c

enddo

close(nin)

return

10 format(a20)

end

c

subroutine moveSystem

implicit none

c

double precision rAu, mTerran, g, deltaT, tF

double precision tDay

common /fundConsts/ rAu, mTerran, g, deltaT, tF, tDay

c

integer nMax, n

parameter(nMax = 20)

character*20 name(nMax)

double precision mass(nMax)

double precision rx(nMax), ry(nMax), rz(nMax)

double precision vOldx(nMax), vOldy(nMax), vOldz(nMax)

double precision vAvex(nMax), vAvey(nMax), vAvez(nMax)

common /size/ n

common /system/ name, mass, rx, ry, rz,

& vOldx, vOldy, vOldz, vAvex, vAvey, vAvez

c

integer i

double precision t

double precision vNewx, vNewy, vNewz

double precision ax, ay, az

c

integer nin, nout, nlog

parameter(nin=1, nout=2, nlog=3)

c

open(nout, file='coords.dat', status='unknown')

c

c Start at time zero

c

t = 0.

call writeSystem(t)

c

c Perform initial half-step

c

do i = 1, n

call accelMass(i, g, ax, ay, az)

vOldx(i) = vOldx(i) + (deltaT * ax) / 2.

vOldy(i) = vOldy(i) + (deltaT * ay) / 2.

vOldz(i) = vOldz(i) + (deltaT * az) / 2.

enddo

c

c Perform time steps

c

do while (t < tF)

t = t + deltaT

do i = 1, n

rx(i) = rx(i) + (deltaT * vOldx(i))

ry(i) = ry(i) + (deltaT * vOldy(i))

rz(i) = rz(i) + (deltaT * vOldz(i))

enddo

do i = 1, n

call accelMass(i, g, ax, ay, az)

vNewx = vOldx(i) + (deltaT * ax)

vNewy = vOldy(i) + (deltaT * ay)

vNewz = vOldz(i) + (deltaT * az)

vAvex(i) = (vOldx(i) + vNewx) / 2.

vAvey(i) = (vOldy(i) + vNewy) / 2.

vAvez(i) = (vOldz(i) + vNewz) / 2.

vOldx(i) = vNewx

vOldy(i) = vNewy

vOldz(i) = vNewz

enddo

call writeSystem(t)

enddo

c

c Calculations are done

c

close(nout)

return

end

c

subroutine accelMass(i, g, ax, ay, az)

implicit none

c

integer i

double precision g

double precision ax, ay, az

c

integer nMax, n

parameter(nMax = 20)

character*20 name(nMax)

double precision mass(nMax)

double precision rx(nMax), ry(nMax), rz(nMax)

double precision vOldx(nMax), vOldy(nMax), vOldz(nMax)

double precision vAvex(nMax), vAvey(nMax), vAvez(nMax)

common /size/ n

common /system/ name, mass, rx, ry, rz,

& vOldx, vOldy, vOldz, vAvex, vAvey, vAvez

c

integer j

double precision rsx, rsy, rsz

double precision rox, roy, roz

double precision rosx, rosy, rosz

double precision ms

double precision r2, r, r3

c

c Initialize the acceleration

c

ax = 0.

ay = 0.

az = 0.

c

c Set the coordinates of the observer mass

c

rox = rx(i)

roy = ry(i)

roz = rz(i)

c

c Loop over all other masses

c

do j = 1, n

if (j .ne. i) then

ms = mass(j)

c

rsx = rx(j)

rsy = ry(j)

rsz = rz(j)

c

rosx = rsx - rox

rosy = rsy - roy

rosz = rsz - roz

c

r2 = rosx*rosx + rosy*rosy + rosz*rosz

r = dsqrt(r2)

r3 = r * r * r

c

ax = ax + g * ms * rosx / r3

ay = ay + g * ms * rosy / r3

az = az + g * ms * rosz / r3

endif

enddo

c

return

end

c

subroutine writeSystem(t)

implicit none

c

double precision t

c

integer nMax, n

parameter(nMax = 20)

character*20 name(nMax)

double precision mass(nMax)

double precision rx(nMax), ry(nMax), rz(nMax)

double precision vOldx(nMax), vOldy(nMax), vOldz(nMax)

double precision vAvex(nMax), vAvey(nMax), vAvez(nMax)

common /size/ n

common /system/ name, mass, rx, ry, rz,

& vOldx, vOldy, vOldz, vAvex, vAvey, vAvez

c

integer i

c

integer nin, nout, nlog

parameter(nin=1, nout=2, nlog=3)

c

c Write the time

c

write(nout,200) t

c

c Write the number of masses in this time frame

c

write(nout,100) n

c

c Loop over all the masses and write the coordinates of each one

c

do i = 1, n

write(nout,300) rx(i), ry(i), rz(i)

enddo

c

return

100 format(I4)

200 format(F12.6)

300 format(F12.6,4X,F12.6,4X,F12.6)

end

c

<math>(a_j)_x = \sum_{i \neq j}^n G \frac{M_i}{( (x_i - x_j)^2 + (y_i - y_j)^2 + (z_i - z_j)^2 )^{3/2}} (x_i - x_j)</math>

with the additional relationships

The following file contains the [[physical constants]] that the numerical model requires, as well as the time information.

<math>a_{x} = \frac{dv_{x}}{dt}</math>, <math>v_{x} = \frac{dx}{dt}</math>

likewise for Y and Z.

The former equation (gravitation) may look foreboding, but its calculation is no problem. The latter equations (motion laws) seem simpler, but yet they cannot be calculated. Computers cannot integrate, they cannot work with infinitesimal values, so instead of dt we use Δt and bringing the resulting variable to the left:

<math>\Delta v_x = a_{x} \Delta t </math>, and: <math>\Delta x = v_{x} \Delta t </math>

149.5985e09

5.977e24

6.673e-11

1.

1000.

Remember that '''a''' is still a function of time. The simplest way to solve these is just the [[Euler]] algorithm, which in essence is the linear addition described above. Limiting ourselves to 1 dimension only in some general computer language:

The file is in the form of:

a.old = gravitationfunction(x.old)

x.new = x.old + v.old * dt

v.new = v.old + a.old * dt

As in essence the acceleration used for the whole duration of the timestep, is the one as it was in the beginning of the timestep, this simple method has no high accuracy. Much better results are achieved by taking a mean acceleration, the average between the beginning value and the expected (unperturbed) end value:

[1.0 AU, in m]

[1.0 mTerran, in kg]

[G, in m^3 kg^-1 s^-2]

[Time step, in d]

[Total time interval of simulation, in d]

a.old = gravitationfunction(x.old)

x.expect = x.old + v.old * dt

a.expect = gravitationfunction(x.expect)

v.new = v.old + (a.old + a.expect) * 0.5 * dt

x.new = x.old + (v.new + v.old) * 0.5 * dt

Of course still better results can be expected by taking intermediate values. This is what happens when using the [[Runge-Kutta]] method, especially the one of grade 4 or 5 are most useful. The most common method used is the [[leapfrog method]] due to its good long term energy conservation.

The next file contains contains the information for a simplified solar system.

A completely different method is the use of [[Taylor series]]. In that case we write: <math>r = r_0 + r'_0 t + r''_0 \frac{t^2}{2!} + ... </math>

but rather than developing up to some higher derivative in r only, one can develop in r and v (that is r') by writing <math>r = f r_0 + g r'_0</math>and then write out the factors ''f'' and ''g'' in a series.

==Approximations==

To calculate the accelerations the gravitational attraction of each body on each other body is to be taken into account. As a consequence the amount of calculation in the simulation goes up with the square of the number of bodies: Doubling the number of bodies increases the work with a factor four. To increase the accuracy of the simulation not only more decimals are to be taken but also smaller timesteps, again quickly increasing the amount of work. Evidently tricks are to be applied to reduce the amount of work. Some of these tricks are given here.

By far the most important trick is the use of a proper integration method, as already outlined above.

3

Sun

332943

0.031254 90. 180.

0. 0. 0.079

Jupiter

2000

5.203 88.694 0.

0. 0. 13.07

Earth

1.

1. 90. 0.

0. 0. 29.78

The choice of units is important. Rather than to work in [[SI units]], which would make some values extremely small and some extremely large, all units are to be scaled such that they are in the neighbourhood of 1. For example, for distances in the Solar System the [[astronomical unit]] is most straightforward. If this is not done one is almost certain to see a simulation abandoned in the middle of a calculation on a [[floating point]] [[arithmetic overflow|overflow]] or [[arithmetic underflow|underflow]], and if not that bad, still accuracy is likely to get lost due to [[truncation]] errors.

The file is in the form of:

If N is large (not so much in Solar System simulations, but more in galaxy simulations) it is customary to create dynamic groups of bodies. All bodies in a particular direction and on large distance from the reference body, which is being calculated at that moment, are taken together and their gravitational attraction is averaged over the whole group.

[Number of solar system entities]

[Name of solar system entity 1]

[Mass in mTerran]

[Initial Radius from center, in AU] [Initial Theta, in degrees] [Initial Phi, in degrees]

[Initial Radius-direction velocity, in km/s] [Initial Theta-direction velocity, in km/s] [Initial Phi-direction velocity, in km/s]

The total amount of [[energy]] and [[angular momentum]] of a closed system are conserved quantities. By calculating these amounts after every time step the simulation can be programmed to increase the stepsize Δt if they do not change significantly, and to reduce it if they start to do so. Combining the bodies in groups as in the previous and apply larger and thus less timesteps on the faraway bodies than on the closer ones, is also possible.

To allow for an excessively rapid change of the acceleration when a particular body is close to the reference body, it is customary to introduce a small parameter ''e'' so that

The time info in the first file, and the solar system info in the second file, may be modified.

<math>a = \frac{G M}{r^2 + e}</math>

==Complications==

This article is intended to explain the [[physical theory]] that is required to numerically [[simulate]] solar system [[dynamics]] using [[Newton's Law of Gravitation]], and [[Newton's Second Law of Motion]], and to provide a simple code that serves as an example of how the [[theory]] may be coded. No attempt has been made here to graphically [[display]] or [[animate]] the [[motion]] of the solar system entities.

If the highest possible accuracy is needed, the calculations become much more complex. In the case of comets, nongravitational forces, such as radiation pressure and gas drag, must be taken into account. In the case of Mercury, and other planets for long term calculations, relativistic effects cannot be ignored. Then also the total energy is no longer a constant (because the four vector energy with linear momentum is). The finite speed of light also makes it important to allow for light-time effects, both classical and relativistic. Planets can no longer be considered as particles, but their shape and density must also be considered. For example, the flattening of the Earth causes precession, which causes the axial tilt to change, which affects the long-term movements of all planets.

Long term models, going beyond a few tens of millions of years, are not possible due to the lack of [[stability of the Solar System]].

==See also==

Users may wish to develop more elaborate codes that will do this. They may wish to use the code provided here as the 'physics engine' in such a code.

*[[Ephemeris]]

*[[VSOP (planets)]]

==References==

Readers may note that the code performs some complicated coordinate system conversion calculations, particularly in the calculation of the initial velocity components. These calculations are described in a separate article, on a [[transformation from spherical coordinates to rectangular coordinates]].

*{{Cite book|first=Dan L. |last=Boulet |title=Methods of orbit determination for the microcomputer |publisher=Willmann-Bell, Inc |___location=[[Richmond, Virginia]] |year=1991 |isbn=978-0-943396-34-7 |oclc=23287041}}{{Page needed|date=September 2010}}

{{Solar System models}}

{{DEFAULTSORT:Numerical Model Of Solar System}}

[[Category:Mathematics]]

[[Category:Numerical analysis]]

[[Category:Computational physics]]

[[Category:Dynamical systems]]

[[Category:Dynamics of the Solar System]]

[[Category:Solar System models]]