Asteroid orbit visualizer with Python and PyOpenGL for dataset from ftp://ftp.lowell.edu/pub/elgb/astorb.html

asteroids-glut.py

#!/usr/bin/python

"""
asteroids-glut.py

Loads astorb.dat file and plot the orbit with PyOpenGL.
"""
from __future__ import print_function
import sys, numbers, gc
from math import *
from random import *
from vec3 import *
from datetime import datetime, timedelta

try:
  from OpenGL.GLUT import *
  from OpenGL.GL import *
  from OpenGL.GLU import *
except:
  print('''
ERROR: PyOpenGL not installed properly.  
        ''')


class Asteroid(object):
	num = 0
	name = ""
	hmag = 1
	argument_of_perihelion = 0.
	ascending_node = 0
	inclination = 0
	eccentricity = 0
	semimajor_axis = 0
	epoch = datetime.now()
	mean_anomaly = 0

	def __init__(self,line):
		self.num = int(line[:6])
		self.name = line[7:25]
		self.hmag = float(line[42:47])
		self.argument_of_perihelion = float(line[126:136])
		self.ascending_node = float(line[137:147])
		self.inclination = float(line[149:157])
		self.eccentricity = float(line[159:168])
		self.semimajor_axis = float(line[169:181])
		self.epoch = datetime(int(line[106:110]), int(line[110:112]), int(line[112:114]))
		self.mean_anomaly = float(line[115:125]) * 2 * pi / 360

	def get_eccentric_anomaly(self, time=datetime.now()):
		""" Calculates eccentric anomaly from mean anomaly in first order approximation
		see http://en.wikipedia.org/wiki/Eccentric_anomaly """
		GMsun = 1.327124400e11 # Product of gravitational constant (G) and Sun's mass (Msun)
		AU = 149597871 # Astronomical unit in kilometers
		td = time - self.epoch
		period = 2 * pi * sqrt((self.semimajor_axis * AU) ** 3 / GMsun)
		now_anomaly = self.mean_anomaly + td.total_seconds() * 2 * pi / period
		return now_anomaly + self.eccentricity * sin(now_anomaly)

if len(sys.argv) < 2:
	print("usage: " + sys.argv[0] + " path")
	quit(0)

f = open(sys.argv[1])

asteroids = []
maxdist = 0

for line in f:
	if 1000 < len(asteroids):
		break
	try:
		a = Asteroid(line)
		asteroids.append(a)
		# Print first 100 asteroids' names to see if the file is correctly loaded
		if len(asteroids) < 100:
			print(a.num, a.name, a.epoch, a.mean_anomaly)
		if maxdist < a.semimajor_axis:
			maxdist = a.semimajor_axis * (1. + a.eccentricity)
	except e:
		print("Error on parsing line ", line)

f.close()

def init():
	global quadratic
	glClearColor(0.0, 0.0, 0.0, 0.0)
	glClearDepth(1.0)
	glShadeModel(GL_SMOOTH)
	quadratic = gluNewQuadric()
	gluQuadricNormals(quadratic, GLU_SMOOTH)
	gluQuadricTexture(quadratic, GL_TRUE)
	glEnable(GL_CULL_FACE)
	glEnable(GL_DEPTH_TEST)

dist = maxdist * 1.2
phi = 30
theta = 30
antialias = False
showorbit = False
simtime = datetime.now() # Simulated time

