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 | |
---|
8 | Using input parameters: |
---|
9 | |
---|
10 | Required |
---|
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 | |
---|
17 | Optional |
---|
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 | |
---|
23 | The returned object is a callable okada function that represents |
---|
24 | the initial water displacement generated by a submarine earthuake. |
---|
25 | |
---|
26 | """ |
---|
27 | |
---|
28 | def 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 | |
---|
72 | Using input parameters: |
---|
73 | |
---|
74 | Required |
---|
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 | |
---|
81 | Optional |
---|
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 | |
---|
89 | class 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 | |
---|