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