#Answer key to Math Methods Autumn 2005
#Problem Set 1: Introduction to Python.
#
#Run this script to generate the correct
#output. You can compare the result with
#your own script.

#----------------Problem 1----------------
print "------------Problem 1 output-----------"
powers = [3**i for i in range(20)]
print "First 20 powers:", powers
powers.append(3**20) #Next power is 3**20 because range(...)
                     #goes zero to 19
print "Append another one:",powers
powers.insert(0,1./3.) #Remember to use float for this!
print "Insert 1/3",powers
powers.reverse()
print "Reversed:",powers
print "The index of 81 is:",powers.index(81)
print "11th item:",powers[10] #Recall, list is zero-based
print "Now pop off five items:"
for i in range(5):
    print powers.pop()
print "List after popping:",powers

#----------------Problem 2----------------
#Parameters of this function are globals
def alpha(T):
    if T < T_i:
        return alpha_i
    elif (T >= T_i)&(T<=T_o):
        r = (T-T_o)**2/(T_o-T_i)**2
        return alpha_o + (alpha_i - alpha_o)*r
    else:
        return alpha_o

#Set the values of the globals
T_i = 260.
T_o = 290.
alpha_i = .1
alpha_o = .9

#Make a list of T and compute a list of alpha
TList = [250.+ 5.*i for i in range(15)]
alphaList = [alpha(T) for T in TList]
#
#Print out a table
print "------------Problem 2 output-----------"
print "T     alpha(T)"
for i in range(len(TList)):
    print TList[i],alphaList[i]
#
#The following shows how you can make a plot,
#if you want to:
from ClimateUtilities import *
c = Curve() #Create a curve object to store the data
c.addCurve(TList,'T') #The second argument gives the column a name (optional)
c.addCurve(alphaList,'alpha')
plot(c) #or w = plot(c) if you want to give the plot a handle so you can
        #delete or save it later

#----------------Problem 3----------------
def f(n):
    sum = 0.
    for i in range(1,n+1): # range generates list [1,...n]
        sum = sum + 1./i #Note the decimal point forces floating arithmetic
    return sum

print "------------Problem 3 output-----------"
#Try out some values of the function. You can do this interactively
# in the python shell once you've defined the function, just
#by typing things like f(3). Here, we put in a print statement
#so you can run this as a script and see the results.
print f(1)
print f(10)
print f(1000)
print f(100000)
#
#Hm. The sum seems to diverge as the argument gets large.  How fast?
#Let's compare it to ln(n). 
import math
print "Comparison with ln(n)"
for i in range(1,1000,10): #Do computation for every 10th value
    print i,f(i)-math.log(i)

#Now let's look at the behavior out to much larger numbers.  Here,
#we compute the function on powers of 2, so we can go to big numbers
#without needing to print out too much stuff
print "Comparison taken out to larger n:"
for i in range(23):
    n = 2**i
    print n, f(n)-math.log(n)

#The difference appears to converge to a constant. In fact
#this can be proved.  The constant has a name, which is
#the "Euler-Mascheroni constant"

#Expert trick:  you can define a "one liner" function without
#using a def block, by using the "lambda function" syntax
#This is illustrated in the following:
g = lambda i:f(i)-math.log(i)
#
#You can use g now like any other function; e.g.
#g(10) computes f(10) - math.log(10).  Lambda functions
#can have multiple arguments, and can deal with any kind
#of data type at all. For example,
# h = lambda x,y : x**2 + y**2
#returns the sum of the squares of the two arguments.  You
#can write lambda functions to process and return lists, strings
#or anything. The only restriction is that the operation they
#define must consist of a single statement.

#----------------Problem 4----------------
print "----------------Problem 4 output----------------"
def deriv(x,f,eps):
    return (f(x+eps)-f(x-eps))/(2.*eps)


#I assume you already imported math above, so we don't need
#to do it again

