ceres-sees-all/Uni-Computing
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Update program and add gitignore

Browse source
This commit is contained in:
Ceres 2025-10-12 15:49:48 +01:00
parent d6aa03309c
commit 00db3aea6f
2 changed files with 60 additions and 51 deletions
Download patch file Download diff file Expand all files Collapse all files

1
.gitignore vendored Normal file
View file

@ -0,0 +1 @@
.vscode/settings.json

110
Exercise 1/exercise1.py
View file

@ -1,68 +1,76 @@
from math import *
import numpy as np import numpy as np
from scipy.integrate import solve_ivp from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
def part1(radFactor, orbits):
#Part one Diffeq
def f_part1(t, state, Me, Mm, G):
xm, ym, vx, vy = state
#Part one Diffeq dvxdt = -(Me*G*xm)/((xm**2+ym**2)**(3/2))
def f_part1(t, state, Me, Mm, G): dxmdt = vx
xm, ym, vx, vy = state dvydt = -(Me*G*ym)/((xm**2+ym**2)**(3/2))
dymdt = vy
dvxdt = -(Me*G*xm)/((xm**2+ym**2)**(3/2)) return (dxmdt, dymdt, dvxdt, dvydt)
dxmdt = vx
dvydt = -(Me*G*ym)/((xm**2+ym**2)**(3/2))
dymdt = vy
return (dxmdt, dymdt, dvxdt, dvydt)
# Initial Conditions # Initial Conditions
t_min = 0 Me = 5.97*(10**24)
t_max = 86400*20 Mm = 7.35*(10**22)
numpoints = 200 G = 6.67*(10**-11)
t = np.linspace(t_min, t_max, numpoints)
xm0 = 384400000 t_min = 0
ym0 = 0 t_max = 2360620*int(orbits)
vx0 = 0 numpoints = 2000*int(orbits)
vy0 = 1017.79 t = np.linspace(t_min, t_max, numpoints)
Me = 5.97*(10**24) r = 384400000
Mm = 7.35*(10**22) v = sqrt((G*Me)/r)*float(radFactor)
G = 6.67*(10**-11)
xm0 = r
ym0 = 0
vx0 = 0
vy0 = v
rtol = 1e-5
atol = 1e-8
#Solver #Solver
results = solve_ivp(f_part1, (t_min,t_max), (xm0, ym0, vx0, vy0), args=(Me, Mm, G), t_eval=t) results = solve_ivp(f_part1, (t_min,t_max), (xm0, ym0, vx0, vy0), args=(Me, Mm, G), t_eval=t, atol=atol, rtol=rtol)
#Graph plotting #Graph plotting
ax=plt.axes() # This creates some axes, so that we ax=plt.axes() # This creates some axes, so that we
ax.set_aspect(1) # can set the aspect ratio to 1 i.e. ax.set_aspect(1) # can set the aspect ratio to 1 i.e.
# x and y axes are scaled equally. # x and y axes are scaled equally.
ax.set_xlabel("x coordinate (m)") # Must label axes (with ax.set_xlabel("x coordinate (m)") # Must label axes (with
ax.set_ylabel("y coordinate (m)") # units) and give ax.set_ylabel("y coordinate (m)") # units) and give
ax.set_title("Orbit of moon around earth") # plot title. ax.set_title("Orbit of moon around earth") # plot title.
# #
ax.plot(results.y[0],results.y[1]) # Make the plot ax.plot(results.y[0],results.y[1]) # Make the plot
ax.legend(['Moon']) # and add a key. ax.plot(0,0)
plt.show() ax.legend(['Moon']) # and add a key.
plt.show()
# MyInput = '0' MyInput = '0'
# while MyInput != 'q': while MyInput != 'q':
# MyInput = input('Enter a choice, "1", "2" or "q" to quit: ') MyInput = input('Enter a choice, "1", "2" or "q" to quit: ')
# print('You entered the choice: ',MyInput) print('You entered the choice: ',MyInput)
# if MyInput == '1': if MyInput == '1':
# print('You have chosen part (1): simulation of a lunar orbit') print('You have chosen part (1): simulation of a lunar orbit')
# # radFactor = input(f'Enter the velocity scalar. A value of 1 will result in a circular orbit, and anything else will give an elliptical orbit: ')
# # put your code for part (1) here orbits = input(f'Enter the approximate number of orbits desired: ')
# # part1(radFactor, orbits)
# elif MyInput == '2': elif MyInput == '2':
# print('You have chosen part (2): earth-moon-probe system') print('You have chosen part (2): earth-moon-probe system')
# # #
# # put your code for part (2) here # put your code for part (2) here
# # #
# elif MyInput != 'q': elif MyInput != 'q':
# print('This is not a valid choice') print('This is not a valid choice')
# print('You have chosen to finish - goodbye.') print('You have chosen to finish - goodbye.')