#This is a shell program implementing the basic impulse-momentum analysis
#for a system of two blobs.  One blob is stationary.  The other moves around it
#under the influence of gravity

#setup
from __future__ import division
from visual import *

#constants
G = 6.67e-11

#create objects & initial conditions
#OBJECTS:
movingblob = sphere(pos = vector(1e7,0,0), radius = 1e4, color=color.white)
    #movingblob is a sphere of radius 1000 positioned at the x = 1,000,000 point
stationaryblob = sphere(pos = vector(0,0,0), radius = 8e4, color = color.red)
    #stationaryblob is a sphere of radius 8000 positioned at the origin

#INITIAL CONDITIONS
movingblob.mass = 1e20
movingblob.velocity = vector(0,200,0)
movingblob.p = movingblob.mass * movingblob.velocity
	#This assignes a starting mass, velocity, and momentum (p) to
	# the blob.  Note that this blob is moving diagonally
	# (theta = 53 degrees, magnitude = 50 units)
trail = curve(pos=movingblob.pos)
        # This adds a curve that will trace the path of our moving blob

stationaryblob.mass = 1e22
stationaryblob.velocity = 0
stationaryblob.p = stationaryblob.mass + stationaryblob.velocity
	#This assignes a starting mass, velocity, and momentum (p) to
	# the other blob.  Note that this blob is not moving

radius = movingblob.pos - stationaryblob.pos
rmag = sqrt(radius.x**2 + radius.y**2 + radius.z**2)
rhat = radius/rmag
        #radius is a displacement vector pointing from the stationary
        # to the moving blob.  Rmag is its magnitude, rhat is a
        # unit vector in its direction

fgravmag = G*stationaryblob.mass*movingblob.mass/rmag**2
fgrav = fgravmag*-rhat
	#fgravmag is the magnitude of the gravitational force between the two
        # blobs.  fgrav is the the magnitude multiplied by the direction of the
        # radius vector between the blobs.  fgrav points from stationary to moving,
        # -fgrav would point from moving to stationary

t = 0
dt = 10
	#t is the initial time, dt is the size of the time step.

#calculations
while t<1000000:	#loop until 1,000,000 seconds have passed
    rate(15000)   #This limits program to 50 loops per second, each of which
                # takes 10 s (dt), so it approximates real time
    #Update the blob's position and momentum first
    movingblob.pos = movingblob.pos + (movingblob.p/movingblob.mass)*dt
    trail.append(pos=movingblob.pos, retain = 10000)
	    #current position = prev. position + velocity * deltat
            #trail.append adds the current position of the moving blob
            # to the trail
    movingblob.p = movingblob.p + fgrav*dt
	    #current momentum = prev. momentum + change in momentum
	    #	(change in momentum = force * deltat)

    #Next, update the radius vector based on the moving blob's new position
    radius = movingblob.pos - stationaryblob.pos
    rmag = sqrt(radius.x**2 + radius.y**2 + radius.z**2)
    rhat = radius/rmag

    #Finally, update the gravity force vector
    fgravmag = G*stationaryblob.mass*movingblob.mass/rmag**2
    fgrav = fgravmag*-rhat

    t = t + dt
	    #keep track of the current time