Issue
Intro to the task
I am working on a simulating the orbit of moon Phobos around the planet Mars. I have completed the task of using numerical integration to update the velocity and position of both bodies. My final task is to produce an animated plot of the orbit of Phobos, with Mars centred at its initial position.
Code snippet
# plot the animated orbit
def plot_orbit(self):
# load all the data from the json file
all_data = self.load_data_for_all_bodies()
# get the positions and velocities as 2d vectors
positions, velocities = self.simulate()
# Create a figure and axis
fig = plt.figure()
ax = plt.axes()
mars_initial_pos_x = all_data["Mars"]["initial_position"][0]
mars_initial_pos_y = all_data["Mars"]["initial_position"][1]
print("mars: ", mars_initial_pos_x, mars_initial_pos_y)
phobos_initial_pos_x = all_data["Phobos"]["initial_position"][0]
phobos_initial_pos_y = all_data["Phobos"]["initial_position"][1]
print("phobos: ", phobos_initial_pos_x, phobos_initial_pos_y)
# Set the limits of the axes
ax.set_xlim(-9377300, 9377300)
# mars as a red circle
mars_circle = plt.Circle(
(mars_initial_pos_x, mars_initial_pos_y),
4,
color="red",
label="Mars",
animated=True,
)
ax.add_patch(mars_circle)
# Phobos as a blue circle
phobos_circle = plt.Circle(
(phobos_initial_pos_x, phobos_initial_pos_y),
80,
color="blue",
label="Phobos",
animated=True,
)
ax.add_patch(phobos_circle)
# Function to update the position of Phobos
def animate(frame):
phobos_circle.center = (positions[frame][0], positions[frame][1])
return (phobos_circle,)
# Create the animation
ani = FuncAnimation(
fig, animate, frames=self.iteration_counts, interval=20, blit=True
)
plt.title(f"Orbit of {self.name}")
plt.xlabel("x (m)")
plt.ylabel("y (m)")
plt.legend()
plt.show()
Data file(JSON)
{
"Mars": {
"mass": 6.4185e+23,
"initial_position": [
0,
0
],
"initial_velocity": [
0,
0
],
"orbital_radius": null
},
"Phobos": {
"mass": 1.06e+16,
"orbital_radius": 9377300.0,
"initial_position": [
9377300.0,
0
],
"initial_velocity": [
0,
2137.326983980658
]
}
}
Output
In the above output I can see that the positions and velocities array has been correctly filled because of my previous methods. However, the graph that is produced is a complete disaster. I ran the code in Jupyter notebook and VS code, Jupyter notebook does this and VS code prints nothing out on the screen. My first step is to get circles on the graph in the correct positions initially and then think about the animate function. Sadly, I can't see any circles in the output. I used very similar code that was given in the answer to a question about FuncAnimation on stack.
Any help would be appreciated!
Edit
revised code(executable)
# include neccesary libraries
import numpy as np
import matplotlib.pyplot as plt
import json
from matplotlib.animation import FuncAnimation
# define constants
timestep = 0.001
iterations = 1000
G = 6.674 * 10**-11
def main():
# make the modification
modify_initial_velocity_of_phobos()
# Example usage
mars = Planet("Mars")
phobos = Planet("Phobos")
# Display information
mars.display_info()
phobos.display_info()
# mars.acceleration()
# phobos.acceleration()
phobos_orbit = Simulate("Phobos", timestep, iterations)
phobos_orbit.plot_orbit()
# Function to calculate the initial velocity of a moon orbiting a planet
def calculate_orbital_velocity(mass, radius):
# use the formula given to calculate
return (G * mass / radius) ** 0.5
# self explanatory function
def modify_initial_velocity_of_phobos():
# Read data from the JSON file
with open("planets.json", "r+") as file:
celestial_data = json.load(file)
# Extract necessary data for calculation
mass_of_mars = celestial_data["Mars"]["mass"]
orbital_radius_of_phobos = celestial_data["Phobos"]["orbital_radius"]
# Calculate orbital velocity of Phobos
velocity_of_phobos = calculate_orbital_velocity(
mass_of_mars, orbital_radius_of_phobos
)
# Update the initial_velocity for Phobos in the data
celestial_data["Phobos"]["initial_velocity"] = [0, velocity_of_phobos]
# Move the file pointer is back to the start of the file
file.seek(0)
# Write the updated data back to the JSON file
json.dump(celestial_data, file, indent=4)
# Truncate the file to remove any leftover data
file.truncate()
# create a class for the planet which calculates the net gravitational force and acceleration due to it acting on it
class Planet():
def __init__(self, name):
self.name = name
body_data = self.load_data_for_main_body()
# Initialize attributes with data from JSON or default values
self.mass = body_data.get("mass")
self.position = np.array(body_data.get("initial_position"))
self.velocity = np.array(body_data.get("initial_velocity"))
self.orbital_radius = body_data.get("orbital_radius", 0)
# load main planet(body) data from the json file
def load_data_for_main_body(self):
with open("planets.json", "r") as file:
all_data = json.load(file)
body_data = all_data.get(self.name, {})
return body_data
# load all data from the json file
def load_data_for_all_bodies(self):
with open("planets.json", "r") as file:
all_data = json.load(file)
return all_data
# calculate the gravitational force between two bodies
def force(self):
all_bodies = self.load_data_for_all_bodies()
# initialize the total force vector
total_force = np.array([0.0, 0.0])
# iterate over all the bodies
for body_name, body_data in all_bodies.items():
if body_name == self.name:
continue
# get the position of each body
other_body_position = np.array(body_data["initial_position"])
# get the mass of each body
other_body_mass = body_data["mass"]
# calculate distance vector between the two bodies
r = other_body_position - self.position
# Calculate the distance between the two bodies
mag_r = np.linalg.norm(r)
# Normalize the distance vector
r_hat = r / mag_r
# Calculate the force vector between the two bodies
force = (G * self.mass * other_body_mass) / (mag_r**2) * r_hat
# Add the force vector to the total force vector
total_force += force
return total_force
# calculate the acceleration due to the force
def acceleration(self):
# Calculate the force vector
force = self.force()
# Calculate the acceleration vector
acceleration = force / self.mass
return acceleration
# update the position of the body using the velocity and time step
def update_position(self):
self.position = self.position + self.velocity * timestep
return self.position
# update the velocity of the body using the acceleration and time step
def update_velocity(self):
self.velocity = self.velocity + self.acceleration() * timestep
return self.velocity
def display_info(self):
print(f"Name: {self.name}")
print(f"Mass: {self.mass} kg")
print(f"Position: {self.position} m")
print(f"Velocity: {self.velocity} m/s")
if self.orbital_radius is not None:
print(f"Orbital Radius: {self.orbital_radius} m")
else:
print("Orbital Radius: Not applicable")
# define class to simulate the orbit of the moon around the planet
class Simulate(Planet):
def __init__(self, name, delta_t, iteration_counts):
super().__init__(name)
self.delta_t = delta_t
self.iteration_counts = iteration_counts
def simulate(self):
global timestep
# initialize the arrays to store the positions and velocities as vectors
positions = np.zeros((self.iteration_counts, 2))
velocities = np.zeros((self.iteration_counts, 2))
# iterate over the number of iterations
for i in range(self.iteration_counts):
# update the position
positions[i] = self.update_position()
# update the velocity
velocities[i] = self.update_velocity()
# update time
timestep += self.delta_t
print("pos: ", positions)
print("vel: ", velocities)
return positions, velocities
# plot the animated orbit
def plot_orbit(self):
# load all the data from the json file
all_data = self.load_data_for_all_bodies()
# get the positions and velocities as vectors
positions, velocities = self.simulate()
# debug statements
print("pos: ", positions)
print("vel: ", velocities)
# Create a figure and axis
fig, ax = plt.subplots()
mars_initial_pos_x = all_data["Mars"]["initial_position"][0]
mars_initial_pos_y = all_data["Mars"]["initial_position"][1]
print("mars: ", mars_initial_pos_x, mars_initial_pos_y)
phobos_initial_pos_x = all_data["Phobos"]["initial_position"][0]
phobos_initial_pos_y = all_data["Phobos"]["initial_position"][1]
print("phobos: ", phobos_initial_pos_x, phobos_initial_pos_y)
# Set the limits of the axes
ax.set_xlim(-9377300, 9377300)
# mars as a red circle
mars_circle = plt.Circle(
(mars_initial_pos_x, mars_initial_pos_y),
0.1,
color="red",
label="Mars",
ec = "None"
)
ax.add_patch(mars_circle)
# Phobos as a blue circle
phobos_circle = plt.Circle(
(phobos_initial_pos_x, phobos_initial_pos_y),
0.1,
color="blue",
label="Phobos",
ec = "None"
)
ax.add_patch(phobos_circle)
# Function to update the position of Phobos
def animate(frame):
phobos_circle.center = (positions[frame][0], positions[frame][1])
ax.draw_artist(phobos_circle)
return phobos_circle,
# Create the animation
ani = FuncAnimation(
fig = fig, func = animate, frames=self.iteration_counts, interval=20, blit=True
)
plt.title(f"Orbit of {self.name}")
plt.xlabel("x (m)")
plt.ylabel("y (m)")
plt.legend()
plt.show()
fig.savefig('MyFigure.png', dpi=200)
# plot the orbit
# plt.plot(positions[:, 0], positions[:, 1])
# plt.title(f"Orbit of {self.name}")
# plt.xlabel("x(m)")
# plt.ylabel("y(m)")
# plt.show()
if __name__ == "__main__":
main()
Solution
ok... with the full code it was a lot easier to see where things go wrong :-)
The main problems are:
- your circles are extremely small!
- you set an enormous x-extent but keep the y-extent to 0-1
(since your axis does not preserve the aspect-ratio (e.g.
ax.set_aspect("equal")) this means your circles look like lines... - your motion is extremely slow
- your object gets garbage-collected as soon as your
main()function has resolved so the animation cannot run
To fix your code, this is what you have to do:
Give Mars and Phobos a proper radius
mars_circle = plt.Circle(
(mars_initial_pos_x, mars_initial_pos_y),
1e5,
color="red",
label="Mars",
ec = "None"
)
# Phobos as a blue circle
phobos_circle = plt.Circle(
(phobos_initial_pos_x, phobos_initial_pos_y),
1e5,
color="blue",
label="Phobos",
ec = "None",
)
make sure your axis-limits are properly set
either
ax.set_xlim(-9377300, 9377300)
ax.set_ylim(-9377300, 9377300)
or
ax.set_xlim(-9377300, 9377300)
ax.set_aspect("equal")
reduce the interval to something you can actually see (otherwise your movement will be extremely slow)
self.ani = FuncAnimation(
fig = fig, func = animate, frames=self.iteration_counts, interval=.5, blit=True
)
make sure the animation object is not garbage-collected
# Create the animation
self.ani = FuncAnimation(
fig = fig, func = animate, frames=self.iteration_counts, interval=.5, blit=True
)
and
def main():
...
...
phobos_orbit = Simulate("Phobos", timestep, iterations)
phobos_orbit.plot_orbit()
return phobos_orbit
and
if __name__ == "__main__":
phobos_orbit = main()
... and you apparently don't need an explicit call to ax.draw_artist()...
If you change all that, it seems to do the job nicely:
Answered By - raphael
0 comments:
Post a Comment
Note: Only a member of this blog may post a comment.