Krdytkyu

I am currently trying to get a two body problem to work, that I can then upgrade to more planets, but it is not working. It is outputting me impossible positions. Does anyone know what is causing that?

This is the code I use:

day = 60*60*24

# Constants

G = 6.67408e-11

dt = 0.1*day

au = 1.496e8

t = 0





class CelBody:



    def __init__(self, id, name, x0, y0, z0, vx0, vy0, vz0, mass, vector, ax0, ay0, az0, totalforcex, totalforcey, totalforcez):

        self.ax0 = ax0

        self.ay0 = ay0

        self.az0 = az0



        self.ax = self.ax0

        self.ay = self.ay0

        self.az = self.az0



        # Constants of nature

        # Universal constant of gravitation

        self.G = 6.67408e-11

        # Name of the body (string)

        self.id = id

        self.name = name

        # Initial position of the body (au)

        self.x0 = x0

        self.y0 = y0

        self.z0 = z0

        # Position (au). Set to initial value.

        self.x = self.x0

        self.y = self.y0

        self.z = self.z0

        # Initial velocity of the body (au/s)

        self.vx0 = vx0

        self.vy0 = vy0

        self.vz0 = vz0

        # Velocity (au/s). Set to initial value.

        self.vx = self.vx0

        self.vy = self.vy0

        self.vz = self.vz0

        # Mass of the body (kg)

        self.M = mass

        # Short name

        self.vector = vector



        self.totalforcex = totalforcex

        self.totalforcey = totalforcey

        self.totalforcez = totalforcez



# All Celestial Bodies



forcex = 0

forcey = 0

forcez = 0



Bodies = [

    CelBody(0, 'Sun', 1, 1, 1, 0, 0, 0, 1.989e30, 'sun', 0, 0, 0, 0, 0, 0),

    CelBody(1, 'Mercury', 1*au, 1, 1, 0, 0.25920, 0, 3.3e23, 'earth', 0, 0, 0, 0, 0, 0),

    ]



leftover_bin = 

templistx = 

templisty = 

templistz = 



for v in range(100):

    for n in range(len(Bodies)):

        #Need to initialize the bodies

        ax = Bodies[n].ax0

        ay = Bodies[n].ay0

        az = Bodies[n].az0



        vx = Bodies[n].vx0

        vy = Bodies[n].vy0

        vz = Bodies[n].vz0



        xx = Bodies[n].x0

        xy = Bodies[n].y0

        xz = Bodies[n].z0



        planetinit = Bodies[n]



        for x in range(len(Bodies)):

            # Temporary lists and initial conditions

            planet = Bodies[x]



            if (planet == planetinit):

                pass



            else:

                rx = Bodies[x].x - Bodies[n].x

                ry = Bodies[x].y - Bodies[n].y

                rz = Bodies[x].z - Bodies[n].z



                r3 = (rx**2+ry**2+rz**2)**1.5

                gravconst = G*Bodies[n].M*Bodies[x].M

                fx = -gravconst*rx/r3

                fy = -gravconst*ry/r3

                fz = -gravconst*rz/r3





                # Make a temporary list of the total forces and then add them to get the resulting force

                templistx.append(fx)

                templisty.append(fy)

                templistz.append(fz)



        forcex = sum(templistx)

        forcey = sum(templisty)

        forcez = sum(templistz)

        templistx.clear()

        templisty.clear()

        templistz.clear()



        x = int(Bodies[n].x) + int(Bodies[n].vx) * dt

        y = int(Bodies[n].y) + int(Bodies[n].vx) * dt

        z = int(Bodies[n].z) + int(Bodies[n].vz) * dt



        Bodies[n].x = x

        Bodies[n].y = y

        Bodies[n].z = z



        vx = int(Bodies[n].vx) + forcex/int(Bodies[n].M)*dt

        vy = int(Bodies[n].vy) + forcey/int(Bodies[n].M)*dt

        vz = int(Bodies[n].vz) + forcez/int(Bodies[n].M)*dt



        Bodies[n].vx = vx

        Bodies[n].vy = vy

        Bodies[n].vz = vz



        t += dt









print(Bodies[0].name)

print(Bodies[0].x)

print(Bodies[0].y)

print(Bodies[0].z)





print(Bodies[1].name)

print(Bodies[1].x)

print(Bodies[1].y)

print(Bodies[1].z)

Note: The problem is probably in the for loops or in the rounding of the sum of the arrays but I am not sure