This repository was archived by the owner on Dec 21, 2021. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathSolarSystemProject.py
More file actions
280 lines (240 loc) · 12.9 KB
/
Copy pathSolarSystemProject.py
File metadata and controls
280 lines (240 loc) · 12.9 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
"""
Combined.py
Date: 23/1/20
Authors:
Seppe Lampe
Petio Todorov
This is a simplified model of the Solar System created
using the Vpython library. The program calculates
the forces between the celestial bodies in the system
in order calculate the position of each body
as a function of time.
Use the Main.py program to run the simulation.
"""
from math import *
from vpython import *
G = 6.67 * 10 ** -11
"""
General Notes:
All mass has units of kg
All radii have units of meters
orbitalPeriod is the orbital period, in Earth days (86400 seconds)
"""
class CelestialBody:
"""CelestialBody is the superclass of all types of celestial bodies"""
def __init__(self, mass, radius, color, name):
self.mass = mass
self.radius = radius
self.color = color
self.name = name
class Star(CelestialBody):
"""
A star is a celestial body without a parent. The properties of the Sun of our Solar System
are set as the default arguments.
"""
def __init__(self, mass=1.989 * 10**30, radius=6.95 * 10**8, color=color.yellow, name = "Sun"):
super(Star, self).__init__(mass, radius, color, name)
class Satellite(CelestialBody):
"""
This is a class that represents a body which revolves around another body.
Parameters:
parent (CelestialBody): is the object that the body revolves around
distanceToParent (double): distance (in meters) between the body and the parent
orbital period (float): is the orbital period of the body around the parent (in Earth days (86400 seconds))
inclination (float): this is the angle (in degrees) between the plane in which
the body revolves and the plane going through the central axes (xy) of the parent
"""
def __init__(self, mass, radius, color, name, parent, distanceToParent, orbitalPeriod,
inclination=0.0): # orbitalPeriod is the orbital period, in Earth days (86400 seconds)
super(Satellite, self).__init__(mass, radius, color, name)
self.parent = parent
self.distanceToParent = distanceToParent
self.orbitalPeriod = orbitalPeriod
self.inclination = radians(inclination)
def gforce(body1, body2):
"""
NOTE:
This function is adapted from the following source:
https://www.glowscript.org/#/user/wlane/folder/Let'sCodePhysics/program/Solar-System-1/edit
Glowscript is a website that has many relevant examples for Vpython.
This is a function that calculates the gravitational force exerted between two objects.
The function computes the force exerted ON body1 BY body2.
General formula for the force of gravity: Fgravity = G*M1*M2/d^2
General formula for the force vector: Fvector = F*r^
Parameters:
body1 (CelestialBody): is a celestial body
body2 (CelestialBody): is a celestial body which is not the same body as body1
"""
dvector = body1.pos - body2.pos # Calculate distance vector between p1 and p2.
d = mag(dvector) # Calculate magnitude of distance vector.
r_hat = dvector / d # Calculate unit vector of distance vector.
Fgravity = G * body1.mass * body2.mass / d ** 2 # Calculate force magnitude.
return Fgravity * r_hat # Calculate force vector.
class RevolvingSystem:
"""This class represents a central body with satellites orbiting around it"""
def __init__(self, centralBody):
"""
The constructor the Revolving System class.
It keeps track of the objects in the system using a list, "self.objects".
The system also has a central body around which all other bodies revolve.
It also creates a new canvas in which objects can be displayed.
(Note no objects are added TO THE CANVAS in the constructor.)
"""
self.canvas = canvas(height=450, width=1000)
self.objects = []
self.centralBody = centralBody
centralBody.momentum = vector(0,0,0)
centralBody.pos = vector(0,0,0)
def addSatellite(self, body):
"""
This is a method that adds a satellite to the RevolvingSystem.
Because every satellite has a parent (central body), the checkParent
function is used to check if the parent (central body) of the satellite
is the same as the central body of the revolving system.
Otherwise, if the parent of the satellite is the same as the central body of the
revolving system then we update the bodies position and momentum and also add
it to the list of objects in the system.
Parameters:
body (CelestialBody): is a celestial body to be added to the system
"""
if not self.checkParent(body):
print("body has the wrong parent")
else:
self.generatePos(body)
self.generateMomentum(body)
self.objects.append(body)
def checkParent(self, body):
"""
This method checks if this particular body can be added to this particular
revolving system. It does so by recursively checking if the parent of the body
is the same as the central of the revolving system. (Imagine multiple layers of
nesting.)
Base Case:
If the parent of the body is the same as the central body of the system,
then the function returns True.
If the parent of the body is a star and it is not the same as the central body of this
revolving system, then the function returns False. For example, because Pluto is a
planet which revolves around the sun (parent) we should not be able to add it to a
revolving system where Jupiter is the central body.
Recursive Call:
In the case where the parent is not the same as the central body of the system
and where the parent is also not a Star, then we run checkParent on the parent.
For example, what happens if we want to add the moon to the Solar System? The parent
of the moon is the Earth and the central body of the system is the Sun. The moon's parent
is not the sun, so the first if condition is False. The moon's parent is also not a Star,
it's a planet, so the second if condition is False. Finally we run checkParent on the moon's
parent (i.e. the Earth). There the first if condition is True because the parent of the earth
is the sun so the function returns True meaning that it's acceptable for this body
to be added to this revolving system.
Parameters:
body (CelestialBody): is a celestial body to be added to the system
"""
if body.parent == self.centralBody:
return True
if isinstance(body.parent, Star):
return False
return self.checkParent(body.parent)
def generatePos(self, body):
"""
This method generates the initial position of each body in the Revolving System.
The position is based on the distance to its parent body and its inclination relative
to the parent body.
Parameters:
body (CelestialBody): is a celestial body
"""
if self.centralBody == body.parent:
body.pos = vector(body.distanceToParent * cos(body.inclination), 0,
body.distanceToParent * sin(body.inclination))
else:
if body.parent.pos == vector(0,0,0):
self.generatePos(body.parent)
body.pos = body.parent.pos + vector(body.distanceToParent * cos(body.inclination), 0,
body.distanceToParent * sin(body.inclination))
def generateMomentum(self, body):
"""
This method generates the initial momentum of each body in the Revolving System.
The momentum is based on the distance to its parent body, the current body's mass
and the current body's orbital period around its parent.
Parameters:
body (CelestialBody): is a celestial body
"""
if self.centralBody == body.parent:
body.momentum = vector(0, body.mass * ((body.distanceToParent * 2 * pi) / (body.orbitalPeriod * 86400)), 0)
else:
if body.parent.pos == vector(0, 0, 0):
self.generateMomentum(body.parent)
body.momentum = body.parent.momentum / body.parent.mass * body.mass + vector(0, body.mass * (
(body.distanceToParent * 2 * pi) / (body.orbitalPeriod * 86400)), 0)
def show(self, body):
"""
This method shows the body on the canvas. It does so by creating a Vpython sphere object.
Calling the sphere vpython function automatically creates the sphere physically on the canvas.
Parameters:
body (CelestialBody): is a celestial body
"""
body.visualisation = sphere(pos=body.pos, radius=body.radius, color=body.color,
mass=body.mass, momentum=body.momentum, force=vector(0,0,0), make_trail=True, name=body.name)
def simulate(self, dt):
"""
This is the fundamental method for showing and then updating the position of bodies in the revolving
system on the canvas.
First the method shows every body in the objects list of the revolving system on the canvas.
Then there is an infinite loop (which ends only when the user closes the simulation browser window).
Otherwise, the simulation continues to run indefinitely. In this loop, the forces, momentum, and then
position are updated for every object in the solar system except for the central body which remains
stationary.
An important element of the simulation method is the rate(100) function call. This is a function which
controls the animation speed of the canvas. In our case it will recompute the position only once every
0.01 seconds. The following link is the documentation of the Vpython rate function:
https://www.glowscript.org/docs/VPythonDocs/rate.html
Parameters:
dt (float): is the current time step of the simulation (units Earth days/second)
"""
# This converts the time step from Earth days to seconds.
# Normally this would mean multiplying it by 3600*24 but here we multiply it by 36*24.
# As the loop is later executed 100 times per second, defined by the rate(100).
dt *= 36 * 24
for x in range(len(self.objects)):
self.show(self.objects[x])
self.show(self.centralBody)
while True:
rate(100)
for x in range(len(self.canvas.objects)):
# This if statement is to prevent the sun from calculating a force which results eventually in a slight,
# continuous position shift for the sun
if self.canvas.objects[x].name != self.centralBody.name:
self.updateForce(self.canvas.objects[x])
self.updateMomentum(self.canvas.objects[x], dt)
self.updatePositions(dt)
def updateForce(self, body):
"""
This method computes the force on the input parameter body, by summing the forces on it from
all other bodies of the revolving system which are shown on the canvas, except for the body itself.
self.canvas.objects is a list of all of the objects shown on the Vpython canvas.
Parameters:
body (CelestialBody): is a celestial body
"""
body.force = vector(0, 0, 0)
for x in range(len(self.canvas.objects)):
if self.canvas.objects[x] != body:
body.force += gforce(self.canvas.objects[x], body)
def updateMomentum(self, body, dt):
"""
This method computes the change in momentum on the input parameter body at a particular time step, dt.
This is based on the formula F=ma where a = dv/dt so F = m*(dv/dt) so F*(dt) = m*(dv) = momentum.
Parameters:
body (CelestialBody): is a celestial body
dt (float): is the current time step of the simulation
"""
body.momentum += body.force * dt
def updatePositions(self, dt):
"""
This method computes the updated position of all objects on the canvas of this revolving system
for this current simulation time step, dt.
It does so by adding the change in position, momentum divided by mass multiplied by the time step, to
the current position of the object on the canvas.
Parameters:
dt (float): is the current time step of the simulation
"""
for x in range(len(self.canvas.objects)):
self.canvas.objects[x].pos += self.canvas.objects[x].momentum / self.canvas.objects[x].mass * dt