Help with modeling Alpha Centauri

Dev Shed Forums Sponsor:

October 23rd, 2012, 06:21 PM

ldownes1

Registered User

Join Date: Oct 2012

Posts: 3

Time spent in forums: 35 m 18 sec
Reputation Power: 0

Help with modeling Alpha Centauri

Hello all! I have attempted to model the three stars and one planet in the system Alpha Centauri for a class, and for some reason they don't want to orbit each other with actual values and I'm at a loss for what to do. Here's what I have so far:

Code:

from visual import *
from random import uniform
from math import sqrt, pow

#list constants
a=1.49E11
G=6.67E-11

#Define objects. I used values for position, radius, mass, and color found online:
starA = sphere( pos= (-1.5E11 ,0, 0) , radius= 8.533785E8, mass= 2.188E30, color =color.yellow)
starB = sphere( pos= (1.5E11 ,0, 0) , radius= 6.016075E8, mass= 1.804E30, color =color.orange)
starC = sphere( pos= (1.94E12, 0, 0) , radius= 9.80655E7, mass= 2.446593E29, color =color.magenta)
planet = sphere(pos = (1.44e11, 0 , 0), radius = 3.5E6, mass = 5.9E24, material = materials.earth)


planet.vel = vector(0 , -22000, 0)
starA.vel = vector(0, 22000, 0)
starB.vel = vector(0, -22000, 0)
starC.vel = vector( 0,  +22000,0)
starA.trail = curve(color = starA.color)
starB.trail = curve(color = starB.color)
starC.trail = curve(color = starC.color)
planet.trail = curve(color = color.white)

# draw coordinate axes
L=2E11
xaxis = curve(pos=[(0,0,0), (L,0,0)], color=color.white)
yaxis = curve(pos=[(0,0,0), (0,L,0)], color=color.white)
zaxis = curve(pos=[(0,0,0), (0,0,L)], color=color.white)

#center of mass
CM = sphere(pos=(0,0,0), radius = 2E12, color=color.green)


#Algorithm
dt = 1.0E8
t=0
dt2=1.0E3

while True:
    r1 = mag(starA.pos - starB.pos)
    r2=mag(starA.pos-starC.pos)
    r3=mag(starB.pos-starC.pos)
    r4=mag(starA.pos - planet.pos)
    r5=mag(starB.pos - planet.pos)
    r6=mag(starC.pos - planet.pos)
    
    starA.vel += ((-G*starB.mass*(starB.pos - starA.pos)/r1**3 - G*starC.mass*(starC.pos - starA.pos)/r2**3 - G*planet.mass*(starB.pos - planet.pos)/r4**3)/(starA.mass))*dt
    starB.vel += ((-G*starA.mass*(starA.pos - starB.pos)/r1**3 - G*starC.mass*(starC.pos - starB.pos)/r3**3 - G*planet.mass*(starB.pos - planet.pos)/r5**3)/(starB.mass))*dt
    starC.vel += ((-G*starB.mass*(starB.pos - starC.pos)/r3**3 - G*starA.mass*(starA.pos - starC.pos)/r2**3 - G*planet.mass*(starC.pos - planet.pos)/r6**3)/(starC.mass))*dt
    planet.vel += ((-G*starA.mass*(starA.pos - planet.pos)/r4**3 - G*starB.mass*(starB.pos - planet.pos)/r5**3 - G*starC.mass*(starC.pos - planet.pos))/r6**3/(planet.mass))*dt2


    starA.pos += starA.vel*dt
    starB.pos += starB.vel*dt
    starC.pos += starC.vel*dt
    planet.pos += planet.vel*dt

    starA.trail.append(pos = starA.pos)
    starB.trail.append(pos = starB.pos)
    starC.trail.append(pos = starC.pos)
    planet.trail.append(pos = planet.pos)

    #M.pos = (starA.mass*starA.pos + starB.mass*starB.pos + starC.mass*starC.pos + planet.mass*planet.pos)/(starA.mass + starB.mass + starC.mass + planet.mass)


    rate(500)

    t += dt

Any help would be appreciated! Thank you!

October 23rd, 2012, 07:58 PM

Lux Perpetua

Contributing User

Dev Shed Intermediate (1500 - 1999 posts)

Join Date: Feb 2004

Location: San Francisco Bay

Posts: 1,936

Time spent in forums: 1 Month 1 Week 2 h 12 m 42 sec
Reputation Power: 1312

