source: anuga_core/source/anuga/shallow_water/eqf.py @ 7675

Last change on this file since 7675 was 6157, checked in by rwilson, 15 years ago

Change Numeric imports to general form - ready to change to NumPy?.

File size: 4.0 KB
Line 
1#
2# earthquake_tsunami function
3#
4
5"""This function returns a callable object representing an initial water
6   displacement generated by a submarine earthqauke.
7   
8Using input parameters:
9
10Required
11 length  along-stike length of rupture area
12 width   down-dip width of rupture area
13 strike  azimuth (degrees, measured from north) of fault axis
14 dip     angle of fault dip in degrees w.r.t. horizontal
15 depth   depth to base of rupture area
16 
17Optional
18 x0      x origin (0)
19 y0      y origin (0)
20 slip    metres of fault slip (1)
21 rake    angle of slip (w.r.t. horizontal) in fault plane (90 degrees)
22
23The returned object is a callable okada function that represents
24the initial water displacement generated by a submarine earthuake.
25
26"""
27
28
29def earthquake_tsunami(length, width, strike, depth, \
30                       dip, x0=0.0, y0=0.0, slip=1.0, rake=90.,\
31                       domain=None, verbose=False):
32
33    from math import sin, radians
34
35    if domain is not None:
36        xllcorner = domain.geo_reference.get_xllcorner()
37        yllcorner = domain.geo_reference.get_yllcorner()
38        x0 = x0 - xllcorner  # fault origin (relative)
39        y0 = y0 - yllcorner
40   
41    #a few temporary print statements
42    if verbose is True:
43        print '\nThe Earthquake ...'
44        print '\tLength: ', length
45        print '\tDepth: ', depth
46        print '\tStrike: ', strike
47        print '\tWidth: ', width
48        print '\tDip: ', dip
49        print '\tSlip: ', slip
50        print '\tx0: ', x0
51        print '\ty0: ', y0
52
53    # warning test
54    test = width*1000.0*sin(radians(dip)) - depth
55
56    if verbose is True:
57        if test > 0.0:
58            print 'Earthquake source not located below seafloor'
59            print 'Please check depth'
60
61    return Okada_func(length=length, width=width, dip=dip, \
62                      x0=x0, y0=y0, strike=strike, depth=depth, \
63                      slip=slip, rake=rake, test=test)
64
65
66#
67# Okada class
68#
69
70"""This is a callable class representing the initial water displacment
71   generated by an earthquake.
72
73Using input parameters:
74
75Required
76 length  along-stike length of rupture area
77 width   down-dip width of rupture area
78 strike  azimuth (degrees, measured from north) of fault axis
79 dip     angle of fault dip in degrees w.r.t. horizontal
80 depth   depth to base of rupture area
81 
82Optional
83 x0      x origin (0)
84 y0      y origin (0)
85 slip    metres of fault slip (1)
86 rake    angle of slip (w.r.t. horizontal) in fault plane (90 degrees)
87
88"""
89
90class Okada_func:
91
92    def __init__(self, length, width, dip, x0, y0, strike, \
93                 depth, slip, rake, test):
94        self.dip = dip
95        self.length = length
96        self.width = width
97        self.x0 = x0
98        self.y0 = y0
99        self.strike = strike
100        self.depth = depth
101        self.slip = slip
102        self.rake = rake
103        self.test = test
104
105
106    def __call__(self, x, y):
107        """Make Okada_func a callable object.
108
109        If called as a function, this object returns z values representing
110        the initial 3D distribution of water heights at the points (x,y)
111        produced by a submarine mass failure.
112        """
113
114        from math import sin, cos, radians, exp, cosh
115        from okada import okadatest
116
117        #ensure vectors x and y have the same length
118        N = len(x)
119        assert N == len(y)
120
121        depth = self.depth
122        dip = self.dip
123        length = self.length
124        width = self.width
125        x0 = self.x0
126        y0 = self.y0
127        strike = self.strike
128        dip = self.dip
129        rake = self.rake
130        slip = self.slip
131       
132        #double Gaussian calculation assumes water displacement is oriented
133        #E-W, so, for displacement at some angle alpha clockwise from the E-W
134        #direction, rotate (x,y) coordinates anti-clockwise by alpha
135
136        cosa = cos(radians(strike))
137        sina = sin(radians(strike))
138
139        xr = ( (x-x0) * sina + (y-y0) * cosa)
140        yr = (-(x-x0) * cosa + (y-y0) * sina) 
141
142        z = okada(xr,yr,depth,length,width,dip,rake,slip)
143   
144               
145        return z
146
147   
Note: See TracBrowser for help on using the repository browser.