Python Forum

I'm trying to model the solar system, with the gravitational effects of each planet interacting with every other planet (including the Sun). So far I have all of the planets orbiting the Sun, but I can't figure out how to do the same for Mercury. My plan was to do exactly what I did for the sun, but I kind of cheated with the sun by excluding it from the loop and using a separate variable for the sun when getting my r value. I had to do this because I get a 'division by zero' error when I leave it in.

Is there a way to include the sun in 'i' for the loop, exclude it while calculating r for the suns gravity, then include it again when calculating Mercuries gravity?

from visual import *

myuniv = display(range=vector(7*10**12,7*10**12,7*10**12))

G = 6.674*10**-11  #Gravitational constant

#Masses of each solar body in order starting with the Sun
mass = [1.989*10**30,3.285*10**23, 4.867*10**24, 
           5.972*10**24, 6.39*10**23, 1.898*10**27,
           5.68*10**26, 8.68*10**25, 1.024*10**26]

#Positions of each solar body starting with the Sun           
posit = [vector(0,0,0),vector(5.791*10**10,0,0),vector(1.082*10**11,0,0),
          vector(1.496*10**11,0,0),vector(2.279*10**11,0,0),vector(7.785*10**11,0,0),
          vector(1.429*10**12,0,0),vector(2.871*10**12,0,0),vector(4.498*10**12,0,0)]

#Velocities of each solar body starting with the Sun
vel = [vector(0,0,0),vector(0,4.87*10**4,0),vector(0,3.502*10**4,0),
          vector(0,2.98*10**4,0),vector(0,2.401*10**4,0),vector(0,1.307*10**4,0),
          vector(0,9.69*10**3,0),vector(0,6.81*10**3,0),vector(0,5.43*10**3,0)]

#Acceleration of each solar body starting with the Sun
acc = [vector(0,0,0),vector(0,0,0),vector(0,0,0),
          vector(0,0,0),vector(0,0,0),vector(0,0,0),vector(0,0,0),
          vector(0,0,0),vector(0,0,0),vector(0,0,0),vector(0,0,0)]

Msun = 1.989*10**30
Ssun = vector(0,0,0)

t = 0.0
dt = 150000

sun = sphere(pos=posit[0],radius=9.9*10**10,color=color.yellow,make_trail=True)
mercury = sphere(pos=posit[1],radius=2.44*10**6,color=color.red,make_trail=True)
venus = sphere(pos=posit[2],radius=6.052*10**6,make_trail=True)
earth = sphere(pos=posit[3],radius=6.371*10**6,color=color.green,make_trail=True)
mars = sphere(pos=posit[4],radius=3.39*10**6,color=color.red,make_trail=True)
jupiter = sphere(pos=posit[5],radius=6.99*10**7,color=color.blue,make_trail=True)
saturn = sphere(pos=posit[6],radius=5.823*10**7,color=color.yellow,make_trail=True)
uranus = sphere(pos=posit[7],radius=2.536*10**7,color=color.cyan,make_trail=True)
neptune = sphere(pos=posit[8],radius=2.462*10**7,make_trail=True)

while True:
    for i in arange(1,9,1):
        rate(1000000)
        t=t+dt
        r = Ssun-posit[i]
        Fnet= ((G*Msun*mass[i])/(mag(r)**2))* norm(r) #Calculating the force of gravity from the sun
        
        acc[i]=Fnet/mass[i]
        vel[i]=vel[i] + acc[i]*dt
        posit[i]=posit[i] + vel[i]*dt
        
        sun.pos=posit[0]
        mercury.pos=posit[1]
        venus.pos=posit[2]
        earth.pos=posit[3]
        mars.pos=posit[4]
        jupiter.pos=posit[5]
        saturn.pos=posit[6]
        uranus.pos=posit[7]
        neptune.pos=posit[8]

Hopefully this made sense. For the most part I know what I'm doing, but I'm still very new to Python.

Thanks!

---

By the way, this is for Python v.2.7.9, and I don't have numpy. (I can't use it for this particular class)

