from ClimateUtilities import *
import math

#Set up a function that computes the slope
#as a function of the dependent and independent
#variables.
#
#Here, the proportionality constant "a" is a parameter.
#In this example, we do things the simplest
#way, and treat a as a global, rather than including it
#in the argument list. Any number of parameters
#can be passed this way, and in fact things (like physical constants)
#that you don't need to change are most conveniently passed to the
#function as globals. 

#Here Y is the dependent variable and t is the independent variable
#The independent variable comes first in the argument list
def dYdt(t,Y):
    return -a*Y

#Now create an instance of the integrator for the function dYdt

tStart = 0. #Starting value of independent variable
YStart = 1. #Starting value of dependent variable
dt = .1 #The increment to use. Smaller is more accurate
m = integrator(dYdt,tStart,YStart,dt)

#Set the parameter value. Setting it outside the function
#definition automatically makes it available inside the function,
#so long as the same symbol isn't defined within the function
a = 2.


#The following loop carries out the integration to t = 10.,
#and just prints out the results together with the analytic
#solution for comparison. I use a "for" loop here for simplicity,
#but a better way would be to use a "while" loop to test when
#to stop.

nsteps = int(10./dt)
for i in range(nsteps):
    t,Y = m.next() #The next() method advances the values by one time step
    print t,Y,math.exp(-a*t)


#Now let's do it all over again, but this time saving the results
#in a list and plotting them out. Also, for the sake of variety,
#we'll use a while loop instead of a for loop, just to show how
#it's done.

m = integrator(dYdt,tStart,YStart,dt)
a = 2. #We don't really need to reset a, but what the heck. 

#Lists to store the results
tList = [tStart]
YList = [YStart]

#Carry out the integration
t = 0.
while t< 10.:
    t,Y = m.next() #The next() method advances the values by one time step
    tList.append(t) #Save the new value of t
    YList.append(Y) #Save the new value of Y

#Make a list of the analytic solution for comparison
Yexact = [math.exp(-a*tt) for tt in tList]

#Now put the lists in a Curve() object for output or plotting
c = Curve()
c.addCurve(tList,'t')
c.addCurve(YList,'Y')
c.addCurve(Yexact,'ExactSolution')
c.YLabel = 'Y(t)'
c.XLabel = 't'
c.YlogAxis = True #Use a logarithmic axis so you can see the small values better
plot(c)

#Note than unless you make dt a lot larger, you probably won't be
#able to see the difference between the exact and approximate
#solution on the plot. Experiment with some different values of
#dt.

#**ToDo:
#        Example of using globals to pass parameters
#        Example passing multiple parameters as an array
#        Example passing multiple parameters as an object