source: anuga_work/production/sydney_2006/find_max.py @ 5001

Last change on this file since 5001 was 3190, checked in by sexton, 18 years ago

MOST and ANUGA comparisons and updates

File size: 1.5 KB
Line 
1"""
2Determining maximum of SMF slump surface elevation function.
3Jane Sexton, GA 2006
4"""
5
6from math import exp, cosh
7import Numeric
8
9
10# Grilli and Watts 2005 example
11dx = 2.0
12x0 = 10.0
13y0 = 0.0
14w = 2.0
15kappa = 3
16g = 9.8
17d = 0.259
18T = 0.052
19theta = 15.0 
20lam0 = 5.0 
21   
22def func(x,y,kappad):
23    numerator = exp(-((x-x0)/lam0)**2.0) - kappad*exp(-((x-dx-x0)/lam0)**2.0)
24    denominator = cosh(kappa*(y-y0)/(w+lam0))**2.0
25    return -numerator/denominator
26
27def dfuncdx(x,kappad):
28    dfdx = (x-x0)*exp(-((x-x0)/lam0)**2.0)-kappad*(x-dx-x0)*exp(-((x-dx-x0)/lam0)**2.0)
29    return dfdx
30
31step = 0.001
32x = x0
33deriv = 10
34tol = 0.001
35c = 0
36direction = 'negative'
37xstar = x0
38print 'start of loop', deriv
39while c < 100000 and deriv > 0:
40    deriv = dfuncdx(x,0.83)
41    if deriv < 0:
42        xstar = x
43        print 'hello', deriv, xstar/lam0, c
44    if direction == 'positive': x += step
45    if direction == 'negative': x -= step
46    c += 1
47   
48
49print 'location of maximum of surface elevation function: xstar = %f' % (xstar/lam0)
50const = 1.0 #a3D = ? #sydney = 86? and kappad=1
51
52#x = arange(-20,25,0.001)
53#y = arange(-10,10,0.1)
54#X,Y = meshgrid(x,y)
55
56#test = 0.0
57#for xi in x:
58#    testi = func(xi,0.0,0.83)
59#    if direction == 'positive':
60#        if testi > test:
61#            test = testi
62#            xstar = xi
63#    else:
64#        if testi < test:
65#            test = testi
66#            xstar = xi
67   
68#print 'check: xstar = %f and eta(xstar) = %f' %(xstar/lam0, test)
Note: See TracBrowser for help on using the repository browser.