the sun and mercury should be like any other body in the system. you know each is affected by all the others. just add up the vectors in each cycle and apply it (minus the vector of what it applies to) to each and repeat. another hint is to sort your list by body mass, so you are adding in the order small to large. that helps avoid some rounding errors.

Nice task you have there! I have played with orbital mechanics simulations before and I expect you are not going to like results after many iterations due to numerical errors.
If you are ready to delve into more advanced stuff, an improvement would be to use a higher order Runge-Kutta solver, for example. Even better improvement would be using a symplectic integrator, which preserves energy in the system over time.

Very interesting indeed.
I changed your script a bit to make it more easy to understand and follow

from visual import * # I HATE this
 
myuniv = display(range=vector(7*10**12,7*10**12,7*10**12))
 
G = 6.674*10**-11  #Gravitational constant

# define objects
sun = {'r':9.9*10**10, 'mass':1.989*10**30, 'pos':vector(0,0,0), 'vel':vector(0,0,0), 'acc':vector(0,0,0), 'color':color.yellow}
planets = (('mercury', {'r':2.44*10**6, 'mass':3.285*10**23,'pos':vector(5.791*10**10,0,0), 'vel': vector(0,4.87*10**4,0), 'acc':vector(0,0,0), 'color':color.red}),
           ('venus', {'r':6.052*10**6, 'mass':4.867*10**24, 'pos':vector(1.082*10**11,0,0), 'vel':vector(0,3.502*10**4,0), 'acc':vector(0,0,0), 'color':color.yellow}), # I added color here
           ('earth', {'r':6.371*10**6, 'mass':5.972*10**24, 'pos':vector(1.496*10**11,0,0), 'vel':vector(0,2.98*10**4,0), 'acc':vector(0,0,0), 'color':color.green}),
           ('mars', {'r':3.39*10**6, 'mass':6.39*10**23, 'pos':vector(2.279*10**11,0,0), 'vel':vector(0,2.401*10**4,0), 'acc':vector(0,0,0), 'color':color.red}),
           ('jupiter', {'r':6.99*10**7, 'mass':1.898*10**27, 'pos':vector(7.785*10**11,0,0), 'vel':vector(0,1.307*10**4,0), 'acc':vector(0,0,0), 'color':color.blue}),
           ('saturn', {'r':5.823*10**7, 'mass':5.68*10**26, 'pos':vector(1.429*10**12,0,0), 'vel':vector(0,9.69*10**3,0), 'acc':vector(0,0,0), 'color':color.yellow}),
           ('uranus', {'r':2.536*10**7, 'mass':8.68*10**25, 'pos':vector(2.871*10**12,0,0), 'vel':vector(0,6.81*10**3,0), 'acc':vector(0,0,0), 'color':color.cyan}),
           ('neptune', {'r':2.462*10**7, 'mass':1.024*10**26, 'pos':vector(4.498*10**12,0,00), 'vel':vector(0,5.43*10**3,0), 'acc':vector(0,0,0), 'color':color.blue})) # I added color here

#visualize objects
solar_system = {'sun': sphere(pos=sun['pos'], radius=sun['r'], color=sun['color'], make_trail=True)}
for name, planet in planets:
    solar_system[name] = sphere(pos=planet['pos'], radius=planet['r'], color=planet['color'], make_trail=True)

t = 0.0
dt = 150000

while True:
    for name, planet in planets:
        rate(1000000)
        t=t+dt
        r =  sun['pos'] - planet['pos'] # also possible  r = solar_system['sun'].pos - solar_system[name].pos
        Fnet= ((G*sun['mass']*planet['mass'])/(mag(r)**2))* norm(r) #Calculating the force of gravity from the sun
        planet['acc'] = Fnet/planet['mass']
        planet['vel'] += planet['acc']*dt
        planet['pos'] += planet['vel']*dt
        solar_system[name].pos = planet['pos']

In order to preserve the order of the planets, I used tuple of tuples to define the planets, but you can also use OrderedDict from collections.

(Jun-13-2017, 08:48 AM)buran Wrote: [ -> ]Very interesting indeed.
I changed your script a bit to make it more easy to understand and follow

from visual import * # I HATE this
 
