#!/usr/bin/env python2.5
"""
pendulum.py

The program assumes the form

theta_dot = omega
omega_dot = - (g/L) sin(theta)

for the equation of motion of a pendulum.  We integrate
the equation using a forward Euler scheme.

We estimate the period by doubling the time between
successive crossings of theta == 0.  To find a crossing,
we look for where theta and theta_old have opposite signs:

theta_k * theta_(k-1) < 0.

We then estimate the period using the associated times:

T_k = 2*(time_k - time_(k-1)).

The average period is

T_avg = (1/Ntot) (T_1 + T_2 + ... + T_Ntot),

with Ntot the total number of T evaluations.  The error bar
is given by error = s/sqrt(Ntot), with

s = sqrt( (1/(Ntot - 1) * ( (T_1 - Tavg)^2 + ... + (T_Ntot - T_a

We also check the energy per unit mass per unit length
squared at the initial and final times:

e/mL^2 = 0.5 theta_dot^2 - g/L * cos(theta) .
"""

import sys, math, string
import pylab

# how to use the program
usage = 'Usage: %s [options] [outfile]' % sys.argv[0]

# default values for options

# read options from command line
while len(sys.argv) > 1:
    option = sys.argv[1]
    del sys.argv[1]
    if (option == '-o'):                             # output file name
        outfilename = sys.argv[1]
        del sys.argv[1]
    elif (option == '-th0') or (option == '-theta'): # initial angle in degrees
        th0 = float(sys.argv[1]) * pi/180
        del sys.argv[1]
    elif (option == '-om0') or (option == '-omega'): # initial angular speed
        om0 = float(sys.argv[1])
        del sys.argv[1]
    elif (option == '-gl'):                          # g/L
        gl = float(sys.argv[1])
        del sys.argv[1]
    elif (option == '-dt') or (option == '-deltat'): # timestep
        dt = float(sys.argv[1])
        del sys.argv[1]
    elif (option == '-n') or (option == '-nsteps'):  # number of timesteps
        nsteps = float(sys.argv[1])
        del sys.argv[1]
    else:
        print sys.argv[0], ': invalid option', option
        sys.exit(1)


#################################
#                               #
# Integrate Equations of Motion #
#                               #
#################################

# create and initialize data lists

# principal integration loop
while time <= dt*nsteps:
    # increment time
    
    # calculate acceleration
    
    # Euler algorithm
    
    # test if pendulum has passed through theta == 0
    # if yes, estimate period
    if (thisTheta*lastTheta < 0):
       # if first crossing, just record time
       # else double time b/w this crossing and last
       # increment reversal counter

##################################
#                                #
# Estimate Period of Oscillation #
#                                #
##################################

# estimate average period including error bar
if nperiod > 1:
    # take average of periods calculated above
    
    # calculate std dev from above average

################################
#                              #
# Estimate Energy Conservation #
#                              #
################################

# find initial and final energies, and relative error b/w them


######################
#                    #
# Print Data to File #
#                    #
######################


#############
#           #
# Plot Data #
#           #
#############

# plot angle and angular velocity as functions of time
pylab.plot(times, theta, label='Theta(t)')
pylab.plot(times, omega, label='Omega(t)')
pylab.legend()
pylab.xlabel("Time, t")
#pylab.show()
pylab.savefig(outfilename + '.pdf')

pylab.clf() # clear figure for next plot

# creat phase portrait
pylab.plot(theta, omega)
pylab.title('Phase Portrait')
pylab.xlabel('Theta')
pylab.ylabel('Omega')
#pylab.show()
pylab.savefig(outfilename + '_phase.pdf')