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 | |
---|

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 |
---|