myuniv = display(range=vector(7*10**12,7*10**12,7*10**12))
 
G = 6.674*10**-11  #Gravitational constant

# define objects
sun = {'r':9.9*10**10, 'mass':1.989*10**30, 'pos':vector(0,0,0), 'vel':vector(0,0,0), 'acc':vector(0,0,0), 'color':color.yellow}
planets = (('mercury', {'r':2.44*10**6, 'mass':3.285*10**23,'pos':vector(5.791*10**10,0,0), 'vel': vector(0,4.87*10**4,0), 'acc':vector(0,0,0), 'color':color.red}),
           ('venus', {'r':6.052*10**6, 'mass':4.867*10**24, 'pos':vector(1.082*10**11,0,0), 'vel':vector(0,3.502*10**4,0), 'acc':vector(0,0,0), 'color':color.yellow}), # I added color here
           ('earth', {'r':6.371*10**6, 'mass':5.972*10**24, 'pos':vector(1.496*10**11,0,0), 'vel':vector(0,2.98*10**4,0), 'acc':vector(0,0,0), 'color':color.green}),
           ('mars', {'r':3.39*10**6, 'mass':6.39*10**23, 'pos':vector(2.279*10**11,0,0), 'vel':vector(0,2.401*10**4,0), 'acc':vector(0,0,0), 'color':color.red}),
           ('jupiter', {'r':6.99*10**7, 'mass':1.898*10**27, 'pos':vector(7.785*10**11,0,0), 'vel':vector(0,1.307*10**4,0), 'acc':vector(0,0,0), 'color':color.blue}),
           ('saturn', {'r':5.823*10**7, 'mass':5.68*10**26, 'pos':vector(1.429*10**12,0,0), 'vel':vector(0,9.69*10**3,0), 'acc':vector(0,0,0), 'color':color.yellow}),
           ('uranus', {'r':2.536*10**7, 'mass':8.68*10**25, 'pos':vector(2.871*10**12,0,0), 'vel':vector(0,6.81*10**3,0), 'acc':vector(0,0,0), 'color':color.cyan}),
           ('neptune', {'r':2.462*10**7, 'mass':1.024*10**26, 'pos':vector(4.498*10**12,0,00), 'vel':vector(0,5.43*10**3,0), 'acc':vector(0,0,0), 'color':color.blue})) # I added color here

#visualize objects
solar_system = {'sun': sphere(pos=sun['pos'], radius=sun['r'], color=sun['color'], make_trail=True)}
for name, planet in planets:
    solar_system[name] = sphere(pos=planet['pos'], radius=planet['r'], color=planet['color'], make_trail=True)

t = 0.0
dt = 150000

while True:
    for name, planet in planets:
        rate(1000000)
        t=t+dt
        r =  sun['pos'] - planet['pos'] # also possible  r = solar_system['sun'].pos - solar_system[name].pos
        Fnet= ((G*sun['mass']*planet['mass'])/(mag(r)**2))* norm(r) #Calculating the force of gravity from the sun
        planet['acc'] = Fnet/planet['mass']
        planet['vel'] += planet['acc']*dt
        planet['pos'] += planet['vel']*dt
        solar_system[name].pos = planet['pos']

In order to preserve the order of the planets, I used tuple of tuples to define the planets, but you can also use OrderedDict from collections.

That's a much cleaner way to organize the planets info!

I'm still a little bit confused as to how I calculate r for the rest of the planets.

When I try this

rmercury = solar_system['mercury'].pos-solar_system[name].pos

I get the same 'float division by zero' error I had before. I need a way to call the position of every planet except for mercury, which is where I'm stuck.
Also, I've never used tuples before. Would the mass of mercury be something like planet['mercury'].mass?

sorry, maybe I don't understand fully what you try to achieve. My understanding was you want to calculate the gravity impact of the sun over each planet.
In your set-up you had the sun in each of the different lists (mass, posit, etc.), thus you had to exclude it when loop over lists to calculate the gravity impact.

In my code I do just what you do in your code, but in more straightforward (I hope) way. Planets holds only the planets and the solar_system dict holds only the visual representation (sphere objects) of sun and planets. Note that you need to keep position info held in both structures in sync. In more advanced approach, using class, this can be done automatically.

