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()

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: