Here is a tutorial explaining how to model a two dimensional 2-body planetary system in Excel. It uses the Euler method of integration.

The tutorial starts with explaining the simple Newtonian laws acting on the two planets.

There are essentially just two forces acting on each body at any time: the inertia and the gravitational attraction. During each small time step, from the distance between bodies we can calculate the gravitational forces and from the Newton’s second law we can further calculate the acceleration. From acceleration and the initial speed (the speed at the beginning of the time step) we can calculate the speed at the end of each time step and from all the above known quantities we can find the position at the end of the time interval.

This iterative process is (dynamically) repeated over and over and the solution is plotted on a 2D scatter chart. A detailed step-by-step Excel-VBA implementation is further shown.

Recommending this site to your friends would be highly appreciated. Thanks for your support!

### Creating a basic planetary system in Excel

by George Lungu

* For the model presented here it is recommended that you use Excel 2003 or older. The 2007 will be either excessively slow.

We want to create the following model:

### Some basic theory from high school:

u=r/|r|

F = m*a – Newton’s 2nd law

F = u *k*m1*m2 / r^2 – Newton’s universal attraction formula

That was not a very rigorous. Pay attention to a few details:

1.- The forces acting on each body have the same direction and absolute value but opposite sense

2.- We are talking attraction not repulsion!

3.- The acceleration has the same direction and sense with the force or gravity.

4. Of course because this problem is in 2-dimenssional we need to solve it separately along x and y axes

### 2D vector decomposition

F2=-F1

F1x = F1*cos(a)

F1y = F1*sin(a)

x1 x2

F1 = k*m1*m2 / ((y2-y1)^2+(x2-x1)^2)

cos(a) = (x2-x1) / sqrt((y2-y1)^2+(x2-x1)^2)

sin(a) = (y2-y1) / sqrt((y2-y1)^2+(x2-x1)^2)

From Newton’s second law (F = m*a) and the decomposition formulas we can therefore write the x and y components of the accelerations of both bodies as:

a2=-a1

a1x = k*m2* (x2-x1) / ((y2-y1)^2+(x2-x1)^2)^(3/2)

a1y = k*m2* (y2-y1) / (( y2-y1)^2+(x2-x1)^2)^(3/2)

a2x = k*m1* (x1-x2) / ((y2-y1)^2+(x2-x1)^2 )^(3/2)

a2y = k*m1* (y1-y2) / ((y2-y1)^2+(x2-x1)^2 )^(3/2)

We also know the definition of acceleration as the derivative of speed with respect to time:

ax = dvx/dt or ay = dvy/dt

for a very small time increment we can assume that the acceleration is constant and re-write :

ax_current =(vx_current – vx_previous)/dt

ay_current =(vy_current – vy_previous)/dt

From here we can derive the following:

vx_current = vx_previous + ax_current *dt

vy_current = vy_previous + ay_current *dt <- Very important

Similarily:

vx = dx/dt or vy = dy/dt – for a very small time increment we can assume that the speed is constant and re-write :

vx-current =(xcurrent – xprevious)/dt

vy-current =(ycurrent – yprevious)/dt

From here we can derive the following: xcurrent = xprevious + vx-current *dt <- Very important

ycurrent = yprevious + vy-current *dt

### Outline:

a1x = k*m2* (x2-x1) / ((y2-y1)^2+(x2-x1)^2)^(3/2)

a1y = k*m2* (y2-y1) / ((y2-y1)^2+(x2-x1)^2 )^(3/2)

a1x = k*m1* (x1-x2) / ((y2-y1)^2+(x2-x1)^2 )^(3/2)

a1y = k*m1* (y1-y2) / ((y2-y1)^2+(x2-x1)^2 )^(3/2)

vx_current = vx_previous + ax_current *dt

vy_current = vy_previous + ay_current *dt

xcurrent = xprevious + vx-current *dt

ycurrent = yprevious + vy-current *dt

The above formulas are all we need to build a planetary system.

### So how is the differential equation system solved?

In the previous page you can see 3 different categories of formulas:

– Acceleration/Force calculations -> yellow formulas (any time: a=F/m or F=ma)

– Speed calculation -> green formulas

– Coordinate calculations -> blue formulas

### A pure spread sheet type solution:

The time in this case time will increase down the column. Each row will represent a discrete moment of time.

#### Advantages of the pure spreadsheet solution:

1. Requires practically no VBA (Visual Basic for Applications).

2. It is very fast

#### There are 2 disadvantage of the tabular solution:

1. The number of time steps is limited

2. The files are large

As we mentioned before, time advances vertically along the column.

The arrows and their colors relate to each type of equation:

#### Initial Conditions:

t[0] Speeds Coordinates

(Constants) ——>

time step – dt

t[1]

Forces Speeds Coordinates

t[2] Forces Speeds Coordinates

t[3] etc etc etc

