# # 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 (km) width down-dip width of rupture area (km) 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 (km) ns the really used number of rectangular sources NSMAX the upper limit of ns Optional x0 x origin (0) y0 y origin (0) z0 z origin 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(ns,NSMAX,length, width, strike, depth,\ dip, xi, yi,z0, slip, rake,\ domain=None, verbose=False): from anuga.abstract_2d_finite_volumes.quantity import Quantity from math import sin, radians #zrec0 = Quantity(domain) zrec0 = Quantity(domain) zrec0.set_values(z0) zrec=zrec0.get_vertex_values(xy=True) x0= num.zeros(ns,num.float) y0= num.zeros(ns,num.float) if ns ==1: x0[0]=xi y0[0]=yi else: x0=xi y0=yi if domain is not None: xllcorner = domain.geo_reference.get_xllcorner() yllcorner = domain.geo_reference.get_yllcorner() for i in range(0,ns): x0[i] = x0[i] - xllcorner # fault origin (relative) y0[i] = y0[i] - yllcorner #a few temporary print statements if verbose is True: print '\nThe Earthquake ...' print '\tns: ', ns print '\tNSMAX: ', NSMAX 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(ns=ns,NSMAX=NSMAX,length=length, width=width, dip=dip, \ x0=x0, y0=y0, strike=strike, depth=depth, \ slip=slip, rake=rake, zrec=zrec) # # 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, ns,NSMAX,length, width, dip, x0, y0, strike, \ depth, slip, rake,zrec): 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 self.ns=ns self.zrec=zrec 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,z) produced by a submarine mass failure. """ from string import replace,strip from math import sin, cos, radians, exp, cosh #ensure vectors x and y have the same length N = len(x) assert N == len(y) #nrec=N*N depth = self.depth dip = self.dip length = self.length width = self.width y0 = self.x0 x0 = self.y0 strike = self.strike dip = self.dip rake = self.rake dislocation = self.slip ns=self.ns zrec=self.zrec #initialization disp0=num.zeros(3,num.float) strain0=num.zeros(6,num.float) tilt0 = num.zeros(2,num.float) dislocations=num.zeros(ns,num.float) depths=num.zeros(ns,num.float) strikes= num.zeros(ns,num.float) lengths= num.zeros(ns,num.float) slips= num.zeros(ns,num.float) rakes= num.zeros(ns,num.float) widths= num.zeros(ns,num.float) dips= num.zeros(ns,num.float) strikes= num.zeros(ns,num.float) strikes= num.zeros(ns,num.float) strain = num.zeros((N,6),num.float) disp = num.zeros((N,3),num.float) tilt = num.zeros((N,2),num.float) xs =num.zeros(ns,num.float) ys =num.zeros(ns,num.float) z=[] if ns==1: dislocations[0]=dislocation depths[0]=depth strikes[0]=strike lengths[0]=length rakes[0]=rake widths[0]=width dips[0]=dip try: xs[0]=x0 ys[0]=y0 except: xs[0]=x0[0] ys[0]=y0[0] else: dislocations=dislocation strikes=strike lengths=length rakes=rake widths=width dips=dip xs=x0 ys=y0 depths=depth #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 ALPHA = 0.5 POT3 = 0.0 POT4 = 0.0 DISL3 = 0.0 AL1 = 0.0 AW2 = 0.0 #rad=1.745329252e-2 #zrec = domain.get_quantity('elevation').get_values(interpolation_points=[[x, y]]) #get Zrec... has to be constant!!! #zrec=zreci.get_values(interpolation_points=[[x, y]]) #zrec=0 #Z=-zrec eps = 1.0e-6 # for irec in range(0,N): xrec=y yrec=x for i in range(0,len(zrec[0])): if (num.alltrue(zrec[0][i]==yrec) and num.alltrue(zrec[1][i]==xrec)): Z=zrec[2][i] Z=0.001*Z break else: continue #zrec=zreci.get_values(interpolation_points=[[xrec[irec],yrec[irec]]],location='edges') for ist in range(0,ns): #st = radians(strikes[ist]) st=radians(strikes[ist]) csst = cos(st) ssst = sin(st) cs2st = cos(2.0*st) ss2st = sin(2.0*st) # #di = radians(dips[ist]) di=radians(dips[ist]) csdi =cos(di) ssdi=sin(di) # #ra= radians(rakes[ist]) ra=radians(rakes[ist]) csra=cos(ra) ssra=sin(ra) # transform from Aki's to Okada's system #first attribute x to north and y to east to match OKADA axis #X=0.001*(x[irec]-xs[ist])*ssst+0.001*(y[irec]-ys[ist])*csst #Y=-0.001*(x[irec]-xs[ist])*csst+0.001*(y[irec]-ys[ist])*ssst X=0.001*((xrec[irec]-xs[ist])*csst+(yrec[irec]-ys[ist])*ssst) Y=0.001*((xrec[irec]-xs[ist])*ssst-(yrec[irec]-ys[ist])*csst) DEPTH=depths[ist] DIP=dips[ist] if lengths[ist]==0 and widths[ist]== 0 : # # point source # POT1=dislocations[ist]*csra POT2=dislocations[ist]*ssra IRET=1 self.DC3D0(ALPHA,X,Y,Z,DEPTH,DIP,POT1,POT2,POT3,POT4) UX=self.UX UY=self.UY UZ=self.UZ UXX=self.UXX UYX=self.UYX UZX=self.UZX UXY=self.UXY UYY=self.UYY UZY=self.UZY UXZ=self.UXZ UYZ=self.UYZ UZZ=self.UZZ IRET=self.IRET if IRET==1: print ' There is a problem in Okada subroutine!' break else: # finite source AL2=lengths[ist] AW1=-widths[ist] DISL1=dislocations[ist]*csra DISL2=dislocations[ist]*ssra #print 'lenghts[',ist,']',lenghts[ist] if lengths[ist]==0: #print 'beh alors' AL2=widths[ist]*eps DISL1=DISL1/AL2 DISL2=DISL2/AL2 elif widths[ist]==0.0: AW1=-lengths[ist]*eps DISL1=DISL1/(-AW1) DISL2=DISL2/(-AW1) IRET=1 self.DC3D(ALPHA,X,Y,Z,DEPTH,DIP,AL1,AL2,AW1,AW2,\ DISL1,DISL2,DISL3) UX=self.UX UY=self.UY UZ=self.UZ UXX=self.UXX UYX=self.UYX UZX=self.UZX UXY=self.UXY UYY=self.UYY UZY=self.UZY UXZ=self.UXZ UYZ=self.UYZ UZZ=self.UZZ IRET=self.IRET #if X==-6 and Y==-1: #print UY #print 'hello' if IRET==1: print ' There is a problem in Okada subroutine!' break # # transform from Okada's to Aki's system # disp0[0]=UX*csst+UY*ssst disp0[1]=UX*ssst-UY*csst disp0[2]=-UZ tilt0[0]=-(UXZ*csst+UYZ*ssst) tilt0[1]=-(UXZ*ssst-UYZ*csst) # strain0[0]=(UXX*csst*csst+UYY*ssst*ssst +0.5*(UXY+UYX)*ss2st) strain0[1]=(UXX*ssst*ssst+UYY*csst*csst -0.5*(UXY+UYX)*ss2st) strain0[2]=(UZZ) strain0[3]=(0.5*(UXX-UYY)*ss2st -0.5*(UXY+UYX)*cs2st) strain0[4]=(-0.5*(UZX+UXZ)*ssst +0.5*(UYZ+UZY)*csst) strain0[5]=(-0.5*(UZX+UXZ)*csst -0.5*(UYZ+UZY)*ssst) for j in range(0,3): disp[irec][j]= disp[irec][j] + disp0[j] for j in range(0,2): tilt[irec][j]= tilt[irec][j] + tilt0[j] for j in range(0,6): strain[irec][j]= strain[irec][j] + strain0[j] z.append(-disp[irec][2]) return z # works on nautilus when have already done # f2py -c okada.f -m okada #z1 = okada(xr,yr,depth,length,width,dip,rake,slip) def DC3D0(self,ALPHA,X,Y,Z,DEPTH,DIP,POT1,POT2,POT3,POT4): """******************************************************************** ***** ***** #***** DISPLACEMENT AND STRAIN AT DEPTH ***** #***** DUE TO BURIED POINT SOURCE IN A SEMIINFINITE MEDIUM ***** #***** CODED BY Y.OKADA ... SEP.1991 ***** #***** REVISED Y.OKADA ... NOV.1991 ***** #***** ***** #******************************************************************** #***** INPUT #***** ALPHA : MEDIUM CONSTANT (LAMBDA+MYU)/(LAMBDA+2*MYU) #***** X,Y,Z : COORDINATE OF OBSERVING POINT #***** DEPTH : SOURCE DEPTH #***** DIP : DIP-ANGLE (DEGREE) #***** POT1-POT4 : STRIKE-, DIP-, TENSILE- AND INFLATE-POTENCY #***** POTENCY=( MOMENT OF DOUBLE-COUPLE )/MYU FOR POT1,2 #***** POTENCY=(INTENSITY OF ISOTROPIC PART)/LAMBDA FOR POT3 #***** POTENCY=(INTENSITY OF LINEAR DIPOLE )/MYU FOR POT4 #***** OUTPUT #***** UX, UY, UZ : DISPLACEMENT ( UNIT=(UNIT OF POTENCY) / #***** : (UNIT OF X,Y,Z,DEPTH)**2 ) #***** UXX,UYX,UZX : X-DERIVATIVE ( UNIT= UNIT OF POTENCY) / #***** UXY,UYY,UZY : Y-DERIVATIVE (UNIT OF X,Y,Z,DEPTH)**3 ) #***** UXZ,UYZ,UZZ : Z-DERIVATIVE #***** IRET : RETURN CODE ( =0....NORMAL, =1....SINGULAR )""" # COMMON /C1/DUMMY(8),R,dumm(15) # DIMENSION U(12),DUA(12),DUB(12),DUC(12) F0=0.0 U=num.zeros((12,1),num.float) DUA=num.zeros((12,1),num.float) DUB=num.zeros((12,1),num.float) DUC=num.zeros((12,1),num.float) if Z>0: print'(''0** POSITIVE Z WAS GIVEN IN SUB-DC3D0'')' for I in range(0,12): U[I]=F0 DUA[I]=F0 DUB[I]=F0 DUC[I]=F0 AALPHA=ALPHA DDIP=DIP self.DCC0N0(AALPHA,DDIP) #====================================== #===== REAL-SOURCE CONTRIBUTION ===== #====================================== XX=X YY=Y ZZ=Z DD=DEPTH+Z self.DCCON1(XX,YY,DD) R=self.R if R==F0: UX=F0 UY=F0 UZ=F0 UXX=F0 UYX=F0 UZX=F0 UXY=F0 UYY=F0 UZY=F0 UXZ=F0 UYZ=F0 UZZ=F0 IRET=1 return else: PP1=POT1 PP2=POT2 PP3=POT3 PP4=POT4 self.UA0(XX,YY,DD,PP1,PP2,PP3,PP4) DUA=self.DUA for I in range(0,12): if I<10: U[I]=U[I]-DUA[I] if I>=10: U[I]=U[I]+DUA[I] #======================================= #===== IMAGE-SOURCE CONTRIBUTION ===== #======================================= DD=DEPTH-Z self.DCCON1(XX,YY,DD) self.UA0(XX,YY,DD,PP1,PP2,PP3,PP4) DUA=self.DUA self.UB0(XX,YY,DD,ZZ,PP1,PP2,PP3,PP4) DUB=self.DUB self.UC0(XX,YY,DD,ZZ,PP1,PP2,PP3,PP4) DUC=self.DUC #----- for I in range(0,12): DU=DUA[I]+DUB[I]+ZZ*DUC[I] if I>=9: DU=DU+DUC[I-9] U[I]=U[I]+DU #===== self.UX=U[0] self.UY=U[1] self.UZ=U[2] self.UXX=U[3] self.UYX=U[4] self.UZX=U[5] self.UXY=U[6] self.UYY=U[7] self.UZY=U[8] self.UXZ=U[9] self.UYZ=U[10] self.UZZ=U[11] self.IRET=0 def UA0(self,X,Y,D,POT1,POT2,POT3,POT4): #SUBROUTINE UA0(X,Y,D,POT1,POT2,POT3,POT4,U) # IMPLICIT REAL*8 (A-H,O-Z) # DIMENSION U(12),DU(12) """******************************************************************** C***** DISPLACEMENT AND STRAIN AT DEPTH (PART-A) ***** C***** DUE TO BURIED POINT SOURCE IN A SEMIINFINITE MEDIUM ***** C******************************************************************** C***** INPUT C***** X,Y,D : STATION COORDINATES IN FAULT SYSTEM C***** POT1-POT4 : STRIKE-, DIP-, TENSILE- AND INFLATE-POTENCY C***** OUTPUT C***** U(12) : DISPLACEMENT AND THEIR DERIVATIVES """ # # COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D # COMMON /C1/P,Q,S,T,XY,X2,Y2,D2,R,R2,R3,R5,QR,QRX,A3,A5,B3,C3, # * UY,VY,WY,UZ,VZ,WZ F0 = 0.0 F1 = 1.0 F3 = 3.0 PI2=6.283185307179586 ALP1=self.ALP1 ALP2=self.ALP2 ALP3=self.ALP3 ALP4=self.ALP4 ALP5=self.ALP4 SD=self.SD CD=self.CD SDSD=self.SDSD CDCD=self.CDCD SDCD=self.SDCD S2D=self.S2D C2D=self.C2D P=self.P Q=self.Q S=self.S T=self.T XY=self.XY X2=self.X2 Y2=self.Y2 D2=self.D2 R=self.R R2=self.R2 R3=self.R3 R5=self.R5 QR=self.QR QRX=self.QRX A3=self.A3 A5=self.A5 B3=self.B3 C3=self.C3 UY=self.UY VY=self.VY WY=self.WY UZ=self.UZ VZ=self.VZ WZ=self.WZ DUA=num.zeros((12,1),num.float) DU=num.zeros((12,1),num.float) U=num.zeros((12,1),num.float) #----- for I in range(0,12): U[I]=F0 #====================================== #===== STRIKE-SLIP CONTRIBUTION ===== #====================================== if POT1 != F0: DU[0]= ALP1*Q/R3 +ALP2*X2*QR DU[1]= ALP1*X/R3*SD +ALP2*XY*QR DU[2]=-ALP1*X/R3*CD +ALP2*X*D*QR DU[3]= X*QR*(-ALP1 +ALP2*(F1+A5) ) DU[4]= ALP1*A3/R3*SD +ALP2*Y*QR*A5 DU[5]=-ALP1*A3/R3*CD +ALP2*D*QR*A5 DU[6]= ALP1*(SD/R3-Y*QR) +ALP2*F3*X2/R5*UY DU[7]= F3*X/R5*(-ALP1*Y*SD +ALP2*(Y*UY+Q) ) DU[8]= F3*X/R5*( ALP1*Y*CD +ALP2*D*UY ) DU[9]= ALP1*(CD/R3+D*QR) +ALP2*F3*X2/R5*UZ DU[10]= F3*X/R5*( ALP1*D*SD +ALP2*Y*UZ ) DU[11]= F3*X/R5*(-ALP1*D*CD +ALP2*(D*UZ-Q) ) for I in range(0,12): U[I]=U[I]+POT1/PI2*DU[I] #=================================== #===== DIP-SLIP CONTRIBUTION ===== #=================================== if POT2 != F0: DU[0]= ALP2*X*P*QR DU[1]= ALP1*S/R3 +ALP2*Y*P*QR DU[2]=-ALP1*T/R3 +ALP2*D*P*QR DU[3]= ALP2*P*QR*A5 DU[4]=-ALP1*F3*X*S/R5 -ALP2*Y*P*QRX DU[5]= ALP1*F3*X*T/R5 -ALP2*D*P*QRX DU[6]= ALP2*F3*X/R5*VY DU[7]= ALP1*(S2D/R3-F3*Y*S/R5) +ALP2*(F3*Y/R5*VY+P*QR) DU[8]=-ALP1*(C2D/R3-F3*Y*T/R5) +ALP2*F3*D/R5*VY DU[9]= ALP2*F3*X/R5*VZ DU[10]= ALP1*(C2D/R3+F3*D*S/R5) +ALP2*F3*Y/R5*VZ DU[11]= ALP1*(S2D/R3-F3*D*T/R5) +ALP2*(F3*D/R5*VZ-P*QR) for I in range(0,12): U[I]=U[I]+POT2/PI2*DU[I] #======================================== #===== TENSILE-FAULT CONTRIBUTION ===== #======================================== if POT3!=F0: DU[0]= ALP1*X/R3 -ALP2*X*Q*QR DU[1]= ALP1*T/R3 -ALP2*Y*Q*QR DU[2]= ALP1*S/R3 -ALP2*D*Q*QR DU[3]= ALP1*A3/R3 -ALP2*Q*QR*A5 DU[4]=-ALP1*F3*X*T/R5 +ALP2*Y*Q*QRX DU[5]=-ALP1*F3*X*S/R5 +ALP2*D*Q*QRX DU[6]=-ALP1*F3*XY/R5 -ALP2*X*QR*WY DU[7]= ALP1*(C2D/R3-F3*Y*T/R5) -ALP2*(Y*WY+Q)*QR DU[8]= ALP1*(S2D/R3-F3*Y*S/R5) -ALP2*D*QR*WY DU[9]= ALP1*F3*X*D/R5 -ALP2*X*QR*WZ DU[10]=-ALP1*(S2D/R3-F3*D*T/R5) -ALP2*Y*QR*WZ DU[11]= ALP1*(C2D/R3+F3*D*S/R5) -ALP2*(D*WZ-Q)*QR for I in range(0,12): U[I]=U[I]+POT3/PI2*DU[I] #========================================= #===== INFLATE SOURCE CONTRIBUTION ===== #========================================= if POT4 != F0: DU[0]=-ALP1*X/R3 DU[1]=-ALP1*Y/R3 DU[2]=-ALP1*D/R3 DU[3]=-ALP1*A3/R3 DU[4]= ALP1*F3*XY/R5 DU[5]= ALP1*F3*X*D/R5 DU[6]= DU[4] DU[7]=-ALP1*B3/R3 DU[8]= ALP1*F3*Y*D/R5 DU[9]=-DU[5] DU[10]=-DU[8] DU[11]= ALP1*C3/R3 for I in range(0,12): U[I]=U[I]+POT4/PI2*DU[I] #for I in range(0,12): #DUA[I]=U[I] self.DUA=U def UB0(self,X,Y,D,Z,POT1,POT2,POT3,POT4): # SUBROUTINE UB0(X,Y,D,Z,POT1,POT2,POT3,POT4,U) # IMPLICIT REAL*8 (A-H,O-Z) # DIMENSION U(12),DU(12) # """******************************************************************** C***** DISPLACEMENT AND STRAIN AT DEPTH (PART-B) ***** C***** DUE TO BURIED POINT SOURCE IN A SEMIINFINITE MEDIUM ***** C******************************************************************** C C***** INPUT C***** X,Y,D,Z : STATION COORDINATES IN FAULT SYSTEM C***** POT1-POT4 : STRIKE-, DIP-, TENSILE- AND INFLATE-POTENCY C***** OUTPUT C***** U(12) : DISPLACEMENT AND THEIR DERIVATIVES """ # # COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D # COMMON /C1/P,Q,S,T,XY,X2,Y2,D2,R,R2,R3,R5,QR,QRX,A3,A5,B3,C3, # * UY,VY,WY,UZ,VZ,WZ # F0,F1,F2,F3,F4,F5,F8,F9 # * /0.D0,1.D0,2.D0,3.D0,4.D0,5.D0,8.D0,9.D0/ # DATA PI2/6.283185307179586D0/ DUB=num.zeros((12,1),num.float) DU=num.zeros((12,1),num.float) U=num.zeros((12,1),num.float) 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 ALP1=self.ALP1 ALP2=self.ALP2 ALP3=self.ALP3 ALP4=self.ALP4 ALP5=self.ALP4 SD=self.SD CD=self.CD SDSD=self.SDSD CDCD=self.CDCD SDCD=self.SDCD S2D=self.S2D C2D=self.C2D P=self.P Q=self.Q S=self.S T=self.T XY=self.XY X2=self.X2 Y2=self.Y2 D2=self.D2 R=self.R R2=self.R2 R3=self.R3 R5=self.R5 QR=self.QR QRX=self.QRX A3=self.A3 A5=self.A5 B3=self.B3 C3=self.C3 UY=self.UY VY=self.VY WY=self.WY UZ=self.UZ VZ=self.VZ WZ=self.WZ #----- C=D+Z RD=R+D D12=F1/(R*RD*RD) D32=D12*(F2*R+D)/R2 D33=D12*(F3*R+D)/(R2*RD) D53=D12*(F8*R2+F9*R*D+F3*D2)/(R2*R2*RD) D54=D12*(F5*R2+F4*R*D+D2)/R3*D12 #----- FI1= Y*(D12-X2*D33) FI2= X*(D12-Y2*D33) FI3= X/R3-FI2 FI4=-XY*D32 FI5= F1/(R*RD)-X2*D32 FJ1=-F3*XY*(D33-X2*D54) FJ2= F1/R3-F3*D12+F3*X2*Y2*D54 FJ3= A3/R3-FJ2 FJ4=-F3*XY/R5-FJ1 FK1=-Y*(D32-X2*D53) FK2=-X*(D32-Y2*D53) FK3=-F3*X*D/R5-FK2 #----- for I in range(0,12): U[I]=F0 #====================================== #===== STRIKE-SLIP CONTRIBUTION ===== #====================================== if POT1!=F0: DU[0]=-X2*QR -ALP3*FI1*SD DU[1]=-XY*QR -ALP3*FI2*SD DU[2]=-C*X*QR -ALP3*FI4*SD DU[3]=-X*QR*(F1+A5) -ALP3*FJ1*SD DU[4]=-Y*QR*A5 -ALP3*FJ2*SD DU[5]=-C*QR*A5 -ALP3*FK1*SD DU[6]=-F3*X2/R5*UY -ALP3*FJ2*SD DU[7]=-F3*XY/R5*UY-X*QR -ALP3*FJ4*SD DU[8]=-F3*C*X/R5*UY -ALP3*FK2*SD DU[9]=-F3*X2/R5*UZ +ALP3*FK1*SD DU[10]=-F3*XY/R5*UZ +ALP3*FK2*SD DU[11]= F3*X/R5*(-C*UZ +ALP3*Y*SD) for I in range(0,12): U[I]=U[I]+POT1/PI2*DU[I] #=================================== #===== DIP-SLIP CONTRIBUTION ===== #=================================== if POT2!=F0: DU[0]=-X*P*QR +ALP3*FI3*SDCD DU[1]=-Y*P*QR +ALP3*FI1*SDCD DU[2]=-C*P*QR +ALP3*FI5*SDCD DU[3]=-P*QR*A5 +ALP3*FJ3*SDCD DU[4]= Y*P*QRX +ALP3*FJ1*SDCD DU[5]= C*P*QRX +ALP3*FK3*SDCD DU[6]=-F3*X/R5*VY +ALP3*FJ1*SDCD DU[7]=-F3*Y/R5*VY-P*QR +ALP3*FJ2*SDCD DU[8]=-F3*C/R5*VY +ALP3*FK1*SDCD DU[9]=-F3*X/R5*VZ -ALP3*FK3*SDCD DU[10]=-F3*Y/R5*VZ -ALP3*FK1*SDCD DU[11]=-F3*C/R5*VZ +ALP3*A3/R3*SDCD for I in range(0,12): U[I]=U[I]+POT2/PI2*DU[I] #======================================== #===== TENSILE-FAULT CONTRIBUTION ===== #======================================== if POT3!=F0: DU[0]= X*Q*QR -ALP3*FI3*SDSD DU[1]= Y*Q*QR -ALP3*FI1*SDSD DU[2]= C*Q*QR -ALP3*FI5*SDSD DU[3]= Q*QR*A5 -ALP3*FJ3*SDSD DU[4]=-Y*Q*QRX -ALP3*FJ1*SDSD DU[5]=-C*Q*QRX -ALP3*FK3*SDSD DU[6]= X*QR*WY -ALP3*FJ1*SDSD DU[7]= QR*(Y*WY+Q) -ALP3*FJ2*SDSD DU[8]= C*QR*WY -ALP3*FK1*SDSD DU[9]= X*QR*WZ +ALP3*FK3*SDSD DU[10]= Y*QR*WZ +ALP3*FK1*SDSD DU[11]= C*QR*WZ -ALP3*A3/R3*SDSD for I in range(0,12): U[I]=U[I]+POT3/PI2*DU[I] #========================================= #===== INFLATE SOURCE CONTRIBUTION ===== #========================================= if POT4!=F0: DU[0]= ALP3*X/R3 DU[1]= ALP3*Y/R3 DU[2]= ALP3*D/R3 DU[3]= ALP3*A3/R3 DU[4]=-ALP3*F3*XY/R5 DU[5]=-ALP3*F3*X*D/R5 DU[6]= DU[4] DU[7]= ALP3*B3/R3 DU[8]=-ALP3*F3*Y*D/R5 DU[9]=-DU[5] DU[10]=-DU[8] DU[11]=-ALP3*C3/R3 for I in range(0,12): U[I]=U[I]+POT4/PI2*DU[I] #for I in range(0,12): #DUB[I]=U[I] self.DUB=U def UC0(self,X,Y,D,Z,POT1,POT2,POT3,POT4): # SUBROUTINE UC0(X,Y,D,Z,POT1,POT2,POT3,POT4,U) # IMPLICIT REAL*8 (A-H,O-Z) # DIMENSION U(12),DU(12) """******************************************************************** C***** DISPLACEMENT AND STRAIN AT DEPTH (PART-B) ***** C***** DUE TO BURIED POINT SOURCE IN A SEMIINFINITE MEDIUM ***** C******************************************************************** C***** INPUT C***** X,Y,D,Z : STATION COORDINATES IN FAULT SYSTEM C***** POT1-POT4 : STRIKE-, DIP-, TENSILE- AND INFLATE-POTENCY C***** OUTPUT C***** U(12) : DISPLACEMENT AND THEIR DERIVATIVES""" # COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D # COMMON /C1/P,Q,S,T,XY,X2,Y2,D2,R,R2,R3,R5,QR,QRX,A3,A5,B3,C3,um(6) # DATA F0,F1,F2,F3,F5,F7,F10,F15 # * /0.D0,1.D0,2.D0,3.D0,5.D0,7.D0,10.D0,15.D0/ # DATA PI2/6.283185307179586D0/ DUC=num.zeros((12,1),num.float) DU=num.zeros((12,1),num.float) U=num.zeros((12,1),num.float) F0=0.0 F1=1.0 F2=2.0 F3=3.0 F5=5.0 F7=7.0 F10=10.0 F15=15.0 PI2=6.283185307179586 ALP1=self.ALP1 ALP2=self.ALP2 ALP3=self.ALP3 ALP4=self.ALP4 ALP5=self.ALP4 SD=self.SD CD=self.CD SDSD=self.SDSD CDCD=self.CDCD SDCD=self.SDCD S2D=self.S2D C2D=self.C2D P=self.P Q=self.Q S=self.S T=self.T XY=self.XY X2=self.X2 Y2=self.Y2 D2=self.D2 R=self.R R2=self.R2 R3=self.R3 R5=self.R5 QR=self.QR QRX=self.QRX A3=self.A3 A5=self.A5 B3=self.B3 C3=self.C3 UY=self.UY VY=self.VY WY=self.WY UZ=self.UZ VZ=self.VZ WZ=self.WZ #----- C=D+Z Q2=Q*Q R7=R5*R2 A7=F1-F7*X2/R2 B5=F1-F5*Y2/R2 B7=F1-F7*Y2/R2 C5=F1-F5*D2/R2 C7=F1-F7*D2/R2 D7=F2-F7*Q2/R2 QR5=F5*Q/R2 QR7=F7*Q/R2 DR5=F5*D/R2 #----- for I in range(0,12): U[I]=F0 #====================================== #===== STRIKE-SLIP CONTRIBUTION ===== #====================================== if POT1!=F0: DU[0]=-ALP4*A3/R3*CD +ALP5*C*QR*A5 DU[1]= F3*X/R5*( ALP4*Y*CD +ALP5*C*(SD-Y*QR5) ) DU[2]= F3*X/R5*(-ALP4*Y*SD +ALP5*C*(CD+D*QR5) ) DU[3]= ALP4*F3*X/R5*(F2+A5)*CD -ALP5*C*QRX*(F2+A7) DU[4]= F3/R5*( ALP4*Y*A5*CD +ALP5*C*(A5*SD-Y*QR5*A7) ) DU[5]= F3/R5*(-ALP4*Y*A5*SD +ALP5*C*(A5*CD+D*QR5*A7) ) DU[6]= DU[4] DU[7]= F3*X/R5*( ALP4*B5*CD -ALP5*F5*C/R2*(F2*Y*SD+Q*B7) ) DU[8]= F3*X/R5*(-ALP4*B5*SD +ALP5*F5*C/R2*(D*B7*SD-Y*C7*CD) ) DU[9]= F3/R5* (-ALP4*D*A5*CD +ALP5*C*(A5*CD+D*QR5*A7) ) DU[10]= F15*X/R7*( ALP4*Y*D*CD +ALP5*C*(D*B7*SD-Y*C7*CD) ) DU[11]= F15*X/R7*(-ALP4*Y*D*SD +ALP5*C*(F2*D*CD-Q*C7) ) for I in range(0,12): U[I]=U[I]+POT1/PI2*DU[I] #=================================== #===== DIP-SLIP CONTRIBUTION ===== #=================================== if POT2!=F0: DU[0]= ALP4*F3*X*T/R5 -ALP5*C*P*QRX DU[1]=-ALP4/R3*(C2D-F3*Y*T/R2) +ALP5*F3*C/R5*(S-Y*P*QR5) DU[2]=-ALP4*A3/R3*SDCD +ALP5*F3*C/R5*(T+D*P*QR5) DU[3]= ALP4*F3*T/R5*A5 -ALP5*F5*C*P*QR/R2*A7 DU[4]= F3*X/R5*(ALP4*(C2D-F5*Y*T/R2)-ALP5*F5*C/R2*(S-Y*P*QR7)) DU[5]= F3*X/R5*(ALP4*(F2+A5)*SDCD -ALP5*F5*C/R2*(T+D*P*QR7)) DU[6]= DU[4] DU[7]= (F3/R5*(ALP4*(F2*Y*C2D+T*B5) +ALP5*C*(S2D-F10*Y*S/R2-P*QR5*B7))) DU[8]= F3/R5*(ALP4*Y*A5*SDCD-ALP5*C*((F3+A5)*C2D+Y*P*DR5*QR7)) DU[9]= F3*X/R5*(-ALP4*(S2D-T*DR5) -ALP5*F5*C/R2*(T+D*P*QR7)) DU[10]= (F3/R5*(-ALP4*(D*B5*C2D+Y*C5*S2D) -ALP5*C*((F3+A5)*C2D+Y*P*DR5*QR7))) DU[11]= F3/R5*(-ALP4*D*A5*SDCD-ALP5*C*(S2D-F10*D*T/R2+P*QR5*C7)) for I in range(0,12): U[I]=U[I]+POT2/PI2*DU[I] #======================================== #===== TENSILE-FAULT CONTRIBUTION ===== #======================================== if POT3!=F0: DU[0]= F3*X/R5*(-ALP4*S +ALP5*(C*Q*QR5-Z)) DU[1]= ALP4/R3*(S2D-F3*Y*S/R2)+ALP5*F3/R5*(C*(T-Y+Y*Q*QR5)-Y*Z) DU[2]=-ALP4/R3*(F1-A3*SDSD) -ALP5*F3/R5*(C*(S-D+D*Q*QR5)-D*Z) DU[3]=-ALP4*F3*S/R5*A5 +ALP5*(C*QR*QR5*A7-F3*Z/R5*A5) DU[4]= (F3*X/R5*(-ALP4*(S2D-F5*Y*S/R2) -ALP5*F5/R2*(C*(T-Y+Y*Q*QR7)-Y*Z))) DU[5]= (F3*X/R5*( ALP4*(F1-(F2+A5)*SDSD) +ALP5*F5/R2*(C*(S-D+D*Q*QR7)-D*Z))) DU[6]= DU[4] DU[7]= (F3/R5*(-ALP4*(F2*Y*S2D+S*B5) -ALP5*(C*(F2*SDSD+F10*Y*(T-Y)/R2-Q*QR5*B7)+Z*B5))) DU[8]= (F3/R5*( ALP4*Y*(F1-A5*SDSD) +ALP5*(C*(F3+A5)*S2D-Y*DR5*(C*D7+Z)))) DU[9]= (F3*X/R5*(-ALP4*(C2D+S*DR5) +ALP5*(F5*C/R2*(S-D+D*Q*QR7)-F1-Z*DR5))) DU[10]= (F3/R5*( ALP4*(D*B5*S2D-Y*C5*C2D) +ALP5*(C*((F3+A5)*S2D-Y*DR5*D7)-Y*(F1+Z*DR5)))) DU[11]= (F3/R5*(-ALP4*D*(F1-A5*SDSD) -ALP5*(C*(C2D+F10*D*(S-D)/R2-Q*QR5*C7)+Z*(F1+C5)))) for I in range(0,12): U[I]=U[I]+POT3/PI2*DU[I] #========================================= #===== INFLATE SOURCE CONTRIBUTION ===== #========================================= if POT4!=F0: DU[0]= ALP4*F3*X*D/R5 DU[1]= ALP4*F3*Y*D/R5 DU[2]= ALP4*C3/R3 DU[3]= ALP4*F3*D/R5*A5 DU[4]=-ALP4*F15*XY*D/R7 DU[5]=-ALP4*F3*X/R5*C5 DU[6]= DU[4] DU[7]= ALP4*F3*D/R5*B5 DU[8]=-ALP4*F3*Y/R5*C5 DU[9]= DU[5] DU[10]= DU[8] DU[11]= ALP4*F3*D/R5*(F2+C5) for I in range(0,12): U[I]=U[I]+POT4/PI2*DU[I] #for I in range(0,12): #DUC[I]=U[I] self.DUC=U def DC3D(self,ALPHA,X,Y,Z,DEPTH,DIP,AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3): # SUBROUTINE DC3D(ALPHA,X,Y,Z,DEPTH,DIP, # * AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3, # * UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,IRET) # IMPLICIT REAL*8 (A-H,O-Z) # REAL*4 ALPHA,X,Y,Z,DEPTH,DIP,AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3, # * UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,EPS """******************************************************************** C***** ***** C***** DISPLACEMENT AND STRAIN AT DEPTH ***** C***** DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM ***** C***** CODED BY Y.OKADA ... SEP.1991 ***** C***** REVISED ... NOV.1991, APR.1992, MAY.1993, ***** C***** JUL.1993 ***** C******************************************************************** C C***** INPUT C***** ALPHA : MEDIUM CONSTANT (LAMBDA+MYU)/(LAMBDA+2*MYU) C***** X,Y,Z : COORDINATE OF OBSERVING POINT C***** DEPTH : DEPTH OF REFERENCE POINT C***** DIP : DIP-ANGLE (DEGREE) C***** AL1,AL2 : FAULT LENGTH RANGE C***** AW1,AW2 : FAULT WIDTH RANGE C***** DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS C C***** OUTPUT C***** UX, UY, UZ : DISPLACEMENT ( UNIT=(UNIT OF DISL) C***** UXX,UYX,UZX : X-DERIVATIVE ( UNIT=(UNIT OF DISL) / C***** UXY,UYY,UZY : Y-DERIVATIVE (UNIT OF X,Y,Z,DEPTH,AL,AW) ) C***** UXZ,UYZ,UZZ : Z-DERIVATIVE C***** IRET : RETURN CODE ( =0....NORMAL, =1....SINGULAR ) """ # # COMMON /C0/DUMMY(5),SD,CD,dumm(5) # DIMENSION XI(2),ET(2),KXI(2),KET(2) # DIMENSION U(12),DU(12),DUA(12),DUB(12),DUC(12) from math import sqrt F0 =0.0 EPS=0.000001 XI=num.zeros(2,num.float) ET=num.zeros(2,num.float) KXI=num.zeros(2,num.float) KET=num.zeros(2,num.float) U=num.zeros(12,num.float) DU=num.zeros(12,num.float) DUA=num.zeros(12,num.float) DUB=num.zeros(12,num.float) DUC=num.zeros(12,num.float) #----- if Z>0: print '** POSITIVE Z WAS GIVEN IN SUB-DC3D' for I in range(0,12): U[I]=F0 DUA[I]=F0 DUB[I]=F0 DUC[I]=F0 AALPHA=ALPHA DDIP=DIP self.DCC0N0(AALPHA,DDIP) CD=self.CD SD=self.SD #----- ZZ=Z DD1=DISL1 DD2=DISL2 DD3=DISL3 XI[0]=X-AL1 XI[1]=X-AL2 #if X==-6 and Y==-1: #print AL1 #print AL2 #print 'hello' if abs(XI[0])F0: SD= F1 if SD