source: branches/numpy/anuga/shallow_water/eqf_v2.py @ 7198

Last change on this file since 7198 was 6304, checked in by rwilson, 16 years ago

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