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

from scipy.special import jn
from scipy import pi, sqrt, linspace, pi, sin, cos
from config import g
from Numeric import Float
from numpy import zeros
from rootsearch import *
from bisect import *

def j0(x):
    return jn(0.0, x)

def j1(x):
    return jn(1.0, x)

def j2(x):
    return jn(2.0, x)

def j3(x):
    return jn(3.0, x)

def jm1(x):
    return jn(-1.0, x)

def jm2(x):
    return jn(-2.0, x)


"""
##==========================================================================##
#DIMENSIONAL PARAMETERS
L = 5e4                 # Length of channel (m)
h_0 = 5e2               # Height at origin when the water is still
N = 550                 # Number of computational cells????????????????
cell_len = 1.1*L/N      # Origin = 0.0
Tp = 15.0*60.0                           # Period of oscillation
a = 1.0                                 # Amplitude at origin
X = linspace(-0.5*cell_len, 1.1*L+0.5*cell_len, N+2)             # Discretized spatial domain
##=========================================================================##
"""
#DIMENSIONAL PARAMETERS
L = 5e4                 # Length of channel (m)
h_0 = 5e2               # Height at origin when the water is still
numb = 550                 # initial Number of computational cells
cell_len = 100.0#100.0 is actually 1.1*L/numb
Tp = 15.0*60.0                           # Period of oscillation
a = 1.0                                  # Amplitude at origin
               
points = zeros(numb+2,Float)
for i in range(numb+2):
        points[i] = i*cell_len - 0.5*cell_len
N = len(points)-1       # Number of computational cells


#DIMENSIONLESS PARAMETERS
eps = a/h_0
T = Tp*sqrt(g*h_0)/L
A = eps/j0(4.0*pi/T)


def w_at_O_JOHNS(t):
    return eps*cos(2.0*pi*t/T)

def u_at_O_JOHNS(t):
    a = -1.01#-1.0
    b = 1.01#1.0
    dx = 0.01
    w = w_at_O_JOHNS(t)
    def fun(u):
        return u + A*j1(4.0*pi/T*(1.0+w)**0.5)*sin(2.0*pi/T*(t+u))/(1+w)**0.5
    while 1:
        x1,x2 = rootsearch(fun,a,b,dx)
        if x1 != None:
            a = x2
            root = bisect(fun,x1,x2,1)
        else:
            break
    return root


def prescribe_at_O_JOHNS(t):
    w = w_at_O_JOHNS(t)
    u = u_at_O_JOHNS(t)
    return w, u