I'm not a physicist or astronomer, but: are your initial conditions (position and velocity) realistic? They look made up to me, not that I would know. Another obvious thing to consider is decreasing your timestep. If you get noticeably different results by halving your time step, then it's probably way too big, and you should decrease it until decreasing it more doesn't help. (Keep in mind that nature's time step is infinitesimal, while yours is 1e8 in some unknown unit.)

Edit: Wait, are you using a different time step for your planet? And why does starA's acceleration depend on starB.pos - planet.pos? The more I read, the less sense it makes...

Last edited by Lux Perpetua : October 23rd, 2012 at 08:07 PM.

October 23rd, 2012, 08:43 PM

ldownes1

Registered User

Join Date: Oct 2012

Posts: 3

Time spent in forums: 35 m 18 sec
Reputation Power: 0

No, those numbers are not made up. I found them online in astronomical units (AU) and converted them to meters so that's why they're so huge. I was using a different time step for the planet since it shot off much faster than the other bodies.
I tried making the time step the same for each of them, and changed it to 1e2. They don't shoot as quickly now, but they're still definitely moving away in straight lines, not curves.
I defined star A's velocity based on acceleration due to gravity, so v=G*m*r/(r^3)*dt is what I used to define that. Gravitational acceleration is dependent on distance, so I used starB.pos - planet.pos to define the distance between each one.
Now I'm thinking maybe I'll switch back to AU, so that the distance between each of them aren't such huge numbers and maybe that will make a difference?

October 23rd, 2012, 09:02 PM

b49P23TIvg

Contributing User

Join Date: Aug 2011

Posts: 3,372

Time spent in forums: 1 Month 2 Weeks 3 Days 11 h 39 m 38 sec
Reputation Power: 383

I am a physicist, and I'll spend some time thinking about this, but based on Lux Perpetua's observation then perhaps you ought to update the planet's position many times for each time you change that of the stars?

My improvement is likely to be
"write a system of coupled equations and use a differential equation solver".

__________________
[code]Code tags[/code] are essential for python code!

October 24th, 2012, 05:35 AM

Lux Perpetua

Contributing User

Join Date: Feb 2004

Location: San Francisco Bay

Posts: 1,936

Time spent in forums: 1 Month 1 Week 2 h 12 m 42 sec
Reputation Power: 1312

It's not the numbers that look made up to me, but the starting positions and directions of motion. You have all the stars and planet lined up on a single line and then taking off perpendicularly to that line, and they all have the same velocity up to a factor of +/-1.

October 25th, 2012, 01:35 AM

ldownes1

Registered User

Join Date: Oct 2012

Posts: 3

Time spent in forums: 35 m 18 sec
Reputation Power: 0

I did arbitrarily choose to line them up on the x-axis for simplicity's sake, since their orbits line up, just like the planets do, occasionally. I made them take off linearly to that line because that is the direction of their tangential velocity, while the gravitational force(should) be causing their orbits to loop around.
I found the velocities for A and B online, and I couldn't find it for C, but just set it equal to A and B, which shouldn't make that much of a difference since it's so much further away.
b49P23TIvg, how would I go about writing a system of coupled equations? any suggestions? also how would that help?
and what do you mean by update the planet's position? I thought that's what the line "planet.pos += planet.vel*dt" was doing.

October 25th, 2012, 03:17 PM

Lux Perpetua

Contributing User

Join Date: Feb 2004

Location: San Francisco Bay

Posts: 1,936

Time spent in forums: 1 Month 1 Week 2 h 12 m 42 sec
Reputation Power: 1312

Quote:

Originally Posted by ldownes1

I don't know if it matters in this particular case, but initial conditions certainly can have an effect on the long-term behavior of a system of differential equations.

Quote:

Originally Posted by ldownes1

b49P23TIvg, how would I go about writing a system of coupled equations? any suggestions? also how would that help?
and what do you mean by update the planet's position? I thought that's what the line "planet.pos += planet.vel*dt" was doing.

I really recommend that you think through the mathematics and your approximation a bit. You may not know it, but you are numerically solving a system of differential equations, and you're not doing it right. The method you're implementing is called "Euler's method," and it's not particularly accurate, but it will generally get the job done if your step size is sufficiently small (depending completely on the problem). I'm pointing this out in case you want to go deeper into the type of computation you're doing.

Dusting off what little physics I know:

