-

Python_Programming_in_OpenGL_Blank/2body.py

@@ -0,0 +1,179 @@
+# 2Body.py
+# a 2 star system based on Piet Hut
+# and Jun Makino's MSA text with
+# Modifications by Stan Blank for
+# use in Python OpenGL/GLUT
+
+from OpenGL.GL import *
+from OpenGL.GLU import *
+from OpenGL.GLUT import *
+from numpy import *
+import sys
+
+#  Set the width and height of the window
+global width
+global height
+
+#  Initial values
+width = 500
+height = 500
+
+# initial values for position, velocity components, and time increment
+global vx1, vy1, vz1, x1, y1, z1, r2, r3, ax1, ay1, az1, dt
+global vx2, vy2, vz2, x2, y2, z2, ax2, ay2, az2, G
+
+# initial x,y,z positions for both stars
+x1 = 1.0
+y1 = 0.0
+z1 = 0.0
+x2 = -1.0
+y2 = 0.0
+z2 = 0.0
+
+# initial vx,vy,vz velocities for both stars
+vx1 = 0.0
+vy1 = -0.128571428
+vz1 = 0.0
+vx2 = 0.0
+vy2 = 0.3
+vz2 = 0.0
+
+# initial masses for both stars
+m1 = 0.7
+m2 = 0.3
+rad1 = 0.1*m1
+rad2 = 0.1*m2
+
+# arbitrary "Big G" gravitational constant
+G = 1.0
+
+# calculate distance and r**3 denominator for universal gravitation
+r2 = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2)
+r3 = r2*sqrt(r2)
+
+# calculate acceleration components along x,y,z axes
+ax1 = -G*(x1-x2)*m2/r3
+ay1 = -G*(y1-y2)*m2/r3
+az1 = -G*(z1-z2)*m2/r3
+ax2 = -G*(x2-x1)*m1/r3
+ay2 = -G*(y2-y1)*m1/r3
+az2 = -G*(z2-z1)*m1/r3
+
+#This value keeps a smooth orbit on my workstation
+#Smaller values slow down the orbit, higher values speed things up
+dt = 0.001
+
+def init():
+	glClearColor(0.0, 0.0, 0.0, 1.0)
+	glEnable(GL_DEPTH_TEST)
+
+def plotFunc():
+	glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT)
+
+	# plot the first star (m1) position
+	glPushMatrix()
+	glTranslatef(x1,y1,z1)
+	glColor3ub(245, 230, 100)
+	#glColor3ub(245, 150, 30)
+	glutSolidSphere(rad1, 10, 10)
+	glPopMatrix()
+
+	# plot the second star (m2) position
+	glPushMatrix()
+	glTranslatef(x2,y2,z2)
+	glColor3ub(245, 150, 30)
+	#glColor3ub(245, 230, 100)
+	glutSolidSphere(rad2, 10, 10)
+	glPopMatrix()
+
+	# swap the drawing buffers
+	glutSwapBuffers()
+
+def Reshape(  w,  h):
+
+	# To insure we don't have a zero height
+	if h==0:
+		h = 1
+
+	#  Fill the entire graphics window!
+	glViewport(0, 0, w, h)
+
+	#  Set the projection matrix... our "view"
+	glMatrixMode(GL_PROJECTION)
+	glLoadIdentity()
+
+	gluPerspective(45.0, 1.0, 1.0, 1000.0)	
+	gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0)
+
+	#  Set the matrix for the object we are drawing
+	glMatrixMode(GL_MODELVIEW)
+	glLoadIdentity()
+
+def keyboard(key, x, y):
+	#  Allows us to quit by pressing 'Esc' or 'q'
+	if key == chr(27):
+		sys.exit()
+	if key == "q":
+		sys.exit()
+
+def orbits():
+	global vx1, vy1, vz1, x1, y1, z1, r2, r3, ax1, ay1, az1
+	global vx2, vy2, vz2, x2, y2, z2, ax2, ay2, az2
+
+	# calculate front half of velocity vector components
+	vx1 += 0.5*ax1*dt
+	vy1 += 0.5*ay1*dt
+	vz1 += 0.5*az1*dt
+	vx2 += 0.5*ax2*dt
+	vy2 += 0.5*ay2*dt
+	vz2 += 0.5*az2*dt
+
+	# calculate x,y,z positions for both stars
+	x1 += vx1*dt
+	y1 += vy1*dt
+	z1 += vz1*dt
+	x2 += vx2*dt
+	y2 += vy2*dt
+	z2 += vz2*dt
+
+	# calculate the new r**3 denominator for each star position
+	r2 = (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2)
+	r3 = r2*sqrt(r2)
+
+	# calculate the new acceleration components
+	ax1 = -G*(x1-x2)*m2/r3
+	ay1 = -G*(y1-y2)*m2/r3
+	az1 = -G*(z1-z2)*m2/r3
+	ax2 = -G*(x2-x1)*m1/r3
+	ay2 = -G*(y2-y1)*m1/r3
+	az2 = -G*(z2-z1)*m1/r3
+
+	# calculate the back half velocity components
+	vx1 += 0.5*ax1*dt
+	vy1 += 0.5*ay1*dt
+	vz1 += 0.5*az1*dt
+	vx2 += 0.5*ax2*dt
+	vy2 += 0.5*ay2*dt
+	vz2 += 0.5*az2*dt
+
+	#send calculated x,y,z star positions to the display
+	glutPostRedisplay()
+
+def main():
+	global width
+	global height
+
+	glutInit(sys.argv)
+	glutInitDisplayMode(GLUT_RGB|GLUT_DOUBLE)
+	glutInitWindowPosition(100,100)
+	glutInitWindowSize(width,height)
+	glutCreateWindow("2 Body Problem")
+	glutReshapeFunc(Reshape)
+	glutDisplayFunc(plotFunc)
+	glutKeyboardFunc(keyboard)
+	glutIdleFunc(orbits)
+
+	init()	
+	glutMainLoop()
+
+main()