# # earthquake_tsunami function # """This function returns a callable object representing an initial water displacement generated by a submarine earthqauke. Using input parameters: Required length along-stike length of rupture area width down-dip width of rupture area strike azimuth (degrees, measured from north) of fault axis dip angle of fault dip in degrees w.r.t. horizontal depth depth to base of rupture area Optional x0 x origin (0) y0 y origin (0) slip metres of fault slip (1) rake angle of slip (w.r.t. horizontal) in fault plane (90 degrees) The returned object is a callable okada function that represents the initial water displacement generated by a submarine earthuake. """ import numpy as num def earthquake_tsunami(length, width, strike, depth, \ dip, x0=0.0, y0=0.0, slip=1.0, rake=90.,\ domain=None, verbose=False): from math import sin, radians if domain is not None: xllcorner = domain.geo_reference.get_xllcorner() yllcorner = domain.geo_reference.get_yllcorner() x0 = x0 - xllcorner # fault origin (relative) y0 = y0 - yllcorner #a few temporary print statements if verbose is True: print '\nThe Earthquake ...' print '\tLength: ', length print '\tDepth: ', depth print '\tStrike: ', strike print '\tWidth: ', width print '\tDip: ', dip print '\tSlip: ', slip print '\tx0: ', x0 print '\ty0: ', y0 # warning state # test = width*1000.0*sin(radians(dip)) - depth test = width*1000.0*sin(radians(dip)) - depth*1000 if test > 0.0: msg = 'Earthquake source not located below seafloor - check depth' raise Exception, msg return Okada_func(length=length, width=width, dip=dip, \ x0=x0, y0=y0, strike=strike, depth=depth, \ slip=slip, rake=rake) # # Okada class # """This is a callable class representing the initial water displacment generated by an earthquake. Using input parameters: Required length along-stike length of rupture area width down-dip width of rupture area strike azimuth (degrees, measured from north) of fault axis dip angle of fault dip in degrees w.r.t. horizontal depth depth to base of rupture area Optional x0 x origin (0) y0 y origin (0) slip metres of fault slip (1) rake angle of slip (w.r.t. horizontal) in fault plane (90 degrees) """ class Okada_func: def __init__(self, length, width, dip, x0, y0, strike, \ depth, slip, rake): self.dip = dip self.length = length self.width = width self.x0 = x0 self.y0 = y0 self.strike = strike self.depth = depth self.slip = slip self.rake = rake def __call__(self, x, y): """Make Okada_func a callable object. If called as a function, this object returns z values representing the initial 3D distribution of water heights at the points (x,y) produced by a submarine mass failure. """ from math import sin, cos, radians, exp, cosh #from okada import okadatest #ensure vectors x and y have the same length N = len(x) assert N == len(y) depth = self.depth dip = self.dip length = self.length width = self.width x0 = self.x0 y0 = self.y0 strike = self.strike dip = self.dip rake = self.rake slip = self.slip #double Gaussian calculation assumes water displacement is oriented #E-W, so, for displacement at some angle alpha clockwise from the E-W #direction, rotate (x,y) coordinates anti-clockwise by alpha cosa = cos(radians(strike)) sina = sin(radians(strike)) xr = ( (x-x0) * sina + (y-y0) * cosa) yr = (-(x-x0) * cosa + (y-y0) * sina) # works on nautilus when have already done # f2py -c okada.f -m okada #z1 = okada(xr,yr,depth,length,width,dip,rake,slip) z2 = num.zeros(N, num.float) alp = 0.5 disl3 = 0.0 zero = 0.0 disl1 = slip*cos(radians(rake)) disl2 = slip*sin(radians(rake)) cd = cos(radians(dip)) sd = sin(radians(dip)) for i in range(N-1): self.SRECTF(alp,xr[i]*.001,yr[i]*.001,depth*.001,zero,length,\ zero,width,sd,cd,disl1,disl2,disl3) z2[i] = self.U3 return z2 def SRECTF(self,ALP,X,Y,DEP,AL1,AL2,AW1,AW2,SD,CD,DISL1,DISL2,DISL3): """ SURFACE DISPLACEMENT,STRAIN,TILT DUE TO RECTANGULAR FAULT IN A HALF-SPACE. CODED BY Y.OKADA ... JAN 1985 INPUT ALP : MEDIUM CONSTANT MYU/(LAMDA+MYU)=1./((VP/VS)**2-1) X,Y : COORDINATE OF STATION DEP : SOURCE DEPTH AL1,AL2 : FAULT LENGTH RANGE AW1,AW2 : FAULT WIDTH RANGE SD,CD : SIN,COS OF DIP-ANGLE (CD=0.D0, SD=+/-1.D0 SHOULD BE GIVEN FOR VERTICAL FAULT) DISL1,DISL2,DISL3 : STRIKE-, DIP- AND TENSILE-DISLOCATION OUTPUT U1, U2, U3 : DISPLACEMENT ( UNIT= UNIT OF DISL ) U11,U12,U21,U22 : STRAIN ( UNIT= UNIT OF DISL / U31,U32 : TILT UNIT OF X,Y,,,AW ) SUBROUTINE USED...SRECTG """ U = num.zeros(9, num.float) DU = num.zeros(9, num.float) F0 = 0.0 F1 = 1.0 P = Y*CD + DEP*SD Q = Y*SD - DEP*CD KVEC = [1,2] JVEC = [1,2] for K in KVEC: if K == 1: ET=P-AW1 if K == 2: ET=P-AW2 for J in JVEC: if J == 1: XI=X-AL1 if J == 2: XI=X-AL2 JK=J+K if JK <> 3: SIGN= F1 else: SIGN=-F1 self.SRECTG(ALP,XI,ET,Q,SD,CD,DISL1,DISL2,DISL3) DU[0] = self.DU1 DU[1] = self.DU2 DU[2] = self.DU3 DU[3] = self.DU11 DU[4] = self.DU12 DU[5] = self.DU21 DU[6] = self.DU22 DU[7] = self.DU31 DU[8] = self.DU32 for i in range(len(U)): U[i]=U[i]+SIGN*DU[i] U1 =U[0] U2 =U[1] U3 =U[2] U11=U[3] U12=U[4] U21=U[5] U22=U[6] U31=U[7] U32=U[8] self.U3 = U3 def SRECTG(self,ALP,XI,ET,Q,SD,CD,DISL1,DISL2,DISL3): """ C C********************************************************************* C***** ***** C***** INDEFINITE INTEGRAL OF SURFACE DISPLACEMENT,STRAIN,TILT ***** C***** DUE TO RECTANGULAR FAULT IN A HALF-SPACE ***** C***** CODED BY Y.OKADA ... JAN 1985 ***** C***** ***** C********************************************************************* C C***** INPUT C***** ALP : MEDIUM CONSTANT MYU/(LAMDA+MYU)=1./((VP/VS)**2-1) C***** XI,ET,Q : FAULT COORDINATE C***** SD,CD : SIN,COS OF DIP-ANGLE C***** (CD=0.D0, SD=+/-1.D0 SHOULD BE GIVEN FOR VERTICAL FAULT) C***** DISL1,DISL2,DISL3 : STRIKE-, DIP- AND TENSILE-DISLOCATION C C***** OUTPUT C***** U1, U2, U3 : DISPLACEMENT ( UNIT= UNIT OF DISL ) C***** U11,U12,U21,U22 : STRAIN ( UNIT= UNIT OF DISL / C***** U31,U32 : TILT UNIT OF XI,ET,Q ) C """ from math import sqrt, atan, log, radians F0 = 0.0 F1 = 1.0 F2 = 2.0 PI2=6.283185307179586 XI2=XI*XI ET2=ET*ET Q2=Q*Q R2=XI2+ET2+Q2 R =sqrt(R2) R3=R*R2 D =ET*SD-Q*CD Y =ET*CD+Q*SD RET=R+ET if RET < F0: RET=F0 RD =R+D RRD=F1/(R*RD) if Q <> F0: TT = atan( radians( XI*ET/(Q*R) )) else: TT = F0 if RET <> F0: RE = F1/RET DLE= log(RET) else: RE = F0 DLE=-log(R-ET) RRX=F1/(R*(R+XI)) RRE=RE/R AXI=(F2*R+XI)*RRX*RRX/R AET=(F2*R+ET)*RRE*RRE/R if CD == 0: #C============================== #C===== INCLINED FAULT ===== #C============================== RD2=RD*RD A1=-ALP/F2*XI*Q/RD2 A3= ALP/F2*( ET/RD + Y*Q/RD2 - DLE ) A4=-ALP*Q/RD A5=-ALP*XI*SD/RD B1= ALP/F2* Q /RD2*(F2*XI2*RRD - F1) B2= ALP/F2*XI*SD/RD2*(F2*Q2 *RRD - F1) C1= ALP*XI*Q*RRD/RD C3= ALP*SD/RD*(XI2*RRD - F1) else: #C============================== #C===== VERTICAL FAULT ===== #C============================== TD=SD/CD X =sqrt(XI2+Q2) if XI == F0: A5=F0 else: A5= ALP*F2/CD*atan( radians((ET*(X+Q*CD)+X*(R+X)*SD) / (XI*(R+X)*CD) )) A4= ALP/CD*( log(RD) - SD*DLE ) A3= ALP*(Y/RD/CD - DLE) + TD*A4 A1=-ALP/CD*XI/RD - TD*A5 C1= ALP/CD*XI*(RRD - SD*RRE) C3= ALP/CD*(Q*RRE - Y*RRD) B1= ALP/CD*(XI2*RRD - F1)/RD - TD*C3 B2= ALP/CD*XI*Y*RRD/RD - TD*C1 A2=-ALP*DLE - A3 B3=-ALP*XI*RRE - B2 B4=-ALP*( CD/R + Q*SD*RRE ) - B1 C2= ALP*(-SD/R + Q*CD*RRE ) - C3 U1 =F0 U2 =F0 U3 =F0 U11=F0 U12=F0 U21=F0 U22=F0 U31=F0 U32=F0 if DISL1 <> F0: #C====================================== #C===== STRIKE-SLIP CONTRIBUTION ===== #C====================================== UN=DISL1/PI2 REQ=RRE*Q U1 =U1 - UN*( REQ*XI + TT + A1*SD ) U2 =U2 - UN*( REQ*Y + Q*CD*RE + A2*SD ) U3 =U3 - UN*( REQ*D + Q*SD*RE + A4*SD ) U11=U11+ UN*( XI2*Q*AET - B1*SD ) U12=U12+ UN*( XI2*XI*( D/(ET2+Q2)/R3 - AET*SD ) - B2*SD ) U21=U21+ UN*( XI*Q/R3*CD + (XI*Q2*AET - B2)*SD ) U22=U22+ UN*( Y *Q/R3*CD + (Q*SD*(Q2*AET-F2*RRE) -(XI2+ET2)/R3*CD - B4)*SD ) U31=U31+ UN*(-XI*Q2*AET*CD + (XI*Q/R3 - C1)*SD ) U32=U32+ UN*( D*Q/R3*CD + (XI2*Q*AET*CD - SD/R + Y*Q/R3 - C2)*SD ) if DISL2 <> F0: #C=================================== #C===== DIP-SLIP CONTRIBUTION ===== #C=================================== UN=DISL2/PI2 SDCD=SD*CD U1 =U1 - UN*( Q/R - A3*SDCD ) U2 =U2 - UN*( Y*Q*RRX + CD*TT - A1*SDCD ) U3 =U3 - UN*( D*Q*RRX + SD*TT - A5*SDCD ) U11=U11+ UN*( XI*Q/R3 + B3*SDCD ) U12=U12+ UN*( Y *Q/R3 - SD/R + B1*SDCD ) U21=U21+ UN*( Y *Q/R3 + Q*CD*RRE + B1*SDCD ) U22=U22+ UN*( Y*Y*Q*AXI - (F2*Y*RRX + XI*CD*RRE)*SD + B2*SDCD ) U31=U31+ UN*( D *Q/R3 + Q*SD*RRE + C3*SDCD ) U32=U32+ UN*( Y*D*Q*AXI - (F2*D*RRX + XI*SD*RRE)*SD + C1*SDCD ) if DISL3 <> F0: #C======================================== #C===== TENSILE-FAULT CONTRIBUTION ===== #C======================================== UN=DISL3/PI2 SDSD=SD*SD U1 =U1 + UN*( Q2*RRE - A3*SDSD ) U2 =U2 + UN*(-D*Q*RRX - SD*(XI*Q*RRE - TT) - A1*SDSD ) U3 =U3 + UN*( Y*Q*RRX + CD*(XI*Q*RRE - TT) - A5*SDSD ) U11=U11- UN*( XI*Q2*AET + B3*SDSD ) U12=U12- UN*(-D*Q/R3 - XI2*Q*AET*SD + B1*SDSD ) U21=U21- UN*( Q2*(CD/R3 + Q*AET*SD) + B1*SDSD ) U22=U22- UN*((Y*CD-D*SD)*Q2*AXI - F2*Q*SD*CD*RRX - (XI*Q2*AET - B2)*SDSD ) U31=U31- UN*( Q2*(SD/R3 - Q*AET*CD) + C3*SDSD ) U32=U32- UN*((Y*SD+D*CD)*Q2*AXI + XI*Q2*AET*SD*CD - (F2*Q*RRX - C1)*SDSD ) self.DU1 = U1 self.DU2 = U2 self.DU3 = U3 self.DU11 = U11 self.DU12 = U12 self.DU21 = U21 self.DU22 = U22 self.DU31 = U31 self.DU32 = U32 def spoint(self,alp,x,y,dep,sd,cd,pot1,pot2,pot3): """ Calculate surface displacement, strain, tilt due to buried point source in a semiinfinite medium. Y. Okada Jan 1985 Input: ALP : MEDIUM CONSTANT MYU/(LAMDA+MYU)=1./((VP/VS)**2-1) X,Y : COORDINATE OF STATION DEP : SOURCE DEPTH SD,CD : SIN,COS OF DIP-ANGLE (CD=0.D0, SD=+/-1.D0 SHOULD BE GIVEN FOR VERTICAL FAULT) POT1,POT2,POT3 : STRIKE-, DIP- AND TENSILE-POTENCY POTENCY=( MOMENT OF DOUBLE-COUPLE )/MYU FOR POT1,2 POTENCY=(INTENSITY OF ISOTROPIC PART)/LAMDA FOR POT3 Output: U1, U2, U3 : DISPLACEMENT ( UNIT=(UNIT OF POTENCY) / : (UNIT OF X,Y,D)**2 ) U11,U12,U21,U22 : STRAIN ( UNIT= UNIT OF POTENCY) / U31,U32 : TILT (UNIT OF X,Y,D)**3 ) """ from math import sqrt F0 = 0.0 F1 = 1.0 F2 = 2.0 F3 = 3.0 F4 = 4.0 F5 = 5.0 F8 = 8.0 F9 = 9.0 PI2 = 6.283185307179586 D =DEP P =Y*CD + D*SD Q =Y*SD - D*CD S =P*SD + Q*CD X2=X*X Y2=Y*Y XY=X*Y D2=D*D R2=X2 + Y2 + D2 R =sqrt(R2) R3=R *R2 R5=R3*R2 QR=F3*Q/R5 XR =F5*X2/R2 YR =F5*Y2/R2 XYR=F5*XY/R2 DR =F5*D /R2 RD =R + D R12=F1/(R*RD*RD) R32=R12*(F2*R + D)/ R2 R33=R12*(F3*R + D)/(R2*RD) R53=R12*(F8*R2 + F9*R*D + F3*D2)/(R2*R2*RD) R54=R12*(F5*R2 + F4*R*D + D2)/R3*R12 A1= ALP*Y*(R12-X2*R33) A2= ALP*X*(R12-Y2*R33) A3= ALP*X/R3 - A2 A4=-ALP*XY*R32 A5= ALP*( F1/(R*RD) - X2*R32 ) B1= ALP*(-F3*XY*R33 + F3*X2*XY*R54) B2= ALP*( F1/R3 - F3*R12 + F3*X2*Y2*R54) B3= ALP*( F1/R3 - F3*X2/R5) - B2 B4=-ALP*F3*XY/R5 - B1 C1=-ALP*Y*(R32 - X2*R53) C2=-ALP*X*(R32 - Y2*R53) C3=-ALP*F3*X*D/R5 - C2 U1 =F0 U2 =F0 U3 =F0 U11=F0 U12=F0 U21=F0 U22=F0 U31=F0 U32=F0 #====================================== #===== STRIKE-SLIP CONTRIBUTION ===== #====================================== if POT1 <> F0: UN=POT1/PI2 QRX=QR*X FX=F3*X/R5*SD U1 =U1 - UN*( QRX*X + A1*SD ) U2 =U2 - UN*( QRX*Y + A2*SD ) U3 =U3 - UN*( QRX*D + A4*SD ) U11=U11- UN*( QRX* (F2-XR) + B1*SD ) U12=U12- UN*(-QRX*XYR + FX*X + B2*SD ) U21=U21- UN*( QR*Y*(F1-XR) + B2*SD ) U22=U22- UN*( QRX *(F1-YR) + FX*Y + B4*SD ) U31=U31- UN*( QR*D*(F1-XR) + C1*SD ) U32=U32- UN*(-QRX*DR*Y + FX*D + C2*SD ) #====================================== #===== DIP-SLIP CONTRIBUTION ===== #====================================== if POT2 <> F0: UN=POT2/PI2 SDCD=SD*CD QRP=QR*P FS=F3*S/R5 U1 =U1 - UN*( QRP*X - A3*SDCD ) U2 =U2 - UN*( QRP*Y - A1*SDCD ) U3 =U3 - UN*( QRP*D - A5*SDCD ) U11=U11- UN*( QRP*(F1-XR) - B3*SDCD ) U12=U12- UN*(-QRP*XYR + FS*X - B1*SDCD ) U21=U21- UN*(-QRP*XYR - B1*SDCD ) U22=U22- UN*( QRP*(F1-YR) + FS*Y - B2*SDCD ) U31=U31- UN*(-QRP*DR*X - C3*SDCD ) U32=U32- UN*(-QRP*DR*Y + FS*D - C1*SDCD ) #======================================== #===== TENSILE-FAULT CONTRIBUTION ===== #======================================== if POT3 <> F0: UN=POT3/PI2 SDSD=SD*SD QRQ=QR*Q FQ=F2*QR*SD U1 =U1 + UN*( QRQ*X - A3*SDSD ) U2 =U2 + UN*( QRQ*Y - A1*SDSD ) U3 =U3 + UN*( QRQ*D - A5*SDSD ) U11=U11+ UN*( QRQ*(F1-XR) - B3*SDSD ) U12=U12+ UN*(-QRQ*XYR + FQ*X - B1*SDSD ) U21=U21+ UN*(-QRQ*XYR - B1*SDSD ) U22=U22+ UN*( QRQ*(F1-YR) + FQ*Y - B2*SDSD ) U31=U31+ UN*(-QRQ*DR*X - C3*SDSD ) U32=U32+ UN*(-QRQ*DR*Y + FQ*D - C1*SDSD )