I loop over all planets (line 27) and calculate the impact of the sun over each planet and update the data for planet position in planets tuple (lines 29-34) and in solar_system dict (line 35).
If you elaborate further on your task/model  would suggest changes. Do you want to calculate also the impact between the planets (all or only bigger ones toward smaller ones)?

as to the tuple - it's more or less same as list, but immutable, i.e. you cannot change the tuple once created.

Python 2.7.12 (default, Nov 19 2016, 06:48:10) 
[GCC 5.4.0 20160609] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> l = [1,2,3] # list
>>> t = (1,2,3) # tuple
>>> l[1]
2
>>> t[1]
2
>>> l.insert(0,0) # insert 0 in the beginning of the list
>>> l
[0, 1, 2, 3]
>>> t.insert(0,0) # try to insert 0 in the begining of the tuple
Traceback (most recent call last):
 File "<stdin>", line 1, in <module>
AttributeError: 'tuple' object has no attribute 'insert'

in my code I have a bit more complex structure - tuple of tuples, i.e. each element itself is a tuple with 2 elements - str (the planet name) and dict (the planet data).
Note that while tuple is immutable, you can change the mutable elements, like dict within the tuple:

>>> a= (('a',{'r':1,'m':2}),('b',{'r':3,'m':4})) #tuple of tuples. each element is tuple of string and dict like planets in our case
>>> a[1]
('b', {'r': 3, 'm': 4})
>>> x,y = a[1] # unpack element with index 1
>>> y['m']=5 # change element in the dict
>>> a[1]
('b', {'r': 3, 'm': 5})
>>> a
(('a', {'r': 1, 'm': 2}), ('b', {'r': 3, 'm': 5}))

You can access tuple elements by index (like with list) so print planets[2] would print the tuple element with index 2 (0-based), so that's the tuple for the earth.
using name, planet = planets[2] you can unpack the tuple element with index 2 into its own elements, i.e. name would hold the str and planet would hold the dict.

there should be exactly one list (or object of whatever type) that contains all system objects/bodies. the sun and mercury should not be any different than the others although they will certainly have different mass, position, velocity, colour and name.

(Jun-14-2017, 04:45 AM)Skaperen Wrote: [ -> ]there should be exactly one list (or object of whatever type) that contains all system objects/bodies. the sun and mercury should not be any different than the others although they will certainly have different mass, position, velocity, colour and name.

yes, I was thinking the same

from visual import * # I HATE this

myuniv = display(range=vector(7*10**12,7*10**12,7*10**12))

G = 6.674*10**-11  #Gravitational constant

# define objects
cosmic_objects = (('sun', {'r':9.9*10**10, 'mass':1.989*10**30, 'pos':vector(0,0,0), 'vel':vector(0,0,0), 'acc':vector(0,0,0), 'color':color.yellow}),
           ('mercury', {'r':2.44*10**6, 'mass':3.285*10**23,'pos':vector(5.791*10**10,0,0), 'vel': vector(0,4.87*10**4,0), 'acc':vector(0,0,0), 'color':color.red}),
           ('venus', {'r':6.052*10**6, 'mass':4.867*10**24, 'pos':vector(1.082*10**11,0,0), 'vel':vector(0,3.502*10**4,0), 'acc':vector(0,0,0), 'color':color.yellow}), # I added color here
           ('earth', {'r':6.371*10**6, 'mass':5.972*10**24, 'pos':vector(1.496*10**11,0,0), 'vel':vector(0,2.98*10**4,0), 'acc':vector(0,0,0), 'color':color.green}),
           ('mars', {'r':3.39*10**6, 'mass':6.39*10**23, 'pos':vector(2.279*10**11,0,0), 'vel':vector(0,2.401*10**4,0), 'acc':vector(0,0,0), 'color':color.red}),
           ('jupiter', {'r':6.99*10**7, 'mass':1.898*10**27, 'pos':vector(7.785*10**11,0,0), 'vel':vector(0,1.307*10**4,0), 'acc':vector(0,0,0), 'color':color.blue}),
           ('saturn', {'r':5.823*10**7, 'mass':5.68*10**26, 'pos':vector(1.429*10**12,0,0), 'vel':vector(0,9.69*10**3,0), 'acc':vector(0,0,0), 'color':color.yellow}),
           ('uranus', {'r':2.536*10**7, 'mass':8.68*10**25, 'pos':vector(2.871*10**12,0,0), 'vel':vector(0,6.81*10**3,0), 'acc':vector(0,0,0), 'color':color.cyan}),
           ('neptune', {'r':2.462*10**7, 'mass':1.024*10**26, 'pos':vector(4.498*10**12,0,00), 'vel':vector(0,5.43*10**3,0), 'acc':vector(0,0,0), 'color':color.blue})) # I added color here


