#This scripts illustrates how to use
#the integrator in ClimateUtilities to
#numerically solve the ordinary differential equation
#
#                dY/dx = f(x,Y)
#
#for the function Y(x)

import math
from ClimateUtilities import *

#Define the function that computes the derivatives.
#the third argument allows us to pass parameters to
#the function.  Here, the parameter is a real number,
#but it can be any python object at all (e.g. a list of
#numbers, or even a function).
def f(x,Y,a):
    return -a*x*Y

#Create an instance of the integrator object
dx = .1 #The increment to use.
        #Smaller values give more accurate results
Ystart = 1. #The initial value of Y
xstart = 0. #The initial value of x
#
m = integrator(f,xstart,Ystart,dx)
#
#Set the parameter value
a = 2.
m.setParams(a) #Sets the parameter  to the value of a
#
#Make lists to save the results for plotting
xList = [xstart]
YList = [Ystart]
#
#Step through until x exceeds the desired target
x = xstart
while x < 5.:
    x,Y = m.next() # next() can take a new value of dx as optional argument
    #Put the values in the results list
    xList.append(x)
    YList.append(Y)
#Compute the exact solution for comparison
YexactList = [Ystart*math.exp(-a*x*x/2.) for x in xList]
#Put results in a Curve object and plot
c = Curve()
c.addCurve(xList,'x')
c.addCurve(YList,'Y','Computed')
c.addCurve(YexactList,'Yexact','Exact')
w = plot(c)
#With dx = .1, the approximate solution is so accurate you
#can't see the difference between the exact and approximate
#curves. Try larger values of dx to see how the accuracy
#degrades at larger dx. The error shows up first at the
#larger values of x
#Alternately, we can compute the difference and plot that:
c1 = Curve()
c1.addCurve(xList,'x')
c1['diff'] = c['Y']-c['Yexact'] #creates new variable named 'diff',
                                #which is the difference between
                                #Y and Yexact
w1 = plot(c1)