The force felt by starA due to a planet is (in Newtonian physics, ignoring relativity):
Code:
```
G*starA.mass*planet.mass*(planet.pos - starA.pos)/mag(planet.pos - starA.pos)**3
```
Force is the product of mass times acceleration, so starA's instantaneous acceleration due to the planet is:
Code:
```
G*planet.mass*(planet.pos - starA.pos)/mag(planet.pos - starA.pos)**3
```
Acceleration is the derivative with respect to time of velocity.
Velocity is the derivative with respect to time of position.

Simplifying the system down to two bodies (starA and planet), this gives us the system of differential equations (writing d{X}/dt for derivatives; this is not literally python code):

Code:

d{starA.pos}/dt = starA.vel
d{starA.vel}/dt = G*planet.mass*(planet.pos - starA.pos)/mag(planet.pos - starA.pos)**3
d{planet.pos}/dt = planet.vel
d{planet.vel}/dt = G*starA.mass*(starA.pos - planet.pos)/mag(starA.pos - planet.pos)**3

Now, in Euler's method, we use a finitely small time step dt and use the approximation that d{X}/dt[at time t]= (X[at time t+dt] - X[at time t])/dt. Here I've written the dependence of everything on time explicitly, and you should do the same when you write down equations like this. Using this, we can approximate the above differential equations in terms of computable quantities:

Code:

(starA.pos[at time t+dt] - starA.pos[at time t])/dt == starA.vel # approximate

# rewrite as:

starA.pos[at time t+dt] = starA.pos[at time t] + starA.vel[at time t]*dt # approximate
# do the same for the other 3 equations:
starA.vel[at time t+dt] = starA.vel[at time t] + G*planet.mass*(planet.pos[at time t] - starA.pos[at time t])/mag(planet.pos[at time t] - starA.pos[at time t])**3*dt # approximate
planet.pos[at time t+dt] = planet.pos[at time t] + planet.vel[at time t]*dt # approximate
planet.vel[at time t+dt] = planet.vel[at time t] + G*starA.mass*(starA.pos[at time t] - planet.pos[at time t])/mag(starA.pos[at time t] - planet.pos[at time t])**3*dt # approximate

