source: anuga_core/source/anuga/shallow_water/eqf_v2.py @ 5611

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

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

File size: 25.1 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    test = width*1000.0*sin(radians(dip)) - depth*1000
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#
65# Okada class
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
88class Okada_func:
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
113        #from okada import okadatest
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
134        cosa = cos(radians(strike))
135        sina = sin(radians(strike))
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
141        # f2py -c okada.f -m okada
142        #z1 = okada(xr,yr,depth,length,width,dip,rake,slip)
143
144        z2 = zeros(N, Float)
145        alp = 0.5
146        disl3 = 0.0
147        zero = 0.0
148        disl1 = slip*cos(radians(rake))
149        disl2 = slip*sin(radians(rake))
150        cd = cos(radians(dip))
151        sd = sin(radians(dip))
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
449        R5=R3*R2
450        QR=F3*Q/R5
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.