Arrays in VPython

We need to update our code to handle all of the 'planets' in an array. The stars.py example below does this.

===

from visual import *
from time import clock
from random import random

# Stars interacting gravitationally
# Program uses numpy arrays for high speed computations

win=600

Nstars = 20 # change this to have more or fewer stars

G = 6.7e-11 # Universal gravitational constant

# Typical values
Msun = 2E30
Rsun = 2E9
Rtrail = 2e8
L = 4e10
vsun = 0.8*sqrt(G*Msun/Rsun)

scene = display(title="Stars", width=win, height=win,
range=2*L, forward=(-1,-1,-1))

xaxis = curve(pos=[(0,0,0), (L,0,0)], color=(0.5,0.5,0.5))
yaxis = curve(pos=[(0,0,0), (0,L,0)], color=(0.5,0.5,0.5))
zaxis = curve(pos=[(0,0,0), (0,0,L)], color=(0.5,0.5,0.5))

Stars = []
colors = [color.red, color.green, color.blue,
color.yellow, color.cyan, color.magenta]
poslist = []
plist = []
mlist = []
rlist = []

for i in range(Nstars):
x = -L+2*L*random()
y = -L+2*L*random()
z = -L+2*L*random()
r = Rsun/2+Rsun*random()
Stars = Stars+[sphere(pos=(x,y,z), radius=r, color=colors[i % 6])]
Stars[-1].trail = curve(pos=[Stars[-1].pos], color=colors[i % 6], radius=Rtrail)
mass = Msun*r**3/Rsun**3
px = mass*(-vsun+2*vsun*random())
py = mass*(-vsun+2*vsun*random())
pz = mass*(-vsun+2*vsun*random())
poslist.append((x,y,z))
plist.append((px,py,pz))
mlist.append(mass)
rlist.append(r)

pos = array(poslist)
p = array(plist)
m = array(mlist)
m.shape = (Nstars,1) # Numeric Python: (1 by Nstars) vs. (Nstars by 1)
radius = array(rlist)

vcm = sum(p)/sum(m) # velocity of center of mass
p = p-m*vcm # make total initial momentum equal zero

dt = 1000.0
Nsteps = 0
pos = pos+(p/m)*(dt/2.) # initial half-step
time = clock()
Nhits = 0

while 1:
rate(100)

# Compute all forces on all stars
try: # numpy
r = pos-pos[:,newaxis] # all pairs of star-to-star vectors
for n in range(Nstars):
r[n,n] = 1e6 # otherwise the self-forces are infinite
rmag = sqrt(add.reduce(r*r,-1)) # star-to-star scalar distances
hit = less_equal(rmag,radius+radius[:,newaxis])-identity(Nstars)
hitlist = sort(nonzero(hit.flat)[0]).tolist() # 1,2 encoded as 1*Nstars+2
F = G*m*m[:,newaxis]*r/rmag[:,:,newaxis]**3 # all force pairs
except: # old Numeric
r = pos-pos[:,NewAxis] # all pairs of star-to-star vectors
for n in range(Nstars):
r[n,n] = 1e6 # otherwise the self-forces are infinite
rmag = sqrt(add.reduce(r*r,-1)) # star-to-star scalar distances
hit = less_equal(rmag,radius+radius[:,NewAxis])-identity(Nstars)
hitlist = sort(nonzero(hit.flat)) # 1,2 encoded as 1*Nstars+2
F = G*m*m[:,NewAxis]*r/rmag[:,:,NewAxis]**3 # all force pairs

for n in range(Nstars):
F[n,n] = 0 # no self-forces
p = p+sum(F,1)*dt

# Having updated all momenta, now update all positions
pos = pos+(p/m)*dt

# Update positions of display objects; add trail
for i in range(Nstars):
Stars[i].pos = pos[i]
if Nsteps % 10 == 0:
Stars[i].trail.append(pos=pos[i])

# If any collisions took place, merge those stars
for ij in hitlist:
i, j = divmod(ij,Nstars) # decode star pair
if not Stars[i].visible: continue
if not Stars[j].visible: continue
# m[i] is a one-element list, e.g. [6e30]
# m[i,0] is an ordinary number, e.g. 6e30
newpos = (pos[i]*m[i,0]+pos[j]*m[j,0])/(m[i,0]+m[j,0])
newmass = m[i,0]+m[j,0]
newp = p[i]+p[j]
newradius = Rsun*((newmass/Msun)**(1./3.))
iset, jset = i, j
if radius[j] > radius[i]:
iset, jset = j, i
Stars[iset].radius = newradius
m[iset,0] = newmass
pos[iset] = newpos
p[iset] = newp
Stars[jset].trail.visible = 0
Stars[jset].visible = 0
p[jset] = vector(0,0,0)
m[jset,0] = Msun*1E-30 # give it a tiny mass
Nhits = Nhits+1
pos[jset] = (10.*L*Nhits, 0, 0) # put it far away

