# source:anuga_core/source/anuga/shallow_water/eqf_v2.py@4371

Last change on this file since 4371 was 4371, checked in by nick, 17 years ago

small correction

File size: 25.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 state
53    test = width*1000.0*sin(radians(dip)) - depth
54
55    if test > 0.0:
56        msg = 'Earthquake source not located below seafloor - check depth'
57        raise Exception, msg
58
59    return Okada_func(length=length, width=width, dip=dip, \
60                      x0=x0, y0=y0, strike=strike, depth=depth, \
61                      slip=slip, rake=rake)
62
63#
65#
66
67"""This is a callable class representing the initial water displacment
68   generated by an earthquake.
69
70Using input parameters:
71
72Required
73 length  along-stike length of rupture area
74 width   down-dip width of rupture area
75 strike  azimuth (degrees, measured from north) of fault axis
76 dip     angle of fault dip in degrees w.r.t. horizontal
77 depth   depth to base of rupture area
78
79Optional
80 x0      x origin (0)
81 y0      y origin (0)
82 slip    metres of fault slip (1)
83 rake    angle of slip (w.r.t. horizontal) in fault plane (90 degrees)
84
85"""
86
88
89    def __init__(self, length, width, dip, x0, y0, strike, \
90                 depth, slip, rake):
91        self.dip = dip
92        self.length = length
93        self.width = width
94        self.x0 = x0
95        self.y0 = y0
96        self.strike = strike
97        self.depth = depth
98        self.slip = slip
99        self.rake = rake
100
101
102    def __call__(self, x, y):
103        """Make Okada_func a callable object.
104
105        If called as a function, this object returns z values representing
106        the initial 3D distribution of water heights at the points (x,y)
107        produced by a submarine mass failure.
108        """
109
110        from math import sin, cos, radians, exp, cosh
111        from Numeric import zeros, Float
113
114        #ensure vectors x and y have the same length
115        N = len(x)
116        assert N == len(y)
117
118        depth = self.depth
119        dip = self.dip
120        length = self.length
121        width = self.width
122        x0 = self.x0
123        y0 = self.y0
124        strike = self.strike
125        dip = self.dip
126        rake = self.rake
127        slip = self.slip
128
129        #double Gaussian calculation assumes water displacement is oriented
130        #E-W, so, for displacement at some angle alpha clockwise from the E-W
131        #direction, rotate (x,y) coordinates anti-clockwise by alpha
132
135
136        xr = ( (x-x0) * sina + (y-y0) * cosa)
137        yr = (-(x-x0) * cosa + (y-y0) * sina)
138
139        # works on nautilus when have already done
142
143        z2 = zeros(N, Float)
144        alp = 0.5
145        disl3 = 0.0
146        zero = 0.0
151
152        for i in range(N-1):
153            self.SRECTF(alp,xr[i]*.001,yr[i]*.001,depth*.001,zero,length,\
154                        zero,width,sd,cd,disl1,disl2,disl3)
155
156            z2[i] = self.U3
157
158        return z2
159
160    def SRECTF(self,ALP,X,Y,DEP,AL1,AL2,AW1,AW2,SD,CD,DISL1,DISL2,DISL3):
161        """ SURFACE DISPLACEMENT,STRAIN,TILT DUE TO RECTANGULAR FAULT
162            IN A HALF-SPACE. CODED BY  Y.OKADA ... JAN 1985
163
164
165         INPUT
166           ALP     : MEDIUM CONSTANT  MYU/(LAMDA+MYU)=1./((VP/VS)**2-1)
167           X,Y     : COORDINATE OF STATION
168           DEP     : SOURCE DEPTH
169           AL1,AL2 : FAULT LENGTH RANGE
170           AW1,AW2 : FAULT WIDTH RANGE
171           SD,CD   : SIN,COS OF DIP-ANGLE
172                  (CD=0.D0, SD=+/-1.D0 SHOULD BE GIVEN FOR VERTICAL FAULT)
173           DISL1,DISL2,DISL3 : STRIKE-, DIP- AND TENSILE-DISLOCATION
174
175         OUTPUT
176           U1, U2, U3      : DISPLACEMENT ( UNIT= UNIT OF DISL     )
177           U11,U12,U21,U22 : STRAIN       ( UNIT= UNIT OF DISL /
178           U31,U32         : TILT                 UNIT OF X,Y,,,AW )
179
180         SUBROUTINE USED...SRECTG
181        """
182
183        from Numeric import zeros, Float
184        U = zeros(9, Float)
185        DU = zeros(9, Float)
186
187        F0 = 0.0
188        F1 = 1.0
189
190        P = Y*CD + DEP*SD
191        Q = Y*SD - DEP*CD
192
193        KVEC = [1,2]
194        JVEC = [1,2]
195        for K in KVEC:
196            if K == 1: ET=P-AW1
197            if K == 2: ET=P-AW2
198            for J in JVEC:
199                if J == 1: XI=X-AL1
200                if J == 2: XI=X-AL2
201                JK=J+K
202                if JK <> 3:
203                    SIGN= F1
204                else:
205                    SIGN=-F1
206
207                self.SRECTG(ALP,XI,ET,Q,SD,CD,DISL1,DISL2,DISL3)
208
209                DU[0] = self.DU1
210                DU[1] = self.DU2
211                DU[2] = self.DU3
212                DU[3] = self.DU11
213                DU[4] = self.DU12
214                DU[5] = self.DU21
215                DU[6] = self.DU22
216                DU[7] = self.DU31
217                DU[8] = self.DU32
218
219                for i in range(len(U)):
220                    U[i]=U[i]+SIGN*DU[i]
221
222        U1 =U[0]
223        U2 =U[1]
224        U3 =U[2]
225        U11=U[3]
226        U12=U[4]
227        U21=U[5]
228        U22=U[6]
229        U31=U[7]
230        U32=U[8]
231
232        self.U3 = U3
233
234    def SRECTG(self,ALP,XI,ET,Q,SD,CD,DISL1,DISL2,DISL3):
235        """
236            C
237            C*********************************************************************
238            C*****                                                           *****
239            C*****  INDEFINITE INTEGRAL OF SURFACE DISPLACEMENT,STRAIN,TILT  *****
240            C*****  DUE TO RECTANGULAR FAULT IN A HALF-SPACE                 *****
241            C*****                          CODED BY  Y.OKADA ... JAN 1985   *****
242            C*****                                                           *****
243            C*********************************************************************
244            C
245            C***** INPUT
246            C*****   ALP     : MEDIUM CONSTANT  MYU/(LAMDA+MYU)=1./((VP/VS)**2-1)
247            C*****   XI,ET,Q : FAULT COORDINATE
248            C*****   SD,CD   : SIN,COS OF DIP-ANGLE
249            C*****          (CD=0.D0, SD=+/-1.D0 SHOULD BE GIVEN FOR VERTICAL FAULT)
250            C*****   DISL1,DISL2,DISL3 : STRIKE-, DIP- AND TENSILE-DISLOCATION
251            C
252            C***** OUTPUT
253            C*****   U1, U2, U3      : DISPLACEMENT ( UNIT= UNIT OF DISL    )
254            C*****   U11,U12,U21,U22 : STRAIN       ( UNIT= UNIT OF DISL /
255            C*****   U31,U32         : TILT                 UNIT OF XI,ET,Q )
256            C
257        """
258
259        from math import sqrt, atan, log, radians
260        F0 = 0.0
261        F1 = 1.0
262        F2 = 2.0
263        PI2=6.283185307179586
264
265        XI2=XI*XI
266        ET2=ET*ET
267        Q2=Q*Q
268        R2=XI2+ET2+Q2
269        R =sqrt(R2)
270        R3=R*R2
271        D =ET*SD-Q*CD
272        Y =ET*CD+Q*SD
273        RET=R+ET
274        if RET < F0: RET=F0
275        RD =R+D
276        RRD=F1/(R*RD)
277
278        if Q <> F0:
279            TT = atan( radians( XI*ET/(Q*R) ))
280        else:
281            TT = F0
282
283        if RET <> F0:
284            RE = F1/RET
285            DLE= log(RET)
286        else:
287            RE = F0
288            DLE=-log(R-ET)
289
290        RRX=F1/(R*(R+XI))
291        RRE=RE/R
292        AXI=(F2*R+XI)*RRX*RRX/R
293        AET=(F2*R+ET)*RRE*RRE/R
294
295        if CD == 0:
296            #C==============================
297            #C=====   INCLINED FAULT   =====
298            #C==============================
299            RD2=RD*RD
300            A1=-ALP/F2*XI*Q/RD2
301            A3= ALP/F2*( ET/RD + Y*Q/RD2 - DLE )
302            A4=-ALP*Q/RD
303            A5=-ALP*XI*SD/RD
304            B1= ALP/F2*  Q  /RD2*(F2*XI2*RRD - F1)
305            B2= ALP/F2*XI*SD/RD2*(F2*Q2 *RRD - F1)
306            C1= ALP*XI*Q*RRD/RD
307            C3= ALP*SD/RD*(XI2*RRD - F1)
308        else:
309            #C==============================
310            #C=====   VERTICAL FAULT   =====
311            #C==============================
312            TD=SD/CD
313            X =sqrt(XI2+Q2)
314            if XI == F0:
315                A5=F0
316            else:
317                A5= ALP*F2/CD*atan( radians((ET*(X+Q*CD)+X*(R+X)*SD) / (XI*(R+X)*CD) ))
318            A4= ALP/CD*(  log(RD) - SD*DLE )
319            A3= ALP*(Y/RD/CD - DLE) + TD*A4
320            A1=-ALP/CD*XI/RD        - TD*A5
321            C1= ALP/CD*XI*(RRD - SD*RRE)
322            C3= ALP/CD*(Q*RRE - Y*RRD)
323            B1= ALP/CD*(XI2*RRD - F1)/RD - TD*C3
324            B2= ALP/CD*XI*Y*RRD/RD       - TD*C1
325
326        A2=-ALP*DLE - A3
327        B3=-ALP*XI*RRE - B2
328        B4=-ALP*( CD/R + Q*SD*RRE ) - B1
329        C2= ALP*(-SD/R + Q*CD*RRE ) - C3
330
331        U1 =F0
332        U2 =F0
333        U3 =F0
334        U11=F0
335        U12=F0
336        U21=F0
337        U22=F0
338        U31=F0
339        U32=F0
340
341        if DISL1 <> F0:
342            #C======================================
343            #C=====  STRIKE-SLIP CONTRIBUTION  =====
344            #C======================================
345            UN=DISL1/PI2
346            REQ=RRE*Q
347            U1 =U1 - UN*( REQ*XI +   TT    + A1*SD )
348            U2 =U2 - UN*( REQ*+ Q*CD*RE + A2*SD )
349            U3 =U3 - UN*( REQ*+ Q*SD*RE + A4*SD )
350            U11=U11+ UN*( XI2*Q*AET - B1*SD )
351            U12=U12+ UN*( XI2*XI*( D/(ET2+Q2)/R3 - AET*SD ) - B2*SD )
352            U21=U21+ UN*( XI*Q/R3*CD + (XI*Q2*AET - B2)*SD )
353            U22=U22+ UN*( Y *Q/R3*CD + (Q*SD*(Q2*AET-F2*RRE) -(XI2+ET2)/R3*CD - B4)*SD )
354            U31=U31+ UN*(-XI*Q2*AET*CD + (XI*Q/R3 - C1)*SD )
355            U32=U32+ UN*( D*Q/R3*CD + (XI2*Q*AET*CD - SD/R + Y*Q/R3 - C2)*SD )
356
357        if DISL2 <> F0:
358            #C===================================
359            #C=====  DIP-SLIP CONTRIBUTION  =====
360            #C===================================
361            UN=DISL2/PI2
362            SDCD=SD*CD
363            U1 =U1 - UN*( Q/R             - A3*SDCD )
364            U2 =U2 - UN*( Y*Q*RRX + CD*TT - A1*SDCD )
365            U3 =U3 - UN*( D*Q*RRX + SD*TT - A5*SDCD )
366            U11=U11+ UN*( XI*Q/R3            + B3*SDCD )
367            U12=U12+ UN*( Y *Q/R3 - SD/R     + B1*SDCD )
368            U21=U21+ UN*( Y *Q/R3 + Q*CD*RRE + B1*SDCD )
369            U22=U22+ UN*( Y*Y*Q*AXI - (F2*Y*RRX + XI*CD*RRE)*SD + B2*SDCD )
370            U31=U31+ UN*( D *Q/R3 + Q*SD*RRE + C3*SDCD )
371            U32=U32+ UN*( Y*D*Q*AXI - (F2*D*RRX + XI*SD*RRE)*SD + C1*SDCD )
372
373        if DISL3 <> F0:
374            #C========================================
375            #C=====  TENSILE-FAULT CONTRIBUTION  =====
376            #C========================================
377            UN=DISL3/PI2
378            SDSD=SD*SD
379            U1 =U1 + UN*( Q2*RRE                       - A3*SDSD )
380            U2 =U2 + UN*(-D*Q*RRX - SD*(XI*Q*RRE - TT) - A1*SDSD )
381            U3 =U3 + UN*( Y*Q*RRX + CD*(XI*Q*RRE - TT) - A5*SDSD )
382            U11=U11- UN*( XI*Q2*AET             + B3*SDSD )
383            U12=U12- UN*(-D*Q/R3 - XI2*Q*AET*SD + B1*SDSD )
384            U21=U21- UN*( Q2*(CD/R3 + Q*AET*SD) + B1*SDSD )
385            U22=U22- UN*((Y*CD-D*SD)*Q2*AXI - F2*Q*SD*CD*RRX - (XI*Q2*AET - B2)*SDSD )
386            U31=U31- UN*( Q2*(SD/R3 - Q*AET*CD) + C3*SDSD )
387            U32=U32- UN*((Y*SD+D*CD)*Q2*AXI + XI*Q2*AET*SD*CD - (F2*Q*RRX - C1)*SDSD )
388
389        self.DU1 = U1
390        self.DU2 = U2
391        self.DU3 = U3
392
393        self.DU11 = U11
394        self.DU12 = U12
395
396        self.DU21 = U21
397        self.DU22 = U22
398
399        self.DU31 = U31
400        self.DU32 = U32
401
402    def spoint(self,alp,x,y,dep,sd,cd,pot1,pot2,pot3):
403        """ Calculate surface displacement, strain, tilt due to buried point
404        source in a semiinfinite medium. Y. Okada Jan 1985
405
406        Input:
407
408        ALP   : MEDIUM CONSTANT  MYU/(LAMDA+MYU)=1./((VP/VS)**2-1)
409        X,Y   : COORDINATE OF STATION
410        DEP     : SOURCE DEPTH
411        SD,CD : SIN,COS OF DIP-ANGLE
412               (CD=0.D0, SD=+/-1.D0 SHOULD BE GIVEN FOR VERTICAL FAULT)
413        POT1,POT2,POT3 : STRIKE-, DIP- AND TENSILE-POTENCY
414            POTENCY=(  MOMENT OF DOUBLE-COUPLE  )/MYU    FOR POT1,2
415            POTENCY=(INTENSITY OF ISOTROPIC PART)/LAMDA  FOR POT3
416
417        Output:
418
419        U1, U2, U3      : DISPLACEMENT ( UNIT=(UNIT OF POTENCY) /
420                        :                     (UNIT OF X,Y,D)**2  )
421        U11,U12,U21,U22 : STRAIN       ( UNIT= UNIT OF POTENCY) /
422        U31,U32         : TILT                (UNIT OF X,Y,D)**3  )
423        """
424
425        from math import sqrt
426
427        F0 = 0.0
428        F1 = 1.0
429        F2 = 2.0
430        F3 = 3.0
431        F4 = 4.0
432        F5 = 5.0
433        F8 = 8.0
434        F9 = 9.0
435        PI2 = 6.283185307179586
436
437        D =DEP
438        P =Y*CD + D*SD
439        Q =Y*SD - D*CD
440        S =P*SD + Q*CD
441        X2=X*X
442        Y2=Y*Y
443        XY=X*Y
444        D2=D*D
445        R2=X2 + Y2 + D2
446        R =sqrt(R2)
447        R3=R *R2
448        R5=R3*R2
449        QR=F3*Q/R5
450        XR =F5*X2/R2
451        YR =F5*Y2/R2
452        XYR=F5*XY/R2
453        DR =F5*D /R2
454        RD =R + D
455        R12=F1/(R*RD*RD)
456        R32=R12*(F2*R + D)/ R2
457        R33=R12*(F3*R + D)/(R2*RD)
458        R53=R12*(F8*R2 + F9*R*D + F3*D2)/(R2*R2*RD)
459        R54=R12*(F5*R2 + F4*R*D +    D2)/R3*R12
460
461        A1= ALP*Y*(R12-X2*R33)
462        A2= ALP*X*(R12-Y2*R33)
463        A3= ALP*X/R3 - A2
464        A4=-ALP*XY*R32
465        A5= ALP*( F1/(R*RD) - X2*R32 )
466        B1= ALP*(-F3*XY*R33      + F3*X2*XY*R54)
467        B2= ALP*( F1/R3 - F3*R12 + F3*X2*Y2*R54)
468        B3= ALP*( F1/R3 - F3*X2/R5) - B2
469        B4=-ALP*F3*XY/R5 - B1
470        C1=-ALP*Y*(R32 - X2*R53)
471        C2=-ALP*X*(R32 - Y2*R53)
472        C3=-ALP*F3*X*D/R5 - C2
473
474        U1 =F0
475        U2 =F0
476        U3 =F0
477        U11=F0
478        U12=F0
479        U21=F0
480        U22=F0
481        U31=F0
482        U32=F0
483
484        #======================================
485        #=====  STRIKE-SLIP CONTRIBUTION  =====
486        #======================================
487
488        if POT1 <> F0:
489            UN=POT1/PI2
490            QRX=QR*X
491            FX=F3*X/R5*SD
492            U1 =U1 - UN*( QRX*X + A1*SD )
493            U2 =U2 - UN*( QRX*Y + A2*SD )
494            U3 =U3 - UN*( QRX*D + A4*SD )
495
496            U11=U11- UN*( QRX* (F2-XR)        + B1*SD )
497            U12=U12- UN*(-QRX*XYR      + FX*X + B2*SD )
498            U21=U21- UN*( QR*Y*(F1-XR)        + B2*SD )
499            U22=U22- UN*( QRX *(F1-YR) + FX*Y + B4*SD )
500            U31=U31- UN*( QR*D*(F1-XR)        + C1*SD )
501            U32=U32- UN*(-QRX*DR*Y     + FX*D + C2*SD )
502
503        #======================================
504        #=====    DIP-SLIP CONTRIBUTION   =====
505        #======================================
506
507        if POT2 <> F0:
508            UN=POT2/PI2
509            SDCD=SD*CD
510            QRP=QR*P
511            FS=F3*S/R5
512            U1 =U1 - UN*( QRP*X - A3*SDCD )
513            U2 =U2 - UN*( QRP*Y - A1*SDCD )
514            U3 =U3 - UN*( QRP*D - A5*SDCD )
515            U11=U11- UN*( QRP*(F1-XR)        - B3*SDCD )
516            U12=U12- UN*(-QRP*XYR     + FS*X - B1*SDCD )
517            U21=U21- UN*(-QRP*XYR            - B1*SDCD )
518            U22=U22- UN*( QRP*(F1-YR) + FS*Y - B2*SDCD )
519            U31=U31- UN*(-QRP*DR*X           - C3*SDCD )
520            U32=U32- UN*(-QRP*DR*Y    + FS*D - C1*SDCD )
521
522        #========================================
523        #=====  TENSILE-FAULT CONTRIBUTION  =====
524        #========================================
525
526        if POT3 <> F0:
527            UN=POT3/PI2
528            SDSD=SD*SD
529            QRQ=QR*Q
530            FQ=F2*QR*SD
531            U1 =U1 + UN*( QRQ*X - A3*SDSD )
532            U2 =U2 + UN*( QRQ*Y - A1*SDSD )
533            U3 =U3 + UN*( QRQ*D - A5*SDSD )
534            U11=U11+ UN*( QRQ*(F1-XR)        - B3*SDSD )
535            U12=U12+ UN*(-QRQ*XYR     + FQ*X - B1*SDSD )
536            U21=U21+ UN*(-QRQ*XYR            - B1*SDSD )
537            U22=U22+ UN*( QRQ*(F1-YR) + FQ*Y - B2*SDSD )
538            U31=U31+ UN*(-QRQ*DR*X           - C3*SDSD )
539            U32=U32+ UN*(-QRQ*DR*Y    + FQ*D - C1*SDSD )
540
Note: See TracBrowser for help on using the repository browser.