from visual import *
from ODlib import *
scene.background = color.white
scene.fullscreen =0 

fieldValues = GetElements()
G = 6.674e-11
labels_on = False
###########################################
# Create the sun and setup the solar system
###########################################
sun = CreateSun()
#########################################
# Create the Earth and establish its initial conditions
#########################################
earth = CreateEarth(sun.mass,G)
JD_stop = float(fieldValues[6])
#############################################################
#Move earth to create earth sun vector on appropriate date. Based on JD_Stop
#############################################################
earth = MoveEarth(earth,sun.mass,sun.pos,JD_stop,G)
JD = JD_stop
print "Earth Sun Vector at this point: ",JD,(earth.pos-sun.pos)
##################################################################
# Make the asteroid in its orbit according to the orbital elements
##################################################################
asteroid = CreateAsteroid(fieldValues,G,sun.mass,sun.pos)
EA_vect = arrow(pos = earth.pos, axis = (asteroid.pos - earth.pos),shaftwidth = 2e9)
print "Earth Asteroid Vector at this point: ",JD,(asteroid.pos-earth.pos)
print "Sun Asteroid Vector at this point: ",JD,(asteroid.pos-sun.pos)
asteroid.rho = asteroid.pos - earth.pos
####
# Now find the RA and DEC from the Earth Asteroid Vector
####
# First rotate the sun asteroid vector into equatorial coordinates:
print "rho: ",asteroid.rho
asteroid.rho = rotate(asteroid.rho,angle = (23.439281*pi/180),axis=(1,0,0))
print "rho equatorial: ",asteroid.rho
# Find rhohat
asteroid.rhohat = norm(asteroid.rho)
DEC = asin((asteroid.rhohat.z))
RAcos = acos(asteroid.rhohat.x/cos(DEC))
RAsin = asin(asteroid.rhohat.y/cos(DEC))
print RAcos, RAsin
#Check Angle Ambiguity
if RAsin < 0:
    RAcos = 2*pi - RAcos 
RA = RAcos
print "RA: ",RA*180/pi,"DEC: ",DEC*180/pi
RA = RA*12/pi
DEC = DEC * 180/pi

print "RA (H:MM:SS) ",int(RA),int((RA-int(RA))*60),\
      (((RA-int(RA))*60)-int((RA-int(RA))*60))*60
print "DEC (D:MM:SS) ",int(DEC),int((DEC-int(DEC))*60),\
      (((DEC-int(DEC))*60)-int((DEC-int(DEC))*60))*60


scene.autoscale = 0
#Approximately conserve momentum in the system
sun.p = -earth.p
##################################################################
# Animate the orbit of the asteroid
##################################################################
#Animate_Orbit(earth,asteroid,sun,EA_vect,JD,G,labels_on)