#### Some important points:

The forces at time t[n] are calculated from the coordinates at time t[n-1] -> this is a compromise which gives decent solutions as long as dt is small enough (see the yellow formulas).

The speeds are calculated from the previous speeds (at t[n-1]) and the current accelerations (see the green formulas).

The coordinates are calculated from the previous coordinates and the current speeds (see the blue formulas).

### A sequential type solution:

There is just one row of calculations. A basic VBA macro will create an infinite loop.

#### Advantages of the pure spreadsheet solution:

1. It can run forever (no fixed/limited number of steps). It can be stopped and restarted.

2. The Excel files are small.

The only disadvantage of the sequential solution is that it is slower than the previous (instantaneous) version.

In our case this does not matter since the simulation will be display limited ( around 20- 40 frames per second depending on the computer used and the chart size).

The macro will copy the speed and coordinate values from the active row (t[n]) to the following row.

The arrows and their colors relate to each type of equation:

t[n]

Forces Speeds Coordinates

time step – dt

t[n-1] (Constants) ——> Speeds Coordinates

The speeds and coordinates at t[n-1] are copied from the speeds and coordinates at t[n]. Therefore the formulas are using their old results to calculate the results in an infinite loop.

### So, how is this whole thing implemented in Excel?

The following slides are cookbook-style.

There is no need to understand the previous sections. Just follow the recipe!

Create the input data area

Cell “A1”: Initial Conditions

Cell “A2”: x1 Cell “A3”: y1

Cell “A4”: x2 Cell “A5”: y2

Cell “A6”: v1x Cell “A7”: v1y

Cell “A8”: v2x Cell “A9”: v2y

Cell “A10”: k

Cell “A11”: dt

Cell “A13”: Zoom

Cell “A14”: m1

Cell “A15”: m2

Cell “A20”: Index

Cell “A21”: Trail Size

Input data area – figures I use font size 16 / bold

The only formula is in the cell “B20”: = L27/B11-1

For the rest of the cells just plug in those numbers as a starter.

Create spin button “k” – This button will change the gravitational constant on-the-fly.

View -> Toolbars -> Control Toolbox -> Spin Button (drag create)

In the control toolbox hit “Design Mode” then right click the button and select “Properties” Change the following:

(Name): k

Min: 0

Max: 200

#### Macro associated to button “k”:

Sub k_Change()

Range(“B10”) = k.Value

End Sub

#### Create spin button “Zoom”:

– This button will change the Zoom on-the-fly

View -> Toolbars -> Control

Toolbox -> Spin Button (drag create)

Right click the button and select “Properties”

Change the following:

(Name): Zoom

Min: 1

Max: 50

#### Macro associated to button “Zoom”:

Private Sub Zoom_Change()

Range(“B13”) = Zoom.Value ^ (1 + Zoom.Value / 100)

With ActiveSheet.ChartObjects(“Chart 1”).Chart

.Axes(xlCategory).MinimumScale = – Range(“B13”)

.Axes(xlCategory).MaximumScale = Range(“B13”)

.Axes(xlValue).MinimumScale = – Range(“B13”)

.Axes(xlValue).MaximumScale = Range(“B13”)

End With

Range(“C17”).Select

End Sub

#### Create spin button:

“Trail_Size”

This button will change the trail size on-the-fly:

View -> Toolbars -> Control

Toolbox -> Spin Button (drag

create)

Right click the button and

select “Properties”

Change the following:

(Name): Trail_Size

Min: 0

Max: 500

#### Macro associated to button “Trail_Size”:

Private Sub Trail_Size_Change()

Range(“B21”) = 10 * Trail_Size.Value

Range(Range(“D29”).Offset(Range(“B21”), 0), _

“L5000”).ClearContents

End Sub

