source: production/sydney_2006/find_max.py @ 3171

"""
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'
while c < 100000 and deriv > 0:
    deriv = dfuncdx(x,0.83)
    if deriv < 0: xstar = x
    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)