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