# source:trunk/anuga_work/development/sudi/sw_1d/periodic_waves/johns/parameter.py@7933

Last change on this file since 7933 was 7933, checked in by mungkasi, 13 years ago

Renaming a module "bisect.py" to "bisect_function.py" since it crashes with "bisect" in pylab.

File size: 2.2 KB
Line
1from scipy.special import jn
2from scipy import pi, sqrt, linspace, pi, sin, cos
3from config import g
4from Numeric import Float
5from numpy import zeros
6from rootsearch import *
7from bisect_function import *
8
9def j0(x):
10    return jn(0.0, x)
11
12def j1(x):
13    return jn(1.0, x)
14
15def j2(x):
16    return jn(2.0, x)
17
18def j3(x):
19    return jn(3.0, x)
20
21def jm1(x):
22    return jn(-1.0, x)
23
24def jm2(x):
25    return jn(-2.0, x)
26
27
28"""
29##==========================================================================##
30#DIMENSIONAL PARAMETERS
31L = 5e4                 # Length of channel (m)
32h_0 = 5e2               # Height at origin when the water is still
33N = 550                 # Number of computational cells????????????????
34cell_len = 1.1*L/N      # Origin = 0.0
35Tp = 15.0*60.0                           # Period of oscillation
36a = 1.0                                 # Amplitude at origin
37X = linspace(-0.5*cell_len, 1.1*L+0.5*cell_len, N+2)             # Discretized spatial domain
38##=========================================================================##
39"""
40#DIMENSIONAL PARAMETERS
41L = 5e4                 # Length of channel (m)
42h_0 = 5e2               # Height at origin when the water is still
43numb = 550                 # initial Number of computational cells
44cell_len = 100.0#100.0 is actually 1.1*L/numb
45Tp = 15.0*60.0                           # Period of oscillation
46a = 1.0                                  # Amplitude at origin
47
48points = zeros(numb+2,Float)
49for i in range(numb+2):
50        points[i] = i*cell_len - 0.5*cell_len
51N = len(points)-1       # Number of computational cells
52
53
54#DIMENSIONLESS PARAMETERS
55eps = a/h_0
56T = Tp*sqrt(g*h_0)/L
57A = eps/j0(4.0*pi/T)
58
59
60def w_at_O_JOHNS(t):
61    return eps*cos(2.0*pi*t/T)
62
63def u_at_O_JOHNS(t):
64    a = -1.01#-1.0
65    b = 1.01#1.0
66    dx = 0.01
67    w = w_at_O_JOHNS(t)
68    def fun(u):
69        return u + A*j1(4.0*pi/T*(1.0+w)**0.5)*sin(2.0*pi/T*(t+u))/(1+w)**0.5
70    while 1:
71        x1,x2 = rootsearch(fun,a,b,dx)
72        if x1 != None:
73            a = x2
74            root = bisect(fun,x1,x2,1)
75        else:
76            break
77    return root
78
79
80def prescribe_at_O_JOHNS(t):
81    w = w_at_O_JOHNS(t)
82    u = u_at_O_JOHNS(t)
83    return w, u
Note: See TracBrowser for help on using the repository browser.