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

Last change on this file since 5539 was 4436, checked in by nick, 16 years ago

clean up and update source function 'eqf however there are still problems with it

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