I hope this has shed some light on what you're doing. I notice a number of differences between my equations and yours; if I've made any errors, I hope the point got across. In particular, I've made the dependence on time explicit everywhere. I hope you can see that replacing `dt' with `dt2' in the right-hand side of the last equation makes zero mathematical sense; you have to replace it on the left-hand side, too!

Code:

planet.vel[at time t+dt2] = planet.vel[at time t] + G*starA.mass*(starA.pos[at time t] - planet.pos[at time t])/mag(starA.pos[at time t] - planet.pos[at time t])**3*dt2 # approximate

...and now you have a problem if you code this up: you have the approximation of planet.vel at time t+dt2 and the approximation of everything else at time t+dt. But now you can't continue the computation to the next step t+2*dt, because you need the values of all the quantities at the same time value. You can theoretically use a smaller time step for planet.vel, but then you have to do multiple time steps for planet.vel for each time step of the other quantities, in order to get planet.vel[at time t+dt]. My advice is not to try that until you get the easy version working.

October 25th, 2012, 04:10 PM

b49P23TIvg

Contributing User

Join Date: Aug 2011

Posts: 3,372

Time spent in forums: 1 Month 2 Weeks 3 Days 11 h 39 m 38 sec
Reputation Power: 383

The reference I found showed said the planet orbits very closely to the larger of those three stars, a good deal closer to the star than Mercury is to our good old Sol.

The two larger stars have an orbit separated by something like that of Sun to Neptune (only not at all circular).

The third star is farther away.

None of this matters too much, since in general you want to solve a 4 body planar gravitational problem. Except for the planet, which may be greatly affected.

dt = 1.0E8 # 3 years
dt2=1.0E3 # 15 minutes

If an exoplanet year is 40 earth days, then you'd need to update its position more often than once every 3 earth years, so it seems.

Hamilton's equations of motion are well suited to numerical ODE solvers since they produce systems of first order equations. And perhaps it's easier to use a cylindrical coordinate system, you'd have angular momentum about the center of mass, and radial momentum. Since momentum rather then velocity would be a variable, kinetic energy is p^2/(2 m) .

Three body gravity problems are chaotic. A tiny change in starting condition leads to a huge difference in positions a million years later. Orbits are hard!

You've got two spatial degrees of freedom per body, x and y, and you've got unknown position and velocity, so you'll have 16 coupled first order equations. None of this should be big news since it's exactly what you solved. A canned ODE solver will have variable time step with error control, and will already have been tested and debugged.

The gls module is a zip file in post 4.

Lux's recent post looks sharp. Here's code using the gsl to solve planar Newtonian orbital mechanics problems.

Code:

# python3
# additional information
# http://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html

# example runs the system for 15 earth years.
# The 3rd object (small star) showed an obviously curved path.

import math
import gsl
from gsl import matrix
from gsl import ode
vector = matrix.vector

#list constants
a=1.49E11
G=6.67E-11

class CelestialObject:

    def __init__(self,**args):
        self.__dict__.update(args)

objects = [CelestialObject(name=n,x=vector((x,0)),v=vector((0,v)),mass=m) for (n,x,v,m,) in zip(
    'ACA ACB ACC ACBb'.split(),
    ( -1.5E11,  1.5E11,    1.94E12,1.44e11),
    (   22000,  -22000,      22000, -22000),
    (2.188E30,1.804E30,2.446593E29, 5.9E24),
   )]

class F(list):

    '''
        # ydot = f(t,y) = dydt(t, y)

        F() supplies a callable of (t,y) that conforms to the spring example below.
        Let the equations be ordered (positions, velocities)
        Positions are   x1, y1,     x2,  y2, ...
        velocities are vx1, vy1,   vx2, vy2, ...
        where the index runs over the celestial objects.
        Since the time rate of change of position is velocity,
        which is computed by the ODE solver, the last half entries of input y
        are the first half of the dydt result vector.
        Use Newton\'s law to compute the accelerations.
        Why can\'t people be nice and use scipy which is somewhat standard?
    '''

    def __len__(self):
        return 4*list.__len__(self)

    def __call__(self,time,Y,f,params):
        L = len(self)
        h = L//2
        y = vector((0,)*L)
        for i in range(L):
            y[i] = Y[i]
        result = [0]*L
        result[:h] = y[h:] # copy velocities as time derivatives of position
        for (i,o,) in enumerate(self):
            oa = 0              # acceleration
            # om = o.mass # do not need this
            ox = y[i*2:(i+1)*2] # dimensionality is hard coded.  Sorry.
            for (j,p,) in enumerate(self):
                if i == j:
                    continue
                px = y[j*2:(j+1)*2]
                r = ox-px
                pm = p.mass
                oa += r*(pm/(r.length**3))
            oa *= G
            result[h+i*2:h+(i+1)*2] = oa
        for i in range(L):                # copy the result onto f
            f[i] = result[i]
        return gsl.gsl_return_code['SUCCESS']

function = F(objects)

de = ode.IV_ODE(function)

earth_year = math.pi*1e7

# set the initial conditions, all positions followed by all velocities
IV = []
for o in objects:
    IV.extend(o.x)
for o in objects:
    IV.extend(o.v)
IV = vector(IV)

print('#time x y '+' '.join(o.name for o in function)+' followed by velocities')
for (t,y,) in de(IV,0,15*earth_year,earth_year/1000):
    a = str(y)
    print(t,' '.join(a[1+a.index('['):a.index(']')].split(',')))   # joy





#class three_springs_two_masses:
#
#    def __init__(self,m1,m2,k1,L1,K2,L2,k3,L3):
#        self.args = (m1,m2,k1,L1,K2,L2,k3,L3)
#
#    def __len__(self):
#        return len('mass1speed mass2speed mass1acceleration mass2acceleration'.split())
#
#    def __call__(self,t,y,f,params):
#        (m1,m2,k1,L1,k2,L2,k3,L3) = self.args
#        f[0] = y[2]
#        f[1] = y[3]
#        f[2] =  -(-k1*(L1-y[0])+k2*(L2-(y[1]-y[0])))
#        f[3] =     k2*(L2-(y[1]-y[0]))-k3*(y[1]-(L1+L2))
#        return gsl.gsl_return_code['SUCCESS']
#
#springs = three_springs_two_masses(*range(1,9))
#the_ODE = gsl.ode.IV_ODE(springs)
#grrrk = len("array('d', [")
#y = (.1,.1,0,0)
#for (t,t_end,) in pairwise(range(0,102,4)):
#    ODEiter = the_ODE(y,t,t_end,1e-6)
#    for (t,y,) in ODEiter:
#        pass
#    #print(t,str(y)[grrrk:-2].replace(',',' '))
#del the_ODE

Help with modeling Alpha Centauri

Developer Shed Advertisers and Affiliates