Thursday, May 20, 2010

[from James] polar 1-body

The polar 2-body (or really 1-body orbiting the central mass) has been working for a few days now; I won't post the whole code, which is long due to the RKF45 integrator I've been using, but here is the relevant physical modelling part of the code (in python):

# 2-body problem in polar co-ordinates. Assume that G*M = 1 and m = 1.
# "t" is the time, "kin" is a list/array of current orbital values passed
# along from the integrator function.
def g(t,kin):
GM = 1.0
r = kin[0] # radial position
theta = kin[1] # angular position
v = kin[2] # radial velocity
omega = kin[3] # angular velocity
E = kin[4] # total energy
L = kin[5] # angular momentum
dr = v
dtheta = omega
dVr= r*omega*omega - GM/(r*r)
domega = -2.0*omega*v/r # careful not to confuse omega and domega when reading below
dE = GM*v/(r*r) + v*dVr + r*v*omega*omega + r*r*omega*domega
dL = 2*r*omega*dr + r*r*domega

# send computed d/dt values back to the integrator
return np.array((dr,dtheta,dVr,domega,dE,dL))


# Initial conditions. Position in (r, theta) and eccentricity will be
# randomly generated. V (radial velocity) will always initially be zero, ie.
# the particle will always begin at apocenter or pericenter. Other values
# are then calculated from r, theta, and eccentricity.
GM = 1.0
X,Y,Eccent = rand.random(),rand.random(),rand.random()
R = sqrt(X*X + Y*Y)
Theta = atan(Y/X) # radians
V = 0.0
Axis = (R*(1.0 + Eccent*cos(Theta)))/(1.0-Eccent*Eccent) # semi-major axis
Omega = (sqrt(GM*Axis*(1.0-Eccent*Eccent)))/(R*R)
# rad/time
Energy = -0.5*GM/Axis
AngMomentum = R*R*Omega


I'll send the whole file out via e-mail. If you have python you can run it as an executable, in which case it should display a polar plot of the particle's orbit, and also print out a table of data displaying time, r, theta, E and L for each timestep taken.

I know most of you are working in C/C++, but the actual formulas should work just as well regardless of language.

As for the N-body: the random generator of constant number density particle disk is working fine, and I think I have the polar-coordinate physics for the particle-particle gravitational interactions worked out, but I'm running into a wall on the programming side of things (which makes it hard to test the physics!). I hope to have it working before the weekend.

No comments:

Post a Comment