#visualize objects
solar_system = {}
for name, obj in cosmic_objects:
    solar_system[name] = sphere(pos=obj['pos'], radius=obj['r'], color=obj['color'], make_trail=True)

t = 0.0
dt = 150000

while True:
    for pla_name, planet in cosmic_objects[1:]: # just the planets
        for obj_name, obj in sorted(cosmic_objects, key=lambda x:x[1]['mass'], reverse=True): # loop ove all objects, sorted by mass descending
            rate(1000000)
            t=t+dt
            if obj_name != pl_name: # check that obj and planet are no the same cosmic object
                r =  obj['pos'] - planet['pos'] # distance between obj and planet
                Fnet= ((G*obj['mass']*planet['mass'])/(mag(r)**2))* norm(r) #Calculating the force of gravity from the obj
                planet['acc'] = Fnet/planet['mass']
                planet['vel'] += planet['acc']*dt
                planet['pos'] += planet['vel']*dt
                solar_system[pl_name].pos = planet['pos']

Sorry for not making the task clear. I need the program to account for the gravitational effects of each body on every other body. Unfortunately I'm going to have to stick with my original format since I'm not familiar enough with tuple yet, and I need this working by tomorrow.

Here's my current version:

from visual import *

myuniv = display(range=vector(7*10**12,7*10**12,7*10**12))

G = 6.674*10**-11  #Gravitational constant

#Masses of each solar body in order starting with the Sun
mass = [1.989*10**30,3.285*10**23, 4.867*10**24, 
           5.972*10**24, 6.39*10**23, 1.898*10**27,
           5.68*10**26, 8.68*10**25, 1.024*10**26]

#Positions of each solar body starting with the Sun           
posit = [vector(0,0,0),vector(5.791*10**10,0,0),vector(1.082*10**11,0,0),
          vector(1.496*10**11,0,0),vector(2.279*10**11,0,0),vector(7.785*10**11,0,0),
          vector(1.429*10**12,0,0),vector(2.871*10**12,0,0),vector(4.498*10**12,0,0)]

#Velocities of each solar body starting with the Sun
vel = [vector(0,0,0),vector(0,4.87*10**4,0),vector(0,3.502*10**4,0),
          vector(0,2.98*10**4,0),vector(0,2.401*10**4,0),vector(0,1.307*10**4,0),
          vector(0,9.69*10**3,0),vector(0,6.81*10**3,0),vector(0,5.43*10**3,0)]

#Acceleration of each solar body starting with the Sun
acc = [vector(0,0,0),vector(0,0,0),vector(0,0,0),
          vector(0,0,0),vector(0,0,0),vector(0,0,0),
          vector(0,0,0),vector(0,0,0),vector(0,0,0)]


t = 0.0
dt = 150000

sun = sphere(pos=posit[0],radius=9.9*10**10,color=color.yellow,make_trail=True)
mercury = sphere(pos=posit[1],radius=2.44*10**6,color=color.red,make_trail=True)
venus = sphere(pos=posit[2],radius=6.052*10**6,make_trail=True)
earth = sphere(pos=posit[3],radius=6.371*10**6,color=color.green,make_trail=True)
mars = sphere(pos=posit[4],radius=3.39*10**6,color=color.red,make_trail=True)
jupiter = sphere(pos=posit[5],radius=6.99*10**7,color=color.blue,make_trail=True)
saturn = sphere(pos=posit[6],radius=5.823*10**7,color=color.yellow,make_trail=True)
uranus = sphere(pos=posit[7],radius=2.536*10**7,color=color.cyan,make_trail=True)
neptune = sphere(pos=posit[8],radius=2.462*10**7,make_trail=True)

