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#
64# Okada class
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
87class Okada_func:
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
112        #from okada import okadatest
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
133        cosa = cos(radians(strike))
134        sina = sin(radians(strike))
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
140        # f2py -c okada.f -m okada
141        #z1 = okada(xr,yr,depth,length,width,dip,rake,slip)
142
143        z2 = zeros(N, Float)
144        alp = 0.5
145        disl3 = 0.0
146        zero = 0.0
147        disl1 = slip*cos(radians(rake))
148        disl2 = slip*sin(radians(rake))
149        cd = cos(radians(dip))
150        sd = sin(radians(dip))
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.