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

Last change on this file since 5611 was 4320, checked in by sexton, 17 years ago

new earthquake source function (converted fortran to python), plus new figures for slide report

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
28def earthquake_tsunami(length, width, strike, depth, \
29                       dip, x0=0.0, y0=0.0, slip=1.0, rake=90.,\
30                       domain=None, verbose=False):
31
32    from math import sin, radians
33
34    if domain is not None:
35        xllcorner = domain.geo_reference.get_xllcorner()
36        yllcorner = domain.geo_reference.get_yllcorner()
37        x0 = x0 - xllcorner  # fault origin (relative)
38        y0 = y0 - yllcorner
39   
40    #a few temporary print statements
41    if verbose is True:
42        print '\nThe Earthquake ...'
43        print '\tLength: ', length
44        print '\tDepth: ', depth
45        print '\tStrike: ', strike
46        print '\tWidth: ', width
47        print '\tDip: ', dip
48        print '\tSlip: ', slip
49        print '\tx0: ', x0
50        print '\ty0: ', y0
51
52    # warning test
53    test = width*1000.0*sin(radians(dip)) - depth
54
55    if verbose is True:
56        if test > 0.0:
57            print 'Earthquake source not located below seafloor'
58            print 'Please check depth'
59
60    return Okada_func(length=length, width=width, dip=dip, \
61                      x0=x0, y0=y0, strike=strike, depth=depth, \
62                      slip=slip, rake=rake, test=test)
63
64
65#
66# Okada class
67#
68
69"""This is a callable class representing the initial water displacment
70   generated by an earthquake.
71
72Using input parameters:
73
74Required
75 length  along-stike length of rupture area
76 width   down-dip width of rupture area
77 strike  azimuth (degrees, measured from north) of fault axis
78 dip     angle of fault dip in degrees w.r.t. horizontal
79 depth   depth to base of rupture area
80 
81Optional
82 x0      x origin (0)
83 y0      y origin (0)
84 slip    metres of fault slip (1)
85 rake    angle of slip (w.r.t. horizontal) in fault plane (90 degrees)
86
87"""
88
89class Okada_func:
90
91    def __init__(self, length, width, dip, x0, y0, strike, \
92                 depth, slip, rake, test):
93        self.dip = dip
94        self.length = length
95        self.width = width
96        self.x0 = x0
97        self.y0 = y0
98        self.strike = strike
99        self.depth = depth
100        self.slip = slip
101        self.rake = rake
102        self.test = test
103
104
105    def __call__(self, x, y):
106        """Make Okada_func a callable object.
107
108        If called as a function, this object returns z values representing
109        the initial 3D distribution of water heights at the points (x,y)
110        produced by a submarine mass failure.
111        """
112
113        from math import sin, cos, radians, exp, cosh
114        from Numeric import zeros, Float
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.