source: production/sydney_2006/find_max.py @ 3276

"""
Determining maximum of SMF slump surface elevation function.
Jane Sexton, GA 2006
"""

from math import exp, cosh
import Numeric


# Grilli and Watts 2005 example
dx = 2.0
x0 = 10.0
y0 = 0.0
w = 2.0
kappa = 3
g = 9.8
d = 0.259
T = 0.052
theta = 15.0 
lam0 = 5.0 
    
def func(x,y,kappad):
    numerator = exp(-((x-x0)/lam0)**2.0) - kappad*exp(-((x-dx-x0)/lam0)**2.0)
    denominator = cosh(kappa*(y-y0)/(w+lam0))**2.0
    return -numerator/denominator

def dfuncdx(x,kappad):
    dfdx = (x-x0)*exp(-((x-x0)/lam0)**2.0)-kappad*(x-dx-x0)*exp(-((x-dx-x0)/lam0)**2.0)
    return dfdx

step = 0.001
x = x0
deriv = 10
tol = 0.001
c = 0
direction = 'negative'
xstar = x0
print 'start of loop', deriv
while c < 100000 and deriv > 0:
    deriv = dfuncdx(x,0.83)
    if deriv < 0:
        xstar = x
        print 'hello', deriv, xstar/lam0, c
    if direction == 'positive': x += step
    if direction == 'negative': x -= step
    c += 1
    

print 'location of maximum of surface elevation function: xstar = %f' % (xstar/lam0)
const = 1.0 #a3D = ? #sydney = 86? and kappad=1

#x = arange(-20,25,0.001)
#y = arange(-10,10,0.1)
#X,Y = meshgrid(x,y)

#test = 0.0
#for xi in x:
#    testi = func(xi,0.0,0.83)
#    if direction == 'positive':
#        if testi > test:
#            test = testi
#            xstar = xi
#    else:
#        if testi < test:
#            test = testi
#            xstar = xi
    
#print 'check: xstar = %f and eta(xstar) = %f' %(xstar/lam0, test)