Intro to the task
我正在模拟火卫一围绕火星的轨道.我已经完成了使用数值积分来更新两个物体的速度和位置的任务.我的最后一个任务是制作火卫一轨道的动画图，火星以其初始位置为中心.
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
在上面的输出中，我可以看到由于我前面的方法，位置和速度数组已经正确填充.然而，生成的图表是一场彻底的灾难.我在Jupyter笔记本和VS代码中运行了代码，Jupyter笔记本做到了这一点，VS代码没有在屏幕上打印出来.我的第一步是首先在图形上将圆放在正确的位置，然后考虑动画功能.遗憾的是，我在输出中看不到任何圆圈.我使用的代码与在回答有关堆栈上的FuncAnimation的问题时给出的代码非常相似.
如有任何帮助，我们将不胜感激！
Edit
修改后的代码(可执行文件)
# 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()