Just a line break in VBA (the line was split in 2 by using “ _”

#### Create button “Start_Pause”:

This button will start or pause the simulation.

On Draw menu: Auto Shapes -> Basic Shapes -> Rounded.

Rectangle (drag)

On Draw menu: Text box -> Left Click in the Rectangle -> Type “Start-Pause”

Right Click the rectangle -> Assign macro -> Choose “Start-Pause” from menu.

##### Macro associated to button “Start_Pause”:

Public s As Boolean

————————————————–

Sub Start_Pause()

Just a line break in VBA (the

If s = False Then

line was split in 2 by using “ _”

s = True

Do

DoEvents

If s = False Then Exit Do

Range(“D28”, Range(“L28”).Offset(Range(“B21”), 0)) = _

Range(“D27”, Range(“L27”).Offset(Range(“B21”), 0)).Value

Loop

Else

s = False

Exit Sub

End If

End Sub

#### Observation

The original line of code didn’t give a trail:

Range(“D28:L28”)= Range(“D27:L27”).Value

Was replaced with:

Range(“D28”, Range(“L28”).Offset(Range(“B21”), 0)) = _

Range(“D27”, Range(“L27”).Offset(Range(“B21”), 0)).Value

Which creates a trail with programmable length behind each planet

#### Create button “Reset”:

This button will reset the simulation:

On Draw menu: Auto Shapes ->

Basic Shapes -> Rounded

Rectangle (drag)

On Draw menu: Text box ->

Left Click in the Rectangle ->

Type “Reset”

Right Click the rectangle ->

Assign macro -> Choose

“Reset_” from menu

#### Macro associated to button “Reset”:

Sub Reset_()

Range(“D28:L5000”).ClearContents

Range(“D28”) = Range(“B6”).Value

Range(“E28”) = Range(“B7”).Value

Range(“F28”) = Range(“B8”).Value

Range(“G28”) = Range(“B9”).Value

Range(“H28”) = Range(“B2”).Value

Range(“I28”) = Range(“B3”).Value

Range(“J28”) = Range(“B4”).Value

Range(“K28”) = Range(“B5”).Value

s = False

End Sub

#### The calculation area:

Cell “B26”: a1x Cell “C26”: a1y

Cell “A27”: t[n]

Cell “D26”: v1x Cell “E26”: v1y

Cell “A28”: t[n-1]

Cell “F26”: v2x Cell “G26”: v2y

Cell “B25”: accelerations

Cell “H26”: x1 Cell “I26”: y1

Cell “E25”: speeds

Cell “J26”: x2 Cell “K26”: y2

Cell “I25”: coordinates

Cell “L26”: time

Just names up to this point !

#### The calculation area – Active Formulas:

Cell “B27”: =B10*B15*B14*(J28-H28)/((J28-H28)^2+(K28-I28)^2)^(3/2)

Cell “C27”: =B10*B15*B14*(K28-I28)/((J28-H28)^2+(K28-I28)^2)^(3/2)

Cell “D27”: =D28+B27*B11/B14

Cell “E27”: =E28+C27*B11/B14

Cell “F27”: =F28-B27*B11/B15

Cell “G27”: =G28-C27*B11/B15

#### The calculation area – Active Formulas (continuation):

Cell “H27”: =H28+D27*B11

Cell “I27”: =I28+E27*B11

Cell “J27”: =J28+F27*B11

Cell “K27”: =K28+G27*B11

Cell “L27”: =L28+B11

#### Create the chart:

Click on an empty cell -> Insert -> Chart -> XY Scatter -> Finish

##### Create the chart – continuation

Delete Legend and Title.

Initial conditions

45

40

35

30

25

20

Initial conditions

15

10

5

0

-5 0 2 4 6 8 10 12

-10

##### Create the chart – continuation

– Click on any grid lines and delete them then change the charting area color to your taste.

– Change the font on both axes to 1 (using Format Axis).

– Click the white area -> SourceData Series.

25

20

15

– Remove whatever series is there and Add the following four series (type them in) ->

-5 0 2 4 6 8 10 12

Generate chart data series

Generate chart data series

##### Create the chart – continuation

In order to have the Zoom macro work right you have to name the chart “Chart 1”

You do this by selecting the chart and running the following macro:

Sub rename()

With ActiveChart

.Parent.Name = “Chart 1”

End With

End Sub

### Finish:

Click the button “Reset” and then “Zoom”.

This will reset the calculations and resize the chart to the proper zoom level.

Click the button “Trail_Size” up and down once.

Then click the button “Start-Pause”.

Wait few seconds and click “Start-Pause” again.

Now go on the chart and adjust the colors and sizes of the bodies and their trails to your liking.

### Tips:

You can experiment with various combinations of masses and initial conditions.

If dt is too large de simulation won’t converge.

You can change the zoom, trail size and even k during the simulation.

Be careful with k, change it slowly during the run. If you loose the bodies Reset the system and start it with the new k.

If you get stuck ask for help. In this case or if you just want do avoid the pain of creating one I can send you a copy of the original file.

by George Lungu <www.excelunusual.com>

Evgeny, Thanks for reporting. You were right, in 2007 the chart was not updating. I fixed it by adding “DoEvents” statement in 2 more places. The uploaded file is usable now in 2007. In 2007 it’s going to be roughly 10 times slower that in 2003 though (welcome to MS innovation).

Also if someone doesn’t know how to turn the Macros on here it is: 1) Left click the round “Office Button”, 2) Left click “Options” (bottom right), 3) Left click “Trust center” on the left side of the pop up, 4) Left Click “Trust Center Settings”, 5) Left click “Macro setings” 6) Click “Enable all macros….” 7) Exit all menus by clicking OK, OK.

There is some troubles under 2007 excel. When i click once on start-pause button – there is no moving on chart. It have to be added (Application.ScreenUpdating = True) before the first (loop) in VBA code. In that way it looks fine.

Thanks a lot, very interesting website. I just was amazed with roller-coaster))