#!-/usr/bin/python
# Runge-Kutta (4th order) method of Solving the following System
# Harmonic Oscillator
# y''(t) + y(t) = 0
# y(0)=1, 	y'(0)=0

import matplotlib.pyplot as plt
import numpy as np

h = .01		# step size
tmax = 100*np.pi	# max time
#t = np.arange(0, 100*pi, h)		# array of time values


global y, error, t, ynew, yprev, values
y = np.array([0,0])
values = np.array([[],[]])
val_t = np.array([])
val_calc = np.array([])
error = 0	# initialize error
t = 0
ynew = 0
yprev = 1
k1 = 0
k2 = 0
k3 = 0
k4 = 0

def runge_values(step, x,y):
	return 0


def slope(x,y):
	return -np.sin(x)

while t < tmax :
	k1 = h*slope(t,yprev)
	k2 = h*slope(t+h/2, yprev+k1/2)	#looks like there will be no difference for R-K
	k3 = h*slope(t+h/2, yprev+k2/2)
	k4 = h*slope(t+h, yprev+k3)
	
	ynew = yprev + k1/6 + k2/3 + k3/3 + k4/6
	yprev = ynew
	
	val_t = np.append(val_t, t)
	val_calc = np.append(val_calc, ynew)
	error = error + np.absolute(np.cos(t) - ynew)
	
	t+= h

error_avg = error/len(val_t)
error_final = np.absolute(np.cos(val_t[len(val_t)-1]) - val_calc[len(val_calc)-1])

print "Step Size: %f" % h
print "average error: %f" % error_avg
print "final error: %f" % error_final

print "final slope: %f" % slope(t-h, 0)
plt.plot(val_t, val_calc)
plt.show()