while True:
    for i in arange(0,9,1):
        rate(1000000)
        t=t+dt
        rsun = posit[0]-posit[i]
        
        if mag(rsun) == 0:
            continue
        Fsun = ((G*mass[0]*mass[i])/(mag(rsun)**2))* norm(rsun) #Calculating the force of gravity from the sun

        rmerc = posit[1]-posit[i]
        if mag(rmerc)==0:
            continue
        Fmerc = ((G*mass[1]*mass[i])/(mag(rmerc)**2))* norm(rmerc)

        rven = posit[2]-posit[i]
        if mag(rven)==0:
            continue
        Fven = ((G*mass[2]*mass[i])/(mag(rven)**2))* norm(rven)

        rear = posit[3]-posit[i]
        if mag(rear)==0:
            continue
        Fear = ((G*mass[3]*mass[i])/(mag(rear)**2))* norm(rear)

        rmars = posit[4]-posit[i]
        if mag(rmars)==0:
            continue
        Fmars = ((G*mass[4]*mass[i])/(mag(rmars)**2))* norm(rmars)

        rjup = posit[5]-posit[i]
        if mag(rjup)==0:
            continue
        Fjup = ((G*mass[5]*mass[i])/(mag(rjup)**2))* norm(rjup)

        rsat = posit[6]-posit[i]
        if mag(rsat)==0:
            continue
        Fsat = ((G*mass[6]*mass[i])/(mag(rsat)**2))* norm(rsat)

        ruran = posit[7]-posit[i]
        if mag(ruran)==0:
            continue
        Furan = ((G*mass[7]*mass[i])/(mag(ruran)**2))* norm(ruran)

        rnep = posit[8]-posit[i]
        if mag(rnep)==0:
            continue
        Fnep = ((G*mass[8]*mass[i])/(mag(rnep)**2))* norm(rnep)
        
        Fnet = Fsun + Fmerc + Fven + Fear + Fmars + Fjup + Fsat + Furan + Fnep
        acc[i]=Fnet/mass[i]
        vel[i]=vel[i] + acc[i]*dt
        posit[i]=posit[i] + vel[i]*dt

        

        sun.pos=posit[0]
        mercury.pos=posit[1]
        venus.pos=posit[2]
        earth.pos=posit[3]
        mars.pos=posit[4]
        jupiter.pos=posit[5]
        saturn.pos=posit[6]
        uranus.pos=posit[7]
        neptune.pos=posit[8]

I avoided dividing by zero by skipping iterations where r = 0. Now I need to change the range of 'i' for each gravity calculation. For example, when calculating mercury, 'i' needs to be in the range 1-9, rather than 0-9. For venus it needs to be 2-9. Otherwise the forces end up canceling out and the planets don't move.

(Jun-14-2017, 04:45 AM)Skaperen Wrote: [ -> ]there should be exactly one list (or object of whatever type) that contains all system objects/bodies. the sun and mercury should not be any different than the others although they will certainly have different mass, position, velocity, colour and name.

Do I have to have a consolidated list? Or can I use separate lists for mass, position, velocity, and acceleration like I have in the above code? My only issue with changing it at this point is my time restraint. (It's due in the morning. I know I'm a horrible procrastinator...)

(Jun-14-2017, 05:18 AM)JakeWitten Wrote: [ -> ]Do I have to have a consolidated list?  Or can I use separate lists for mass, position, velocity, and acceleration like I have in the above code?  My only issue with changing it at this point is my time restraint. (It's due in the morning. I know I'm a horrible procrastinator...)

yes, you could do it with separate lists, where a common list index is your system body identification. if you already know Fortran (still in common use in physics) then you probably see that as the classic way there.

too bad you don't have the time to do it right. but then, if you did do it right they might not understand it.