def display():
	global simtime
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

	glLoadIdentity()
	# viewing transformation
	lookpos = vec3(cos(phi * pi / 180) * cos(theta * pi / 180), sin(theta * pi / 180), sin(phi * pi / 180) * cos(theta * pi / 180)) * dist
	gluLookAt(lookpos[0], lookpos[1], lookpos[2], 0.0, 0.0, 0.0, 0.0, 1.0, 0.0)

	glBegin(GL_LINES)
	glColor3f (1.0, 0.0, 0.0)
	glVertex3d(-1000, 0, 0)
	glVertex3d( 1000, 0, 0)
	glColor3f (0.0, 1.0, 0.0)
	glVertex3d(0, -1000, 0)
	glVertex3d(0,  1000, 0)
	glColor3f (0.0, 0.0, 1.0)
	glVertex3d(0, 0, -1000)
	glVertex3d(0, 0,  1000)
	glEnd()

	glPushAttrib(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_CURRENT_BIT | GL_ENABLE_BIT | GL_POLYGON_BIT)
	glDepthMask(GL_FALSE)
	if antialias:
		glEnable(GL_LINE_SMOOTH)
		glHint(GL_LINE_SMOOTH_HINT, GL_NICEST)
		glLineWidth(1)
	glEnable(GL_BLEND)
	glBlendFunc(GL_SRC_ALPHA, GL_ONE)

	for a in asteroids:
		# Semi-minor axis
		smia = a.semimajor_axis * (sqrt(1. - a.eccentricity * a.eccentricity) if a.eccentricity != 0. else 1.)

		glPushMatrix()
		glRotated(a.ascending_node, 0, 1, 0)
		glRotated(a.inclination, 1, 0, 0)
		glRotated(a.argument_of_perihelion, 0, 1, 0)
		glScaled(smia, 1, a.semimajor_axis)
		glTranslated(0, 0, a.eccentricity)

		e = a.get_eccentric_anomaly(simtime)
		glColor4f(1.0, 1.0, 1.0, 1)
		glBegin(GL_POINTS)
		glVertex3d(sin(e), 0, cos(e))
		glEnd()

		if showorbit:
			glBegin(GL_LINE_LOOP)
			glColor4f(1.0, 1.0, 1.0, 0.1)
			ellipse_points = 32
			for j in range(ellipse_points):
				jtheta = j * 2 * pi / ellipse_points
				#print sin(jtheta) * p.semimajor_axis, cos(jtheta) * p.semimajor_axis
				glVertex3d(sin(jtheta), 0, cos(jtheta))
			glEnd()
		glPopMatrix()

	glPopAttrib()

	glPushAttrib(GL_LIGHTING_BIT)
	glEnable(GL_LIGHTING)
	glEnable(GL_LIGHT0)
	glLightfv(GL_LIGHT0, GL_POSITION, [1., 2., 1., 0.]);
	glLightfv(GL_LIGHT0, GL_AMBIENT, [0.2, 0.2, 0.2, 1.]);
	glLightfv(GL_LIGHT0, GL_DIFFUSE, [1., 1., 1., 0]);
	glLightfv(GL_LIGHT0, GL_SPECULAR, [0, 0, 1, 0]);
	# Sun's radius measured in AUs
	Rsun = 695800. / 149597871.
	gluSphere(quadratic, Rsun, 20, 20)
	glPopAttrib()

	glFlush()

	glutPostRedisplay()
	print(gc.get_count(), gc.garbage, simtime)

	simtime += timedelta(days=1)

mousestate = False
mousepos = [0,0]

def mouse(button, state, x, y):
	if state:
		mousestate = True
	else:
		mousestate = False
	mousepos[0] = x
	mousepos[1] = y

def motion(x, y):
	global phi, theta
	theta -= mousepos[1] - y
	if theta < -90:
		theta = -90
	if 90 < theta:
		theta = 90
	mousepos[1] = y
	phi -= mousepos[0] - x
	mousepos[0] = x

def keyboard(key, x, y):
	global dist, antialias, showorbit
	if key == '+':
		dist /= 1.1
	if key == '-':
		dist *= 1.1
	if key == chr(27):
		sys.exit(0)
	if key == 'a':
		antialias = not antialias
	if key == 'o':
		showorbit = not showorbit

def reshape (w, h):
	glViewport(0, 0, w, h)
	glMatrixMode(GL_PROJECTION)
	glLoadIdentity()
	glFrustum(-0.1, .1, -.1, .1, 0.1, 5000.0)
	glMatrixMode(GL_MODELVIEW)

glutInit(sys.argv)
glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB)
glutInitWindowSize(500, 500)
glutInitWindowPosition(100, 100)
glutCreateWindow('Asteroid Orbit Plotter')
init()
glutDisplayFunc(display)
glutReshapeFunc(reshape)
glutMouseFunc(mouse)
glutMotionFunc(motion)
glutKeyboardFunc(keyboard)
glutMainLoop()

vec3.py

import numbers
from math import *

class vec3(object):
	""" Basic 3-D vector implementation """
	x = 0
	y = 0
	z = 0
	def __init__(self, x=0, y=0, z=0):
		self.x = x
		self.y = y
		self.z = z
	def __neg__(self):
		return vec3(-self.x, -self.y, -self.z)
	def __add__(self, o):
		return vec3(self.x + o.x, self.y + o.y, self.z + o.z)
	def __sub__(self, o):
		return vec3(self.x - o.x, self.y - o.y, self.z - o.z)
	def __mul__(self, s):
		if not isinstance(s, numbers.Number):
			return NotImplemented
		return vec3(self.x * s, self.y * s, self.z * s)
	def __div__(self, s):
		if not isinstance(s, numbers.Number):
			return NotImplemented
		return vec3(self.x / s, self.y / s, self.z / s)
	def __repr__(self):
		return "(" + str(self.x) + "," + str(self.y) + "," + str(self.z) + ")"

	def __getitem__(self,key):
		if key == 0:
			return self.x
		elif key == 1:
			return self.y
		elif key == 2:
			return self.z
		else:
			raise RangeError

	def __setitem__(self,key,value):
		if key == 0:
			self.x = value
		elif key == 1:
			self.y = value
		elif key == 2:
			self.z = value

	def slen(self):
		return self.x ** 2 + self.y ** 2 + self.z ** 2
	def len(self):
		return sqrt(self.slen())