Nsteps += 1

Vector Table

We found the vector table ephemeris option during our meeting today. Looks like it will solve the remaining issues with our code...

Velocity direction

This paper could help us with our velocity vector direction issue.

http://www.adass.org/adass/proceedings/adass94/aket.html

Planet Positions from JPL

We found a good paper by E. M. Standish at JPL on Keplerian Elements for Approximate Positions of the Major Planets. This should be useful for initializing our program. However, we still don't know how to calculate the initial velocity vector if a planet is not at aphelion or perihelion. (Does this need to be found computationally?)

Adding ELON and ELAT to the code...

Matt and I tried adding the ELON and ELAT values from NASA to the Python code.

Crider3.py

Screen of MaximDL setup for astrometry

Plotting the Planets

Here is a NASA site for predicting the positions of planets on a given date.

http://cohoweb.gsfc.nasa.gov/helios/planet.html


The accuracy is pretty low. We should try to find something more accurate, I think.

New NEO Coming By Soon

We should try to get an image of this!

http://science.slashdot.org/article.pl?sid=09/03/01/208251&from=rss

Matt - Grade Update

Completed Tasks - (Grade):

1. Make list of top ten candidate asteroids and why they're interesting/post on blog - (D-)
2. Take one practice observation of an asteroid/try to get pictures of 1999 AQ10/try every day for a week if necessary - (D)

Simulating the Earth's Orbit

Here's the rough template that will serve as the basis for our Earth-Sun-asteroid orbital computations:

Crider2.py

"Bright" vs "Faint" Asteroids in MPC

The Minor Planet Center has lists of recovery opportunities for bright (V<21) and faint (21custom list.

Spring 2009: Asteroid List

Here is the list I put together of the asteroids to keep an eye on this semester.


List of Interesting Asteroids:
- Toutatis (4179) –
o On September 29, 2004, Toutatis will passed by Earth at a range of four times the distance between the Earth and the Moon, the closest approach of any known asteroid or comet between then and 2060. One consequence of the asteroid's frequent close approaches to Earth is that its trajectory more than several centuries from now cannot be predicted accurately. In fact, of all the Earth-crossing asteroids, the orbit of Toutatis is thought to be one of the most chaotic.
- 1999 AQ10 –
o Will pass 0.0117 AU (4.6 LD) from Earth on February 18th, 2009. Its very close approach to Earth should alter its orbit and make it something interesting to look at/calculate.
- 1994 CC –
o Will pass 0.0170 AU (6.6 LD) from Earth on June 10th, 2009. Its very close approach to Earth should alter its orbit and make it something interesting to look at/calculate.
- Apophis (99942) –
o On Friday, April 13, 2029, Apophis will pass Earth within the orbits of geosynchronous communication satellites. Although its chance of impact has been reduced to a negligible probability, there probably aren’t many more asteroids that more orbital data exists for. This may make it interesting to observe and compare my data with.
- 2007 TU24 –
o Passed closer to Earth last year (1.4 LD) than any asteroid is known to pass Earth until 2027. Its orbit could be of interest because some believe it to be a “rubble pile” asteroid. It could be cool to see how its composition affects its orbit and compares to that of other asteroids.
- 2001 FE90
o Will pass 0.0180 AU (7.0 LD) from Earth on June 28th, 2009. It could be as big as 650 m in diameter.
- 2009 CV
o Will pass 0.0124 AU (64.8 LD) from Earth on February 23th, 2009. It is passing very close to us but is only estimated to be about 90 m in diameter. It may be interesting to look at the orbit of an asteroids of its size, though.
- 1998 OR2
o The largest asteroid that is making a close approach in the near future. It is about 4.4 km in diameter and will miss us by 70 LD on March 12.
- 2003 QO104
o Another very large asteroid (4.2 km diameter( like 1998 OR2, that will miss us by about half the distance (37 LD) on June 9.
- 161989 Cacus
o An asteroid that’s been discovered for some time and should have precise orbital information. It is 2.2 km in diameter, and will miss Earth by 71 LD on March 7.