# PSCB57H3 Lecture Notes - Lecture 8: Celestial Mechanics, Astronomical Object, Fiji

10 pages93 viewsSpring 2017

Department

Physical SciencesCourse Code

PSCB57H3Professor

Hanno ReinLecture

8This

**preview**shows pages 1-3. to view the full**10 pages of the document.**Introduction to Scientiﬁc Computing

Lecture 8

Professor Hanno Rein

Last updated: November 7, 2017

1 N-body integrations

1.1 Newton’s law

We’ll discuss today an important real world application of numerical ODE integration: celestial me-

chanics. Celestial mechanics is the branch of astronomy that deals with the motions of any kind of

celestial object. These motions are governed by Newton’s law of gravitation. In some cases there are

small correction due to general relativity, tides or radiation forces, which we will not include in our

model.

Newtonian gravity says that a bodies of mass mifeels another body of mass mj, separated by the

vector ~r, with the following force

~

Fij =−G mimj

|~rij |2·~rij

|~rij |

where Gis the gravitational constant. The last term is simply telling us in which direction the force

is directed (towards the other body). If there are Nbodies in total, then the force is a sum over all

bodies (except the body on which the force is exerted):

~

Fi=−

N−1

X

j=0

i6=j

G mimj

|~rij |2·~rij

|~rij |

To get an acceleration from a force, we use the well known law

~

F=m~a

so that we end up with the acceleration on the i-particle being

~ai=−

N−1

X

j=0

i6=j

G mj

|~rij |2·~rij

|~rij |

The acceleration is the second time derivate of the position, so we can also write the last equation as

¨

~ri=−

N−1

X

j=0

i6=j

G mj

|~rij |2·~rij

|~rij |

We have a total of Nparticles, thus we end up with a set of diﬀerential equation of second order

¨

~ri=F(~r)

1

###### You're Reading a Preview

Unlock to view full version

Only half of the first page are available for preview. Some parts have been intentionally blurred.

one for each i= 0...N −1 where F(~r) means that Fdepends on the rof all particles. All those

equations are coupled. We already discussed how to solve this equation using a trick to convert it so

ﬁrst order.

r

˙r0

=˙r

F(r)

Thus, we can simply use the Euler or midpoint method which we already know. Let’s try that.

2

###### You're Reading a Preview

Unlock to view full version

Only half of the first page are available for preview. Some parts have been intentionally blurred.

1.2 Initial conditions

As a test case, we use data from NASA’s Horizon system to simulate the Solar system. Besides

the major planets, we also include the comet 67P/Churyumov-Gerasimenko which is currently being

visited by the Rosetta spacecraft. The data is given in the ﬁle init.txt and has one line per body

and the format of name, followed by the x, y and z position and the x, y and z velocity. The units are

km for the position, km/s for the velocity and the mass is given in kg.

sun 1.988544 E+30 .4914033 47836458 E+05

−3.632715592552171 E+05 −1.049148558556447E +04

6.242516505718979E−03 1.163260794114081E−02

−2.475374674040771E−04

mercury 3.302 E+23 5.419423 38 5170247 E +07

−1.625487416441774 E+07 −6.231968220825911E +06

4.381057475129131E+00 4.891120140478765E+01

3.593097527852676E+00

venus 4.8685 E+24 −2.256798710024889E +07

1.046561273591759E+08 2.760223024328955E+06

−3.431401385797417 E+01 −7.710006447409605E +00

1.875067390915588E+00

earth 5.97219 E+24 −1.418287679667581E +08

4.126884716814923E+07 −1.125285256157373E+04

−8.830883522656004 E+00 −2.868395352996171E +01

−1.085239827735510E−05

mars 6.4185E +23 2.40827 9029085545 E+07

2.311290271817935E+08 4.261631465386531E+06

−2.317767379354534 E +01 4.5213 1275238314 1 E+00

6.638119814842468E−01

jupiter 1.89813 E+27 −7.778414202838075 E+08

2.244518699271340E+08 1.647441732604697E+07

−3.784015648894145 E+00 −1.193534649902395E +01

1.342303598962500E−01

saturn 5.68319 E+26 −2.818805400917871 E+08

1.321479657136385E+09 −1.177469869993168E+07

−9.963111372343775 E+00 −2.038111913779875E +00

4.325059188729626E−01

uranus 8.68103 E +25 2.668275 057715974 E+09

−1.368207197729832 E+09 −3.965456607834578E +07

3.057369253435637E+00 5.742539618944625E+00

−1.841246675277081E−02

neptune 1.0241 E+26 3.068964 88 3542690 E +09

−3.290508769241064 E+09 −2.965492893564129E +06

3.938602703572163E+00 3.739029327757723E+00

−1.671116647134993E−01

67P 0. −5.829716977332731 E+08 −3.371344262751943E

+08 3.008099774720960E+07 −1.302925387423112E

+00 −1.173853622210837 E+01 −7.989778520791224E

−01

Code 1: The ﬁle init.txt

The plot below show the position of the planet on March 4th 2004.

3

###### You're Reading a Preview

Unlock to view full version

#### Loved by over 2.2 million students

Over 90% improved by at least one letter grade.