0

I'm trying to make a basic gravity sim using py but for some reason it plots all the lines as straight, I looked at a couple of examples when I got stuck but they all use similar equations / ways of plotting the data so not sure where I went wrong with it

class Body:

    # A class to initialise a body of mass to plot

    def __init__(self, id, mass, coordinates, velocities):
        self.id = id
        self.coordinates = np.array(coordinates, dtype=float)
        self.v = np.array(velocities, dtype=float)
        self.mass = mass
        self.a = np.zeros(3, dtype=float)
        MOTION_LOG.append({"name": self.id, "x":[coordinates[0]], "y": [coordinates[1]], "z": [coordinates[2]]})


    # Procedure for checking gravity effects on body
    def gravity(self):
        self.a = 0
        for body in bodies:
            if body != self:
                dist = body.coordinates - self.coordinates
                r = np.sqrt(np.sum(dist**2))
                self.a += (SETTINGS['G'] * body.mass * dist / r**3) ** 2


# Procedure to plot the new coordinates of the body         
    def move(self):
        self.v += self.a * SETTINGS["deltaT"]
        self.coordinates += self.v * SETTINGS['deltaT']

Then to actually simulate it I did

# For loop to run a simulation for a specific time set
for step in range(int(SETTINGS["tLimit"] / SETTINGS["deltaT"])):
    SETTINGS['elapsedT'] += SETTINGS['deltaT']
    if SETTINGS['elapsedT'] % SETTINGS["frequency"] == 0:
        prog = ((SETTINGS['elapsedT'] / SETTINGS['tLimit'])*100)//1
        for index, location in enumerate(MOTION_LOG):
            location["x"].append(bodies[index].coordinates[0])
            location["y"].append(bodies[index].coordinates[1])
            location["z"].append(bodies[index].coordinates[2])
        print(f"{prog}%")

    for body in bodies:
        body.gravity()
    for body in bodies:
        body.move()

Then finally for plotting it on the graph I did

fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection='3d')

for body in MOTION_LOG:
    ax.plot(body["x"], body["y"], body["z"])

plt.show()

Not sure if I've just done the acceleration wrong or I'm simply plotting the points wrong but other examples I've seen don't do things much differently

Example data

SETTINGS = {
    'G' : 6.67e-11,
    'deltaT' : 172800,
    'elapsedT' : 0,
    'tLimit' : 315360000,
    "frequency": 1,
}

MOTION_LOG = []
bodies = []
set_bodies = {
    "Sun": {
        "id": "Sun",
        "mass": 1e20,
        "coordinates": [0, 0, 0],
        "velocities": [0, 0, 0]
    },
    "Earth": {
        "id": "Earth",
        "mass": 1e3,
        "coordinates": [1e2, -1e2, 0],
        "velocities": [0, 0, 0]
    }
}
for body, body_data in set_bodies.items():
    bodies.append(Body(**body_data))
AJ Biffl
  • 493
  • 3
  • 10

1 Answers1

1

Firstly, your acceleration is squared. The line you have is:

self.a += (SETTINGS['G'] * body.mass * dist / r**3) ** 2

What you want is:

self.a += (SETTINGS['G'] * body.mass * dist / r**3) 

Second, your numbers are all out of whack. Never, never just guess at nice combinations of numbers to use for simulations like this. Either use completely natural units (where the grav constant, solar mass, and AU are all equal to one), or use all the actual numbers. When you don't do that, you're just guessing at what timesteps are appropriate, and it's really difficult to get that right. For instance, with the numbers you're using for the Solar mass, grav constant, and Earth orbital radius, one unit of time is about 7.7 years. To get something that looks decent, then you'd want a timestep of about 7e-4 and run until a final time of about 1.3, and use an initial velocity of [5000, 5000, 0] for Earth (you'll also need to change how you add to MOTION_LOG since small floats don't play well with mod (%)). But don't do that - use sensible numbers. You'll be happier in the end when you don't have to spend hours upon hours converting back and forth into units you actually want.

Third, your "Earth" object doesn't have an initial velocity. It's falling into the Sun, no matter how you fix your equations or your units.

Here are some about-right "real" numbers using SI units that will work:

set_bodies = {
    "Sun": {
        "id": "Sun",
        "mass": 2e30,
        "coordinates": [0, 0, 0],
        "velocities": [0, 0, 0]
    },
    "Earth": {
        "id": "Earth",
        "mass": 6e24,
        "coordinates": [1.5e11, 0, 0],
        "velocities": [0, 3e4, 0]
    }
}
AJ Biffl
  • 493
  • 3
  • 10
  • And a quick note - I generally recommend using natural units in almost any simulation, especially orbital systems, because you run much less risk of hitting machine precision or overflow/underflow issues with your numbers since they're all close to one (instead of being, well, astronomical) – AJ Biffl Dec 19 '21 at 23:27
  • I didn't tend to use the test data I provided anyways, I used data from [here](https://ssd.jpl.nasa.gov/api/horizons.api) which I scraped, but even if I change the values to more natural ones I still get straight lines on the graph – Epicyon Haydeni Dec 19 '21 at 23:33
  • Did you try the numbers I gave in my answer? – AJ Biffl Dec 19 '21 at 23:35
  • Oh, I forgot to mention another error I found. See edits in answer. – AJ Biffl Dec 19 '21 at 23:39
  • oh yeah forgot to remove the square when I pasted it in, I changed around my equation a bit and remember someone mentioning to square it I must of forgot to remove that before adding it into the embed. My actual equation is `self.a += -(SETTINGS["G"] * body.mass) / r * dist` I just can't edit it. even with those numbers I get two straight lines that repel each other – Epicyon Haydeni Dec 19 '21 at 23:47
  • You shouldn't have the negative sign, and the `r` needs to be cubed still: `self.a += (SETTINGS["G"] * body.mass) / (r**3) * dist` – AJ Biffl Dec 19 '21 at 23:48
  • fair enough I was just trying different equations to try to get it to work, similar result though, I did used to use that equation that you posted just didn't seem to work – Epicyon Haydeni Dec 19 '21 at 23:49