| 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 okada import okadatest |
|---|
| 115 | |
|---|
| 116 | #ensure vectors x and y have the same length |
|---|
| 117 | N = len(x) |
|---|
| 118 | assert N == len(y) |
|---|
| 119 | |
|---|
| 120 | depth = self.depth |
|---|
| 121 | dip = self.dip |
|---|
| 122 | length = self.length |
|---|
| 123 | width = self.width |
|---|
| 124 | x0 = self.x0 |
|---|
| 125 | y0 = self.y0 |
|---|
| 126 | strike = self.strike |
|---|
| 127 | dip = self.dip |
|---|
| 128 | rake = self.rake |
|---|
| 129 | slip = self.slip |
|---|
| 130 | |
|---|
| 131 | #double Gaussian calculation assumes water displacement is oriented |
|---|
| 132 | #E-W, so, for displacement at some angle alpha clockwise from the E-W |
|---|
| 133 | #direction, rotate (x,y) coordinates anti-clockwise by alpha |
|---|
| 134 | |
|---|
| 135 | cosa = cos(radians(strike)) |
|---|
| 136 | sina = sin(radians(strike)) |
|---|
| 137 | |
|---|
| 138 | xr = ( (x-x0) * sina + (y-y0) * cosa) |
|---|
| 139 | yr = (-(x-x0) * cosa + (y-y0) * sina) |
|---|
| 140 | |
|---|
| 141 | z = okada(xr,yr,depth,length,width,dip,rake,slip) |
|---|
| 142 | |
|---|
| 143 | |
|---|
| 144 | return z |
|---|
| 145 | |
|---|
| 146 | |
|---|
