#----------------------------------------------------------
#OBSOLETE SCRIPT. The function in this script, or their
#equivalents, are all included in solar.py
#----------------------------------------------------------
#
# Computes the orbit of a planet about a point mass,
# and uses the orbit to compute the solar flux as
# a function of time and latitude. 
# The units are such that G*M = 1
#

import rk 
import math
from solar_flux import *

def f(t,X):
	x = X[0]
	y = X[1]
	u = X[2]
	v = X[3]
	r2 = x*x+y*y
	accelx = - x/(r2*math.sqrt(r2))
	accely = -y/(r2*math.sqrt(r2))
	return [u,v,accelx,accely]


#
# The orbit is initialized at perihelion.
# We get aphelion in terms of eccentricity,
# then figure which initial velocity gives us
# the desired aphelion
#
#      Choose your eccentricity here
e = 0.01
#
r_ap = (1.+2.*e)/(1.-e)
xstart = 1.
vstart = math.sqrt(2./(1.+ (1./r_ap))) 

#First do an integration to get the period
dt = .01
f_i = rk.integrator(f,0.,[xstart,0.,0.,vstart],dt)
ylast =1.e30
y = 1.e30
while not ((ylast < 0) & (y > 0)):
    ylast = y
    result = f_i.next()
    time = result[0]
    x = (result[1])[0]
    y = (result[1])[1]
    if time > 10000.:
	    break
T = (time - dt) + dt*ylast/(y-ylast)
print "Period of orbit = ",T

#Now do a run with results saved to file.
lat = math.pi/4.
obliquity = math.pi/8.
# The phase parameter characterizes the precession
# of the perihelion relative to the tilt season.
# Perihelion is a t = 0, so:
#       phase = 0 implies perihelion occurs at Northern summer
#                         aphelion at Northern winter
#       phase = pi/2 implies perihelion at equinox
#       phase = pi implies perihelion at Northern winter
#                          aphelion at Northern summer
phase = pi

N = 360
n_sub = 10
dt = T/float(n_sub*N)
f_i = rk.integrator(f,0.,[xstart,0.,0.,vstart],dt)

outfile = open("orbit.out","w")
fluxfile = open("flux.out","w")
for i in range(n_sub*N):
    ylast = y
    result = f_i.next()
    time = result[0]
    x = (result[1])[0]
    y = (result[1])[1]
    season = math.atan2(y,x)
    if season < 0.:
	    season += 2.*math.pi
    season += phase
    r = math.sqrt(x*x + y*y)
    L = 1./(r*r)
    if i%n_sub == 0:
	    fluxfile.write("%f %f %f\n"%(time/T,
					 L*solar(lat,season,obliquity),
					 L*solar(-lat,season,obliquity)))
	    outfile.write("%f %f %f %f\n"%(time,x,y,r))
outfile.close()
fluxfile.close()