xx = [i*math.pi/20. for i in range(10)]
print 'x   exact  eps=.1   eps=.01'
for x in xx:
    print x,math.cos(x),deriv(x,math.sin,.1),deriv(x,math.sin,.01)


#----------------Problem 5----------------
print "----------------Problem 5 output----------------"
class SpaceShip:
    def __init__(self,x,y,u,v):
        self.x = x
        self.y = y
        self.u = u
        self.v = v
        self.time = 0.
        self.dt = 1. #Could make dt an argument to __init__
    def move(self):
        #Note: In constructions like the following, you can
        #also use the shorthand self.x += dt*u
        self.x = self.x + self.dt*self.u
        self.y = self.y + self.dt*self.v
    def distance(self):
        return math.sqrt(self.x*self.x + self.y*self.y)

#Create two instances
ship1 = SpaceShip(0.,0.,1,1.)
ship2 = SpaceShip(6.,6.,-1.,2.)
#Make them go:
for i in range(100):
    print ship1.distance(),ship2.distance()
    ship1.move()
    ship2.move()

#For Experts:
#If you wanted to set dt to be something other than the
#default, you can do, for example, ship1.dt = .5,ship2.dt = .5
#The problem with this is that you usually want the time increment
#to be the same for all spaceships, and you have to remember to
#reset the dt for each one.  In a case like this, it would be
#better to make dt a "class variable" common to all instances of
#the class.  You do this by defining the class as follows:
class BetterSpaceShip:
    dt = .1
    def __init__(self,x,y,u,v):
        self.x = x
        self.y = y
        self.u = u
        self.v = v
        self.time = 0.
    def move(self):
        #Note: In constructions like the following, you can
        #also use the shorthand self.x += dt*u
        self.x = self.x + BetterSpaceShip.dt*self.u
        self.y = self.y + BetterSpaceShip.dt*self.v
        self.time = self.time + dt
    def distance(self):
        return math.sqrt(self.x*self.x + self.y*self.y)
    #
    #The next two methods illustrate how to allow an instance
    #to change the value of a class variable
    def setTimeInc(self,dt):
        BetterSpaceShip.dt = dt
    def doubleTimeInc(self):
        BetterSpaceShip.dt = 2.*BetterSpaceShip.dt
    #    
    #The following method causes an instance to print
    #out its current position when you type its name
    def __repr__(self):
        return 'I am at (%f,%f)'%(self.x,self.y)
#Note that in improving the class, we've also added a __repr__
#method, so that an instance of the class will type out it's
#current position when you type its name.
#
#Now, if you do
s1 = BetterSpaceShip(0.,0.,1.,1.)
s2 = BetterSpaceShip(0.,0.,2.,-2.)
#and then
BetterSpaceShip.dt = .5
#you'll find that both s1.dt and s2.dt are updated.
#Caution: If you did s1.dt = .5 instead, you'd be actually
#over-riding the class variable definition, and making
#an instance variable.  If you wanted to allow an instance
#to change the value, you would follow the example in
#the methods setTimeInc and doubleTimeInc. For example,
#s1.setTimeInc(.5) reset the time increment for all instances of
#the BetterSpaceShip class, and s2.doubleTimeInc() would double
#dt for all members of the class. Try it.
#

#----------------Problem 6----------------
print "----------------Problem 6 output----------------"
class deriv:
    def __init__(self,f,eps=.01): #This gives eps a default value,
                                  #and makes it an optional argument
        self.func = f
        self.eps = eps
    def __call__(self,x):
        return (self.func(x+self.eps)-self.func(x-self.eps))/(2.*self.eps)

#Create some instances:
sineDeriv = deriv(math.sin,.001)
# Note, if we were happy with the default value, we could have just
# written deriv(math.sin)
cosDeriv = deriv(math.cos,.001)
#
#Now try them out
#
print sineDeriv(math.pi/4.),math.cos(math.pi/4.)
print cosDeriv(math.pi/4.),-math.sin(math.pi/4.)