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

Improve Comments

Browse source
This commit is contained in:
Ceres 2025-10-20 15:55:14 +01:00
parent ea1b5a04f7
commit bd42f0f7e6
Signed by: ceres-sees-all
GPG key ID: 9814758436430045
1 changed files with 22 additions and 21 deletions
Download patch file Download diff file Expand all files Collapse all files

43
Exercise 1/exercise1.py
View file

@ -4,38 +4,40 @@ from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
def part1(radFactor, orbits): def part1(velocityFactor, orbits):
#Part one Diffeq #Part one Diffeq
def f_part1(t, state, Me, Mm, G): def f_part1(t, state, Me, Mm, G):
xm, ym, vx, vy = state xm, ym, vx, vy = state #all input values
dvxdt = -(Me*G*xm)/((xm**2+ym**2)**(3/2)) dxmdt = vx #functions for each diffeq
dxmdt = vx
dvydt = -(Me*G*ym)/((xm**2+ym**2)**(3/2))
dymdt = vy dymdt = vy
dvxdt = -(Me*G*xm)/((xm**2+ym**2)**(3/2))
dvydt = -(Me*G*ym)/((xm**2+ym**2)**(3/2))
return (dxmdt, dymdt, dvxdt, dvydt)
return (dxmdt, dymdt, dvxdt, dvydt) #return differentiated values
# Initial Conditions # Initial Conditions
Me = 5.97*(10**24) Me = 5.97*(10**24) #Constants
Mm = 7.35*(10**22) Mm = 7.35*(10**22)
G = 6.67*(10**-11) G = 6.67*(10**-11)
t_min = 0 t_min = 0
t_max = 2360620*int(orbits) t_max = 2360620*int(orbits) #Time for one orbit * number of orbits
numpoints = 2000*int(orbits) numpoints = 2000*int(orbits) #To keep good clarity when plotting
t = np.linspace(t_min, t_max, numpoints) t = np.linspace(t_min, t_max, numpoints)
r = 384400000 r = 384400000 #Distance or orbit
v = sqrt((G*Me)/r)*float(radFactor) v = sqrt((G*Me)/r)*float(velocityFactor)#Balancing centrepetal force and gravitational force,
#then multiplying by a scalar to obtain an elliptical orbit
xm0 = r xm0 = r #Initial position and velocity of moon
ym0 = 0 ym0 = 0 #Set so moon is travelling only in the y direction initially
vx0 = 0 vx0 = 0
vy0 = v vy0 = v
rtol = 1e-5 rtol = 1e-5 #Tolerences set to balance computation speed with accuracy
atol = 1e-8 atol = 1e-8
@ -58,7 +60,7 @@ def part1(radFactor, orbits):
ax.axvline(x=0, color='k', linewidth=0.5) ax.axvline(x=0, color='k', linewidth=0.5)
plt.show() plt.show()
def part2(radFactor, orbits): def part2(velocityFactor, orbits):
#Part one Diffeq #Part one Diffeq
def f_part1(t, state, Me, Mm, G): def f_part1(t, state, Me, Mm, G):
xm, ym, vx, vy, xp, yp, vpx, vpy = state xm, ym, vx, vy, xp, yp, vpx, vpy = state
@ -90,7 +92,7 @@ def part2(radFactor, orbits):
t = np.linspace(t_min, t_max, numpoints) t = np.linspace(t_min, t_max, numpoints)
rm = 384400000 rm = 384400000
vm = sqrt((G*Me)/rm)*float(radFactor) vm = sqrt((G*Me)/rm)*float(velocityFactor)
rpm = 10000000 rpm = 10000000
vpm = sqrt((G*Mm)/rpm) vpm = sqrt((G*Mm)/rpm)
@ -108,7 +110,6 @@ def part2(radFactor, orbits):
rtol = 1e-6 rtol = 1e-6
atol = 1e-9 atol = 1e-9
#Solver #Solver
results = solve_ivp(f_part1, (t_min,t_max), (xm0, ym0, vx0, vy0, xp0, yp0, vpx0, vpy0), args=(Me, Mm, G), t_eval=t, atol=atol, rtol=rtol) results = solve_ivp(f_part1, (t_min,t_max), (xm0, ym0, vx0, vy0, xp0, yp0, vpx0, vpy0), args=(Me, Mm, G), t_eval=t, atol=atol, rtol=rtol)
@ -134,14 +135,14 @@ while MyInput != 'q':
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: ') velocityFactor = input(f'Enter the velocity scalar. A value of 1 will result in a circular orbit, and anything else will give an elliptical orbit: ')
orbits = input(f'Enter the approximate number of orbits desired: ') orbits = input(f'Enter the approximate number of orbits desired: ')
part1(radFactor, orbits) part1(velocityFactor, 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')
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: ') velocityFactor = input(f'Enter the velocity scalar. A value of 1 will result in a circular orbit, and anything else will give an elliptical orbit: ')
orbits = input(f'Enter the approximate number of orbits desired: ') orbits = input(f'Enter the approximate number of orbits desired: ')
part2(radFactor, orbits) part2(velocityFactor, orbits)
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.')