Last change on this file since 5421 was 5252, checked in by herve, 15 years ago

new definition for zrec
change sign for displacement

File size: 99.6 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 (km)
12 width   down-dip width of rupture area (km)
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 (km)
16 ns      the really used number of rectangular sources
17 NSMAX   the upper limit of ns
18Optional
19 x0      x origin (0)
20 y0      y origin (0)
21 z0      z origin
22 slip    metres of fault slip (1)
23 rake    angle of slip (w.r.t. horizontal) in fault plane (90 degrees)
24
25
26
27The returned object is a callable okada function that represents
28the initial water displacement generated by a submarine earthuake.
29
30"""
31
32def earthquake_tsunami(ns,NSMAX,length, width, strike, depth,\
33                       dip, xi, yi,z0, slip, rake,\
34                       domain=None, verbose=False):
35
36    from anuga.abstract_2d_finite_volumes.quantity import Quantity
37    from math import sin, radians
38    from Numeric import zeros, Float
39
40    #zrec0 = Quantity(domain)
41    zrec0 = Quantity(domain)
42    zrec0.set_values(z0)
43    zrec=zrec0.get_vertex_values(xy=True)
44
45    x0= zeros(ns,Float)
46    y0= zeros(ns,Float)
47    if ns ==1:
48        x0[0]=xi
49        y0[0]=yi
50    else:
51        x0=xi
52        y0=yi
53    if domain is not None:
54        xllcorner = domain.geo_reference.get_xllcorner()
55        yllcorner = domain.geo_reference.get_yllcorner()
56        for i in range(0,ns):
57            x0[i] = x0[i] - xllcorner  # fault origin (relative)
58            y0[i] = y0[i] - yllcorner
59    #a few temporary print statements
60    if verbose is True:
61        print '\nThe Earthquake ...'
62        print '\tns: ', ns
63        print '\tNSMAX: ', NSMAX
64        print '\tLength: ', length
65        print '\tDepth: ', depth
66        print '\tStrike: ', strike
67        print '\tWidth: ', width
68        print '\tDip: ', dip
69        print '\tSlip: ', slip
70        print '\tx0: ', x0
71        print '\ty0: ', y0
72
73    # warning state
74#    test = width*1000.0*sin(radians(dip)) - depth
75    test = width*1000.0*sin(radians(dip)) - depth*1000
76
77    if test > 0.0:
78        msg = 'Earthquake source not located below seafloor - check depth'
79        raise Exception, msg
80
81    return Okada_func(ns=ns,NSMAX=NSMAX,length=length, width=width, dip=dip, \
82                      x0=x0, y0=y0, strike=strike, depth=depth, \
83                      slip=slip, rake=rake, zrec=zrec)
84
85#
87#
88
89"""This is a callable class representing the initial water displacment
90   generated by an earthquake.
91
92Using input parameters:
93
94Required
95 length  along-stike length of rupture area
96 width   down-dip width of rupture area
97 strike  azimuth (degrees, measured from north) of fault axis
98 dip     angle of fault dip in degrees w.r.t. horizontal
99 depth   depth to base of rupture area
100
101Optional
102 x0      x origin (0)
103 y0      y origin (0)
104 slip    metres of fault slip (1)
105 rake    angle of slip (w.r.t. horizontal) in fault plane (90 degrees)
106
107"""
108
110
111    def __init__(self, ns,NSMAX,length, width, dip, x0, y0, strike, \
112                     depth, slip, rake,zrec):
113        self.dip = dip
114        self.length = length
115        self.width = width
116        self.x0 = x0
117        self.y0 = y0
118        self.strike = strike
119        self.depth = depth
120        self.slip = slip
121        self.rake = rake
122        self.ns=ns
123        self.zrec=zrec
124
125    def __call__(self, x, y):
126        """Make Okada_func a callable object.
127
128        If called as a function, this object returns z values representing
129        the initial 3D distribution of water heights at the points (x,y,z)
130        produced by a submarine mass failure.
131        """
132        from Numeric import zeros, Float
133        from string import replace,strip
134        from math import sin, cos, radians, exp, cosh
135        from Numeric import zeros, Float
136        #ensure vectors x and y have the same length
137        N = len(x)
138        assert N == len(y)
139        #nrec=N*N
140        depth = self.depth
141        dip = self.dip
142        length = self.length
143        width = self.width
144        y0 = self.x0
145        x0 = self.y0
146        strike = self.strike
147        dip = self.dip
148        rake = self.rake
149        dislocation = self.slip
150        ns=self.ns
151        zrec=self.zrec
152        #initialization
153        disp0=zeros(3,Float)
154        strain0=zeros(6,Float)
155        tilt0 = zeros(2,Float)
156        dislocations=zeros(ns,Float)
157        depths=zeros(ns,Float)
158        strikes= zeros(ns,Float)
159        lengths= zeros(ns,Float)
160        slips= zeros(ns,Float)
161        rakes= zeros(ns,Float)
162        widths= zeros(ns,Float)
163        dips= zeros(ns,Float)
164        strikes= zeros(ns,Float)
165        strikes= zeros(ns,Float)
166        strain = zeros((N,6),Float)
167        disp = zeros((N,3),Float)
168        tilt = zeros((N,2),Float)
169        xs =zeros(ns,Float)
170        ys =zeros(ns,Float)
171        z=[]
172        if ns==1:
173            dislocations[0]=dislocation
174            depths[0]=depth
175            strikes[0]=strike
176            lengths[0]=length
177            rakes[0]=rake
178            widths[0]=width
179            dips[0]=dip
180            try:
181                xs[0]=x0
182                ys[0]=y0
183            except:
184                xs[0]=x0[0]
185                ys[0]=y0[0]
186        else:
187            dislocations=dislocation
188            strikes=strike
189            lengths=length
190            rakes=rake
191            widths=width
192            dips=dip
193            xs=x0
194            ys=y0
195            depths=depth
196        #double Gaussian calculation assumes water displacement is oriented
197        #E-W, so, for displacement at some angle alpha clockwise from the E-W
198        #direction, rotate (x,y) coordinates anti-clockwise by alpha
199        ALPHA = 0.5
200        POT3 = 0.0
201        POT4 = 0.0
202        DISL3 = 0.0
203        AL1 = 0.0
204        AW2 = 0.0
206        #zrec = domain.get_quantity('elevation').get_values(interpolation_points=[[x, y]])   #get Zrec... has to be constant!!!
207        #zrec=zreci.get_values(interpolation_points=[[x, y]])
208        #zrec=0
209        #Z=-zrec
210        eps = 1.0e-6
211
212#
213        for irec in range(0,N):
214            xrec=y
215            yrec=x
216            for i in range(0,len(zrec[0])):
217                if zrec[0][i]==yrec and zrec[1][i]==xrec:
218                    Z=zrec[2][i]
219                    Z=0.001*Z
220                    break
221                else: continue
222            #zrec=zreci.get_values(interpolation_points=[[xrec[irec],yrec[irec]]],location='edges')
223
224            for ist in range(0,ns):
227                csst = cos(st)
228                ssst = sin(st)
229                cs2st = cos(2.0*st)
230                ss2st = sin(2.0*st)
231#
234                csdi =cos(di)
235                ssdi=sin(di)
236#
239                csra=cos(ra)
240                ssra=sin(ra)
241                # transform from Aki's to Okada's system
242                #first attribute x to north and y to east to match OKADA axis
243
244                #X=0.001*(x[irec]-xs[ist])*ssst+0.001*(y[irec]-ys[ist])*csst
245                #Y=-0.001*(x[irec]-xs[ist])*csst+0.001*(y[irec]-ys[ist])*ssst
246
247                X=0.001*((xrec[irec]-xs[ist])*csst+(yrec[irec]-ys[ist])*ssst)
248                Y=0.001*((xrec[irec]-xs[ist])*ssst-(yrec[irec]-ys[ist])*csst)
249                DEPTH=depths[ist]
250                DIP=dips[ist]
251                if lengths[ist]==0 and widths[ist]== 0 :
252#
253#             point source
254#
255                    POT1=dislocations[ist]*csra
256                    POT2=dislocations[ist]*ssra
257                    IRET=1
258                    self.DC3D0(ALPHA,X,Y,Z,DEPTH,DIP,POT1,POT2,POT3,POT4)
259                    UX=self.UX
260                    UY=self.UY
261                    UZ=self.UZ
262                    UXX=self.UXX
263                    UYX=self.UYX
264                    UZX=self.UZX
265                    UXY=self.UXY
266                    UYY=self.UYY
267                    UZY=self.UZY
268                    UXZ=self.UXZ
269                    UYZ=self.UYZ
270                    UZZ=self.UZZ
271                    IRET=self.IRET
272                    if IRET==1:
273                        print ' There is a problem in Okada subroutine!'
274                        break
275                else:
276#               finite source
277                    AL2=lengths[ist]
278                    AW1=-widths[ist]
279                    DISL1=dislocations[ist]*csra
280                    DISL2=dislocations[ist]*ssra
281                    #print 'lenghts[',ist,']',lenghts[ist]
282                    if lengths[ist]==0:
283                        #print 'beh alors'
284                        AL2=widths[ist]*eps
285                        DISL1=DISL1/AL2
286                        DISL2=DISL2/AL2
287                    elif widths[ist]==0.0:
288                        AW1=-lengths[ist]*eps
289                        DISL1=DISL1/(-AW1)
290                        DISL2=DISL2/(-AW1)
291                    IRET=1
292                    self.DC3D(ALPHA,X,Y,Z,DEPTH,DIP,AL1,AL2,AW1,AW2,\
293                             DISL1,DISL2,DISL3)
294                    UX=self.UX
295                    UY=self.UY
296                    UZ=self.UZ
297                    UXX=self.UXX
298                    UYX=self.UYX
299                    UZX=self.UZX
300                    UXY=self.UXY
301                    UYY=self.UYY
302                    UZY=self.UZY
303                    UXZ=self.UXZ
304                    UYZ=self.UYZ
305                    UZZ=self.UZZ
306                    IRET=self.IRET
307                    #if X==-6 and Y==-1:
308                        #print UY
309                        #print 'hello'
310                    if IRET==1:
311                        print ' There is a problem in Okada subroutine!'
312                        break
313#
314#           transform from Okada's to Aki's system
315#
316                disp0[0]=UX*csst+UY*ssst
317                disp0[1]=UX*ssst-UY*csst
318                disp0[2]=-UZ
319
320                tilt0[0]=-(UXZ*csst+UYZ*ssst)
321                tilt0[1]=-(UXZ*ssst-UYZ*csst)
322#
323                strain0[0]=(UXX*csst*csst+UYY*ssst*ssst
324                            +0.5*(UXY+UYX)*ss2st)
325                strain0[1]=(UXX*ssst*ssst+UYY*csst*csst
326                             -0.5*(UXY+UYX)*ss2st)
327                strain0[2]=(UZZ)
328                strain0[3]=(0.5*(UXX-UYY)*ss2st
329                             -0.5*(UXY+UYX)*cs2st)
330                strain0[4]=(-0.5*(UZX+UXZ)*ssst
331                             +0.5*(UYZ+UZY)*csst)
332                strain0[5]=(-0.5*(UZX+UXZ)*csst
333                              -0.5*(UYZ+UZY)*ssst)
334                for j in range(0,3):
335                    disp[irec][j]=   disp[irec][j]   + disp0[j]
336                for j in range(0,2):
337                    tilt[irec][j]=   tilt[irec][j]   + tilt0[j]
338                for j in range(0,6):
339                    strain[irec][j]= strain[irec][j] + strain0[j]
340            z.append(-disp[irec][2])
341        return z
342
343
344        # works on nautilus when have already done
347
348
349    def DC3D0(self,ALPHA,X,Y,Z,DEPTH,DIP,POT1,POT2,POT3,POT4):
350
351          """********************************************************************
352          *****                                                          *****
353          #*****    DISPLACEMENT AND STRAIN AT DEPTH                      *****
354          #*****    DUE TO BURIED POINT SOURCE IN A SEMIINFINITE MEDIUM   *****
355          #*****                         CODED BY  Y.OKADA ... SEP.1991   *****
356          #*****                         REVISED   Y.OKADA ... NOV.1991   *****
357          #*****                                                          *****
358          #********************************************************************
359          #***** INPUT
360          #*****   ALPHA : MEDIUM CONSTANT  (LAMBDA+MYU)/(LAMBDA+2*MYU)
361          #*****   X,Y,Z : COORDINATE OF OBSERVING POINT
362          #*****   DEPTH : SOURCE DEPTH
363          #*****   DIP   : DIP-ANGLE (DEGREE)
364          #*****   POT1-POT4 : STRIKE-, DIP-, TENSILE- AND INFLATE-POTENCY
365          #*****       POTENCY=(  MOMENT OF DOUBLE-COUPLE  )/MYU     FOR POT1,2
366          #*****       POTENCY=(INTENSITY OF ISOTROPIC PART)/LAMBDA  FOR POT3
367          #*****       POTENCY=(INTENSITY OF LINEAR DIPOLE )/MYU     FOR POT4
368          #***** OUTPUT
369          #*****   UX, UY, UZ  : DISPLACEMENT ( UNIT=(UNIT OF POTENCY) /
370          #*****               :                     (UNIT OF X,Y,Z,DEPTH)**2  )
371          #*****   UXX,UYX,UZX : X-DERIVATIVE ( UNIT= UNIT OF POTENCY) /
372          #*****   UXY,UYY,UZY : Y-DERIVATIVE        (UNIT OF X,Y,Z,DEPTH)**3  )
373          #*****   UXZ,UYZ,UZZ : Z-DERIVATIVE
374          #*****   IRET        : RETURN CODE  ( =0....NORMAL,   =1....SINGULAR )"""
375
376#        COMMON /C1/DUMMY(8),R,dumm(15)
377#          DIMENSION  U(12),DUA(12),DUB(12),DUC(12)
378          from Numeric import zeros, Float
379
380          F0=0.0
381          U=zeros((12,1),Float)
382          DUA=zeros((12,1),Float)
383          DUB=zeros((12,1),Float)
384          DUC=zeros((12,1),Float)
385
386
387          if Z>0: print'(''0** POSITIVE Z WAS GIVEN IN SUB-DC3D0'')'
388          for I in range(0,12):
389            U[I]=F0
390            DUA[I]=F0
391            DUB[I]=F0
392            DUC[I]=F0
393
394          AALPHA=ALPHA
395          DDIP=DIP
396          self.DCC0N0(AALPHA,DDIP)
397#======================================
398#=====  REAL-SOURCE CONTRIBUTION  =====
399#======================================
400          XX=X
401          YY=Y
402          ZZ=Z
403          DD=DEPTH+Z
404          self.DCCON1(XX,YY,DD)
405          R=self.R
406          if R==F0:
407              UX=F0
408              UY=F0
409              UZ=F0
410              UXX=F0
411              UYX=F0
412              UZX=F0
413              UXY=F0
414              UYY=F0
415              UZY=F0
416              UXZ=F0
417              UYZ=F0
418              UZZ=F0
419              IRET=1
420              return
421          else:
422              PP1=POT1
423              PP2=POT2
424              PP3=POT3
425              PP4=POT4
426
427              self.UA0(XX,YY,DD,PP1,PP2,PP3,PP4)
428              DUA=self.DUA
429
430              for I in range(0,12):
431                if I<10: U[I]=U[I]-DUA[I]
432                if I>=10: U[I]=U[I]+DUA[I]
433#=======================================
434#=====  IMAGE-SOURCE CONTRIBUTION  =====
435#=======================================
436              DD=DEPTH-Z
437              self.DCCON1(XX,YY,DD)
438              self.UA0(XX,YY,DD,PP1,PP2,PP3,PP4)
439              DUA=self.DUA
440              self.UB0(XX,YY,DD,ZZ,PP1,PP2,PP3,PP4)
441              DUB=self.DUB
442              self.UC0(XX,YY,DD,ZZ,PP1,PP2,PP3,PP4)
443              DUC=self.DUC
444
445#-----
446              for I in range(0,12):
447                DU=DUA[I]+DUB[I]+ZZ*DUC[I]
448                if I>=9: DU=DU+DUC[I-9]
449                U[I]=U[I]+DU
450#=====
451              self.UX=U[0]
452              self.UY=U[1]
453              self.UZ=U[2]
454              self.UXX=U[3]
455              self.UYX=U[4]
456              self.UZX=U[5]
457              self.UXY=U[6]
458              self.UYY=U[7]
459              self.UZY=U[8]
460              self.UXZ=U[9]
461              self.UYZ=U[10]
462              self.UZZ=U[11]
463              self.IRET=0
464
465    def UA0(self,X,Y,D,POT1,POT2,POT3,POT4):
466
467#SUBROUTINE  UA0(X,Y,D,POT1,POT2,POT3,POT4,U)
468 #     IMPLICIT REAL*8 (A-H,O-Z)
469 #     DIMENSION U(12),DU(12)
470
471       """********************************************************************
472       C*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-A)             *****
473       C*****    DUE TO BURIED POINT SOURCE IN A SEMIINFINITE MEDIUM   *****
474       C********************************************************************
475       C***** INPUT
476       C*****   X,Y,D : STATION COORDINATES IN FAULT SYSTEM
477       C*****   POT1-POT4 : STRIKE-, DIP-, TENSILE- AND INFLATE-POTENCY
478       C***** OUTPUT
479       C*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES    """
480#
481# COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
482#      COMMON /C1/P,Q,S,T,XY,X2,Y2,D2,R,R2,R3,R5,QR,QRX,A3,A5,B3,C3,
483#    *           UY,VY,WY,UZ,VZ,WZ
484
485       from Numeric import zeros, Float
486
487
488       F0 = 0.0
489       F1 = 1.0
490       F3 = 3.0
491       PI2=6.283185307179586
492
493       ALP1=self.ALP1
494       ALP2=self.ALP2
495       ALP3=self.ALP3
496       ALP4=self.ALP4
497       ALP5=self.ALP4
498       SD=self.SD
499       CD=self.CD
500       SDSD=self.SDSD
501       CDCD=self.CDCD
502       SDCD=self.SDCD
503       S2D=self.S2D
504       C2D=self.C2D
505       P=self.P
506       Q=self.Q
507       S=self.S
508       T=self.T
509       XY=self.XY
510       X2=self.X2
511       Y2=self.Y2
512       D2=self.D2
513       R=self.R
514       R2=self.R2
515       R3=self.R3
516       R5=self.R5
517       QR=self.QR
518       QRX=self.QRX
519       A3=self.A3
520       A5=self.A5
521       B3=self.B3
522       C3=self.C3
523       UY=self.UY
524       VY=self.VY
525       WY=self.WY
526       UZ=self.UZ
527       VZ=self.VZ
528       WZ=self.WZ
529
530       DUA=zeros((12,1),Float)
531       DU=zeros((12,1),Float)
532       U=zeros((12,1),Float)
533#-----
534       for I in range(0,12):
535           U[I]=F0
536#======================================
537#=====  STRIKE-SLIP CONTRIBUTION  =====
538#======================================
539       if POT1 != F0:
540           DU[0]= ALP1*Q/R3    +ALP2*X2*QR
541           DU[1]= ALP1*X/R3*SD +ALP2*XY*QR
542           DU[2]=-ALP1*X/R3*CD +ALP2*X*D*QR
543           DU[3]= X*QR*(-ALP1 +ALP2*(F1+A5) )
544           DU[4]= ALP1*A3/R3*SD +ALP2*Y*QR*A5
545           DU[5]=-ALP1*A3/R3*CD +ALP2*D*QR*A5
546           DU[6]= ALP1*(SD/R3-Y*QR) +ALP2*F3*X2/R5*UY
547           DU[7]= F3*X/R5*(-ALP1*Y*SD +ALP2*(Y*UY+Q) )
548           DU[8]= F3*X/R5*( ALP1*Y*CD +ALP2*D*UY )
549           DU[9]= ALP1*(CD/R3+D*QR) +ALP2*F3*X2/R5*UZ
550           DU[10]= F3*X/R5*( ALP1*D*SD +ALP2*Y*UZ )
551           DU[11]= F3*X/R5*(-ALP1*D*CD +ALP2*(D*UZ-Q) )
552           for I in range(0,12):
553               U[I]=U[I]+POT1/PI2*DU[I]
554
555#===================================
556#=====  DIP-SLIP CONTRIBUTION  =====
557#===================================
558       if POT2 != F0:
559           DU[0]=            ALP2*X*P*QR
560           DU[1]= ALP1*S/R3 +ALP2*Y*P*QR
561           DU[2]=-ALP1*T/R3 +ALP2*D*P*QR
562           DU[3]=                 ALP2*P*QR*A5
563           DU[4]=-ALP1*F3*X*S/R5 -ALP2*Y*P*QRX
564           DU[5]= ALP1*F3*X*T/R5 -ALP2*D*P*QRX
565           DU[6]=                          ALP2*F3*X/R5*VY
566           DU[7]= ALP1*(S2D/R3-F3*Y*S/R5) +ALP2*(F3*Y/R5*VY+P*QR)
567           DU[8]=-ALP1*(C2D/R3-F3*Y*T/R5) +ALP2*F3*D/R5*VY
568           DU[9]=                          ALP2*F3*X/R5*VZ
569           DU[10]= ALP1*(C2D/R3+F3*D*S/R5) +ALP2*F3*Y/R5*VZ
570           DU[11]= ALP1*(S2D/R3-F3*D*T/R5) +ALP2*(F3*D/R5*VZ-P*QR)
571           for I in range(0,12):
572               U[I]=U[I]+POT2/PI2*DU[I]
573
574#========================================
575#=====  TENSILE-FAULT CONTRIBUTION  =====
576#========================================
577       if POT3!=F0:
578           DU[0]= ALP1*X/R3 -ALP2*X*Q*QR
579           DU[1]= ALP1*T/R3 -ALP2*Y*Q*QR
580           DU[2]= ALP1*S/R3 -ALP2*D*Q*QR
581           DU[3]= ALP1*A3/R3     -ALP2*Q*QR*A5
582           DU[4]=-ALP1*F3*X*T/R5 +ALP2*Y*Q*QRX
583           DU[5]=-ALP1*F3*X*S/R5 +ALP2*D*Q*QRX
584           DU[6]=-ALP1*F3*XY/R5           -ALP2*X*QR*WY
585           DU[7]= ALP1*(C2D/R3-F3*Y*T/R5) -ALP2*(Y*WY+Q)*QR
586           DU[8]= ALP1*(S2D/R3-F3*Y*S/R5) -ALP2*D*QR*WY
587           DU[9]= ALP1*F3*X*D/R5          -ALP2*X*QR*WZ
588           DU[10]=-ALP1*(S2D/R3-F3*D*T/R5) -ALP2*Y*QR*WZ
589           DU[11]= ALP1*(C2D/R3+F3*D*S/R5) -ALP2*(D*WZ-Q)*QR
590           for I in range(0,12):
591              U[I]=U[I]+POT3/PI2*DU[I]
592
593#=========================================
594#=====  INFLATE SOURCE CONTRIBUTION  =====
595#=========================================
596       if POT4 != F0:
597           DU[0]=-ALP1*X/R3
598           DU[1]=-ALP1*Y/R3
599           DU[2]=-ALP1*D/R3
600           DU[3]=-ALP1*A3/R3
601           DU[4]= ALP1*F3*XY/R5
602           DU[5]= ALP1*F3*X*D/R5
603           DU[6]= DU[4]
604           DU[7]=-ALP1*B3/R3
605           DU[8]= ALP1*F3*Y*D/R5
606           DU[9]=-DU[5]
607           DU[10]=-DU[8]
608           DU[11]= ALP1*C3/R3
609           for I in range(0,12):
610               U[I]=U[I]+POT4/PI2*DU[I]
611
612       #for I in range(0,12):
613           #DUA[I]=U[I]
614
615       self.DUA=U
616
617    def UB0(self,X,Y,D,Z,POT1,POT2,POT3,POT4):
618 #     SUBROUTINE  UB0(X,Y,D,Z,POT1,POT2,POT3,POT4,U)
619  #    IMPLICIT REAL*8 (A-H,O-Z)
620   #   DIMENSION U(12),DU(12)
621#
622        """********************************************************************
623        C*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-B)             *****
624        C*****    DUE TO BURIED POINT SOURCE IN A SEMIINFINITE MEDIUM   *****
625        C********************************************************************
626        C
627        C***** INPUT
628        C*****   X,Y,D,Z : STATION COORDINATES IN FAULT SYSTEM
629        C*****   POT1-POT4 : STRIKE-, DIP-, TENSILE- AND INFLATE-POTENCY
630        C***** OUTPUT
631        C*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES    """
632#
633 #     COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
634 #     COMMON /C1/P,Q,S,T,XY,X2,Y2,D2,R,R2,R3,R5,QR,QRX,A3,A5,B3,C3,
635 #    *           UY,VY,WY,UZ,VZ,WZ
636# F0,F1,F2,F3,F4,F5,F8,F9
637#     *        /0.D0,1.D0,2.D0,3.D0,4.D0,5.D0,8.D0,9.D0/
638#      DATA PI2/6.283185307179586D0/
639        from Numeric import zeros, Float
640
641        DUB=zeros((12,1),Float)
642        DU=zeros((12,1),Float)
643        U=zeros((12,1),Float)
644
645        F0=0.0
646        F1=1.0
647        F2=2.0
648        F3=3.0
649        F4=4.0
650        F5=5.0
651        F8=8.0
652        F9=9.0
653        PI2=6.283185307179586
654
655        ALP1=self.ALP1
656        ALP2=self.ALP2
657        ALP3=self.ALP3
658        ALP4=self.ALP4
659        ALP5=self.ALP4
660        SD=self.SD
661        CD=self.CD
662        SDSD=self.SDSD
663        CDCD=self.CDCD
664        SDCD=self.SDCD
665        S2D=self.S2D
666        C2D=self.C2D
667        P=self.P
668        Q=self.Q
669        S=self.S
670        T=self.T
671        XY=self.XY
672        X2=self.X2
673        Y2=self.Y2
674        D2=self.D2
675        R=self.R
676        R2=self.R2
677        R3=self.R3
678        R5=self.R5
679        QR=self.QR
680        QRX=self.QRX
681        A3=self.A3
682        A5=self.A5
683        B3=self.B3
684        C3=self.C3
685        UY=self.UY
686        VY=self.VY
687        WY=self.WY
688        UZ=self.UZ
689        VZ=self.VZ
690        WZ=self.WZ
691
692#-----
693        C=D+Z
694        RD=R+D
695        D12=F1/(R*RD*RD)
696        D32=D12*(F2*R+D)/R2
697        D33=D12*(F3*R+D)/(R2*RD)
698        D53=D12*(F8*R2+F9*R*D+F3*D2)/(R2*R2*RD)
699        D54=D12*(F5*R2+F4*R*D+D2)/R3*D12
700#-----
701        FI1= Y*(D12-X2*D33)
702        FI2= X*(D12-Y2*D33)
703        FI3= X/R3-FI2
704        FI4=-XY*D32
705        FI5= F1/(R*RD)-X2*D32
706        FJ1=-F3*XY*(D33-X2*D54)
707        FJ2= F1/R3-F3*D12+F3*X2*Y2*D54
708        FJ3= A3/R3-FJ2
709        FJ4=-F3*XY/R5-FJ1
710        FK1=-Y*(D32-X2*D53)
711        FK2=-X*(D32-Y2*D53)
712        FK3=-F3*X*D/R5-FK2
713#-----
714        for I in range(0,12):
715            U[I]=F0
716#======================================
717#=====  STRIKE-SLIP CONTRIBUTION  =====
718#======================================
719        if POT1!=F0:
720            DU[0]=-X2*QR  -ALP3*FI1*SD
721            DU[1]=-XY*QR  -ALP3*FI2*SD
722            DU[2]=-C*X*QR -ALP3*FI4*SD
723            DU[3]=-X*QR*(F1+A5) -ALP3*FJ1*SD
724            DU[4]=-Y*QR*A5      -ALP3*FJ2*SD
725            DU[5]=-C*QR*A5      -ALP3*FK1*SD
726            DU[6]=-F3*X2/R5*UY      -ALP3*FJ2*SD
727            DU[7]=-F3*XY/R5*UY-X*QR -ALP3*FJ4*SD
728            DU[8]=-F3*C*X/R5*UY     -ALP3*FK2*SD
729            DU[9]=-F3*X2/R5*UZ  +ALP3*FK1*SD
730            DU[10]=-F3*XY/R5*UZ  +ALP3*FK2*SD
731            DU[11]= F3*X/R5*(-C*UZ +ALP3*Y*SD)
732            for I in range(0,12):
733                U[I]=U[I]+POT1/PI2*DU[I]
734#===================================
735#=====  DIP-SLIP CONTRIBUTION  =====
736#===================================
737        if POT2!=F0:
738            DU[0]=-X*P*QR +ALP3*FI3*SDCD
739            DU[1]=-Y*P*QR +ALP3*FI1*SDCD
740            DU[2]=-C*P*QR +ALP3*FI5*SDCD
741            DU[3]=-P*QR*A5 +ALP3*FJ3*SDCD
742            DU[4]= Y*P*QRX +ALP3*FJ1*SDCD
743            DU[5]= C*P*QRX +ALP3*FK3*SDCD
744            DU[6]=-F3*X/R5*VY      +ALP3*FJ1*SDCD
745            DU[7]=-F3*Y/R5*VY-P*QR +ALP3*FJ2*SDCD
746            DU[8]=-F3*C/R5*VY      +ALP3*FK1*SDCD
747            DU[9]=-F3*X/R5*VZ -ALP3*FK3*SDCD
748            DU[10]=-F3*Y/R5*VZ -ALP3*FK1*SDCD
749            DU[11]=-F3*C/R5*VZ +ALP3*A3/R3*SDCD
750            for I in range(0,12):
751                U[I]=U[I]+POT2/PI2*DU[I]
752#========================================
753#=====  TENSILE-FAULT CONTRIBUTION  =====
754#========================================
755        if POT3!=F0:
756            DU[0]= X*Q*QR -ALP3*FI3*SDSD
757            DU[1]= Y*Q*QR -ALP3*FI1*SDSD
758            DU[2]= C*Q*QR -ALP3*FI5*SDSD
759            DU[3]= Q*QR*A5 -ALP3*FJ3*SDSD
760            DU[4]=-Y*Q*QRX -ALP3*FJ1*SDSD
761            DU[5]=-C*Q*QRX -ALP3*FK3*SDSD
762            DU[6]= X*QR*WY     -ALP3*FJ1*SDSD
763            DU[7]= QR*(Y*WY+Q) -ALP3*FJ2*SDSD
764            DU[8]= C*QR*WY     -ALP3*FK1*SDSD
765            DU[9]= X*QR*WZ +ALP3*FK3*SDSD
766            DU[10]= Y*QR*WZ +ALP3*FK1*SDSD
767            DU[11]= C*QR*WZ -ALP3*A3/R3*SDSD
768            for I in range(0,12):
769                U[I]=U[I]+POT3/PI2*DU[I]
770
771#=========================================
772#=====  INFLATE SOURCE CONTRIBUTION  =====
773#=========================================
774        if POT4!=F0:
775            DU[0]= ALP3*X/R3
776            DU[1]= ALP3*Y/R3
777            DU[2]= ALP3*D/R3
778            DU[3]= ALP3*A3/R3
779            DU[4]=-ALP3*F3*XY/R5
780            DU[5]=-ALP3*F3*X*D/R5
781            DU[6]= DU[4]
782            DU[7]= ALP3*B3/R3
783            DU[8]=-ALP3*F3*Y*D/R5
784            DU[9]=-DU[5]
785            DU[10]=-DU[8]
786            DU[11]=-ALP3*C3/R3
787            for I in range(0,12):
788                U[I]=U[I]+POT4/PI2*DU[I]
789
790        #for I in range(0,12):
791            #DUB[I]=U[I]
792        self.DUB=U
793
794    def UC0(self,X,Y,D,Z,POT1,POT2,POT3,POT4):
795  #    SUBROUTINE  UC0(X,Y,D,Z,POT1,POT2,POT3,POT4,U)
796  #    IMPLICIT REAL*8 (A-H,O-Z)
797  #    DIMENSION U(12),DU(12)
798
799      """********************************************************************
800      C*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-B)             *****
801      C*****    DUE TO BURIED POINT SOURCE IN A SEMIINFINITE MEDIUM   *****
802      C********************************************************************
803      C***** INPUT
804      C*****   X,Y,D,Z : STATION COORDINATES IN FAULT SYSTEM
805      C*****   POT1-POT4 : STRIKE-, DIP-, TENSILE- AND INFLATE-POTENCY
806      C***** OUTPUT
807      C*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES"""
808
809#      COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
810 #     COMMON /C1/P,Q,S,T,XY,X2,Y2,D2,R,R2,R3,R5,QR,QRX,A3,A5,B3,C3,um(6)
811#      DATA F0,F1,F2,F3,F5,F7,F10,F15
812#     *        /0.D0,1.D0,2.D0,3.D0,5.D0,7.D0,10.D0,15.D0/
813#      DATA PI2/6.283185307179586D0/
814
815
816      from Numeric import zeros, Float
817
818      DUC=zeros((12,1),Float)
819      DU=zeros((12,1),Float)
820      U=zeros((12,1),Float)
821
822      F0=0.0
823      F1=1.0
824      F2=2.0
825      F3=3.0
826      F5=5.0
827      F7=7.0
828      F10=10.0
829      F15=15.0
830      PI2=6.283185307179586
831
832      ALP1=self.ALP1
833      ALP2=self.ALP2
834      ALP3=self.ALP3
835      ALP4=self.ALP4
836      ALP5=self.ALP4
837      SD=self.SD
838      CD=self.CD
839      SDSD=self.SDSD
840      CDCD=self.CDCD
841      SDCD=self.SDCD
842      S2D=self.S2D
843      C2D=self.C2D
844      P=self.P
845      Q=self.Q
846      S=self.S
847      T=self.T
848      XY=self.XY
849      X2=self.X2
850      Y2=self.Y2
851      D2=self.D2
852      R=self.R
853      R2=self.R2
854      R3=self.R3
855      R5=self.R5
856      QR=self.QR
857      QRX=self.QRX
858      A3=self.A3
859      A5=self.A5
860      B3=self.B3
861      C3=self.C3
862      UY=self.UY
863      VY=self.VY
864      WY=self.WY
865      UZ=self.UZ
866      VZ=self.VZ
867      WZ=self.WZ
868
869#-----
870      C=D+Z
871      Q2=Q*Q
872      R7=R5*R2
873      A7=F1-F7*X2/R2
874      B5=F1-F5*Y2/R2
875      B7=F1-F7*Y2/R2
876      C5=F1-F5*D2/R2
877      C7=F1-F7*D2/R2
878      D7=F2-F7*Q2/R2
879      QR5=F5*Q/R2
880      QR7=F7*Q/R2
881      DR5=F5*D/R2
882#-----
883      for I in range(0,12):
884          U[I]=F0
885#======================================
886#=====  STRIKE-SLIP CONTRIBUTION  =====
887#======================================
888      if  POT1!=F0:
889        DU[0]=-ALP4*A3/R3*CD  +ALP5*C*QR*A5
890        DU[1]= F3*X/R5*( ALP4*Y*CD +ALP5*C*(SD-Y*QR5) )
891        DU[2]= F3*X/R5*(-ALP4*Y*SD +ALP5*C*(CD+D*QR5) )
892        DU[3]= ALP4*F3*X/R5*(F2+A5)*CD   -ALP5*C*QRX*(F2+A7)
893        DU[4]= F3/R5*( ALP4*Y*A5*CD +ALP5*C*(A5*SD-Y*QR5*A7) )
894        DU[5]= F3/R5*(-ALP4*Y*A5*SD +ALP5*C*(A5*CD+D*QR5*A7) )
895        DU[6]= DU[4]
896        DU[7]= F3*X/R5*( ALP4*B5*CD -ALP5*F5*C/R2*(F2*Y*SD+Q*B7) )
897        DU[8]= F3*X/R5*(-ALP4*B5*SD +ALP5*F5*C/R2*(D*B7*SD-Y*C7*CD) )
898        DU[9]= F3/R5*   (-ALP4*D*A5*CD +ALP5*C*(A5*CD+D*QR5*A7) )
899        DU[10]= F15*X/R7*( ALP4*Y*D*CD  +ALP5*C*(D*B7*SD-Y*C7*CD) )
900        DU[11]= F15*X/R7*(-ALP4*Y*D*SD  +ALP5*C*(F2*D*CD-Q*C7) )
901        for I in range(0,12):
902            U[I]=U[I]+POT1/PI2*DU[I]
903#===================================
904#=====  DIP-SLIP CONTRIBUTION  =====
905#===================================
906      if POT2!=F0:
907        DU[0]= ALP4*F3*X*T/R5          -ALP5*C*P*QRX
908        DU[1]=-ALP4/R3*(C2D-F3*Y*T/R2) +ALP5*F3*C/R5*(S-Y*P*QR5)
909        DU[2]=-ALP4*A3/R3*SDCD         +ALP5*F3*C/R5*(T+D*P*QR5)
910        DU[3]= ALP4*F3*T/R5*A5              -ALP5*F5*C*P*QR/R2*A7
911        DU[4]= F3*X/R5*(ALP4*(C2D-F5*Y*T/R2)-ALP5*F5*C/R2*(S-Y*P*QR7))
912        DU[5]= F3*X/R5*(ALP4*(F2+A5)*SDCD   -ALP5*F5*C/R2*(T+D*P*QR7))
913        DU[6]= DU[4]
914        DU[7]= (F3/R5*(ALP4*(F2*Y*C2D+T*B5)
915                               +ALP5*C*(S2D-F10*Y*S/R2-P*QR5*B7)))
916        DU[8]= F3/R5*(ALP4*Y*A5*SDCD-ALP5*C*((F3+A5)*C2D+Y*P*DR5*QR7))
917        DU[9]= F3*X/R5*(-ALP4*(S2D-T*DR5) -ALP5*F5*C/R2*(T+D*P*QR7))
918        DU[10]= (F3/R5*(-ALP4*(D*B5*C2D+Y*C5*S2D)
919                                -ALP5*C*((F3+A5)*C2D+Y*P*DR5*QR7)))
920        DU[11]= F3/R5*(-ALP4*D*A5*SDCD-ALP5*C*(S2D-F10*D*T/R2+P*QR5*C7))
921        for I in range(0,12):
922            U[I]=U[I]+POT2/PI2*DU[I]
923#========================================
924#=====  TENSILE-FAULT CONTRIBUTION  =====
925#========================================
926      if POT3!=F0:
927        DU[0]= F3*X/R5*(-ALP4*S +ALP5*(C*Q*QR5-Z))
928        DU[1]= ALP4/R3*(S2D-F3*Y*S/R2)+ALP5*F3/R5*(C*(T-Y+Y*Q*QR5)-Y*Z)
929        DU[2]=-ALP4/R3*(F1-A3*SDSD)   -ALP5*F3/R5*(C*(S-D+D*Q*QR5)-D*Z)
930        DU[3]=-ALP4*F3*S/R5*A5 +ALP5*(C*QR*QR5*A7-F3*Z/R5*A5)
931        DU[4]= (F3*X/R5*(-ALP4*(S2D-F5*Y*S/R2)
932                               -ALP5*F5/R2*(C*(T-Y+Y*Q*QR7)-Y*Z)))
933        DU[5]= (F3*X/R5*( ALP4*(F1-(F2+A5)*SDSD)
934                              +ALP5*F5/R2*(C*(S-D+D*Q*QR7)-D*Z)))
935        DU[6]= DU[4]
936        DU[7]= (F3/R5*(-ALP4*(F2*Y*S2D+S*B5)
937                     -ALP5*(C*(F2*SDSD+F10*Y*(T-Y)/R2-Q*QR5*B7)+Z*B5)))
938        DU[8]= (F3/R5*( ALP4*Y*(F1-A5*SDSD)
939                     +ALP5*(C*(F3+A5)*S2D-Y*DR5*(C*D7+Z))))
940        DU[9]= (F3*X/R5*(-ALP4*(C2D+S*DR5)
941                    +ALP5*(F5*C/R2*(S-D+D*Q*QR7)-F1-Z*DR5)))
942        DU[10]= (F3/R5*( ALP4*(D*B5*S2D-Y*C5*C2D)
943                   +ALP5*(C*((F3+A5)*S2D-Y*DR5*D7)-Y*(F1+Z*DR5))))
944        DU[11]= (F3/R5*(-ALP4*D*(F1-A5*SDSD)
945                    -ALP5*(C*(C2D+F10*D*(S-D)/R2-Q*QR5*C7)+Z*(F1+C5))))
946        for I in range(0,12):
947            U[I]=U[I]+POT3/PI2*DU[I]
948#=========================================
949#=====  INFLATE SOURCE CONTRIBUTION  =====
950#=========================================
951      if POT4!=F0:
952        DU[0]= ALP4*F3*X*D/R5
953        DU[1]= ALP4*F3*Y*D/R5
954        DU[2]= ALP4*C3/R3
955        DU[3]= ALP4*F3*D/R5*A5
956        DU[4]=-ALP4*F15*XY*D/R7
957        DU[5]=-ALP4*F3*X/R5*C5
958        DU[6]= DU[4]
959        DU[7]= ALP4*F3*D/R5*B5
960        DU[8]=-ALP4*F3*Y/R5*C5
961        DU[9]= DU[5]
962        DU[10]= DU[8]
963        DU[11]= ALP4*F3*D/R5*(F2+C5)
964        for I in range(0,12):
965            U[I]=U[I]+POT4/PI2*DU[I]
966
967      #for I in range(0,12):
968          #DUC[I]=U[I]
969      self.DUC=U
970
971    def DC3D(self,ALPHA,X,Y,Z,DEPTH,DIP,AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3):
972
973#      SUBROUTINE  DC3D(ALPHA,X,Y,Z,DEPTH,DIP,
974#     *              AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3,
975 #    *              UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,IRET)
976#      IMPLICIT REAL*8 (A-H,O-Z)
977#      REAL*4   ALPHA,X,Y,Z,DEPTH,DIP,AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3,
978#     *         UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,EPS
979
980      """********************************************************************
981      C*****                                                          *****
982      C*****    DISPLACEMENT AND STRAIN AT DEPTH                      *****
983      C*****    DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM   *****
984      C*****              CODED BY  Y.OKADA ... SEP.1991              *****
985      C*****              REVISED ... NOV.1991, APR.1992, MAY.1993,   *****
986      C*****                          JUL.1993                        *****
987      C********************************************************************
988      C
989      C***** INPUT
990      C*****   ALPHA : MEDIUM CONSTANT  (LAMBDA+MYU)/(LAMBDA+2*MYU)
991      C*****   X,Y,Z : COORDINATE OF OBSERVING POINT
992      C*****   DEPTH : DEPTH OF REFERENCE POINT
993      C*****   DIP   : DIP-ANGLE (DEGREE)
994      C*****   AL1,AL2   : FAULT LENGTH RANGE
995      C*****   AW1,AW2   : FAULT WIDTH RANGE
996      C*****   DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS
997      C
998      C***** OUTPUT
999      C*****   UX, UY, UZ  : DISPLACEMENT ( UNIT=(UNIT OF DISL)
1000      C*****   UXX,UYX,UZX : X-DERIVATIVE ( UNIT=(UNIT OF DISL) /
1001      C*****   UXY,UYY,UZY : Y-DERIVATIVE        (UNIT OF X,Y,Z,DEPTH,AL,AW) )
1002      C*****   UXZ,UYZ,UZZ : Z-DERIVATIVE
1003      C*****   IRET        : RETURN CODE  ( =0....NORMAL,   =1....SINGULAR )  """
1004#
1005#      COMMON /C0/DUMMY(5),SD,CD,dumm(5)
1006 #     DIMENSION  XI(2),ET(2),KXI(2),KET(2)
1007 #     DIMENSION  U(12),DU(12),DUA(12),DUB(12),DUC(12)
1008      from Numeric import zeros, Float
1009      from math import sqrt
1010      F0 =0.0
1011      EPS=0.000001
1012
1013      XI=zeros(2,Float)
1014      ET=zeros(2,Float)
1015      KXI=zeros(2,Float)
1016      KET=zeros(2,Float)
1017      U=zeros(12,Float)
1018      DU=zeros(12,Float)
1019      DUA=zeros(12,Float)
1020      DUB=zeros(12,Float)
1021      DUC=zeros(12,Float)
1022
1023#-----
1024      if Z>0: print '** POSITIVE Z WAS GIVEN IN SUB-DC3D'
1025
1026      for I in range(0,12):
1027        U[I]=F0
1028        DUA[I]=F0
1029        DUB[I]=F0
1030        DUC[I]=F0
1031
1032      AALPHA=ALPHA
1033      DDIP=DIP
1034      self.DCC0N0(AALPHA,DDIP)
1035      CD=self.CD
1036      SD=self.SD
1037
1038#-----
1039      ZZ=Z
1040      DD1=DISL1
1041      DD2=DISL2
1042      DD3=DISL3
1043      XI[0]=X-AL1
1044      XI[1]=X-AL2
1045      #if X==-6 and Y==-1:
1046          #print AL1
1047          #print AL2
1048          #print 'hello'
1049
1050      if abs(XI[0])<EPS: XI[0]=F0
1051      if abs(XI[1])<EPS: XI[1]=F0
1052#======================================
1053#=====  REAL-SOURCE CONTRIBUTION  =====
1054#======================================
1055      D=DEPTH+Z
1056      P=Y*CD+D*SD
1057      Q=Y*SD-D*CD
1058      ET[0]=P-AW1
1059      ET[1]=P-AW2
1060      if abs(Q)<EPS:  Q=F0
1061      if abs(ET[0])<EPS: ET[0]=F0
1062      if abs(ET[1])<EPS: ET[1]=F0
1063
1064#--------------------------------
1065#----- REJECT SINGULAR CASE -----
1066#--------------------------------
1067#----- ON FAULT EDGE
1068      if (Q==F0 and ((XI[0]*XI[1]<F0 and ET[0]*ET[1]==F0)
1069                   or (ET[0]*ET[1]<F0 and XI[0]*XI[1]==F0) )):
1070          self.UX=F0
1071          self.UY=F0
1072          self.UZ=F0
1073          self.UXX=F0
1074          self.UYX=F0
1075          self.UZX=F0
1076          self.UXY=F0
1077          self.UYY=F0
1078          self.UZY=F0
1079          self.UXZ=F0
1080          self.UYZ=F0
1081          self.UZZ=F0
1082          self.IRET=1
1083          return
1084#----- ON NEGATIVE EXTENSION OF FAULT EDGE
1085      KXI[0]=0
1086      KXI[1]=0
1087      KET[0]=0
1088      KET[1]=0
1089      R12= sqrt(XI[0]*XI[0]+ET[1]*ET[1]+Q*Q)
1090      R21= sqrt(XI[1]*XI[1]+ET[0]*ET[0]+Q*Q)
1091      R22= sqrt(XI[1]*XI[1]+ET[1]*ET[1]+Q*Q)
1092      if XI[0]<F0 and R21+XI[1]<EPS: KXI[0]=1
1093      if XI[0]<F0 and R22+XI[1]<EPS: KXI[1]=1
1094      if ET[0]<F0 and R12+ET[1]<EPS: KET[0]=1
1095      if ET[0]<F0 and R22+ET[1]<EPS: KET[1]=1
1096#=====
1097      for K in range(0,2):
1098          for J in range(0,2):
1099             self.DCCON2(XI[J],ET[K],Q,SD,CD,KXI[K],KET[J])
1100             self.UA(XI[J],ET[K],Q,DD1,DD2,DD3)
1101             DUA=self.DUA
1102             for I in range(0,10,3):
1103                  DU[I]  =-DUA[I]
1104                  DU[I+1]=-DUA[I+1]*CD+DUA[I+2]*SD
1105                  DU[I+2]=-DUA[I+1]*SD-DUA[I+2]*CD
1106                  if I==9:
1107                      DU[I]  =-DU[I]
1108                      DU[I+1]=-DU[I+1]
1109                      DU[I+2]=-DU[I+2]
1110                  else: continue
1111             for I in range(0,12):
1112                  if (J+K)!=1: U[I]=U[I]+DU[I]
1113                  if (J+K)==1: U[I]=U[I]-DU[I]
1114
1115#=======================================
1116#=====  IMAGE-SOURCE CONTRIBUTION  =====
1117#=======================================
1118      D=DEPTH-Z
1119      P=Y*CD+D*SD
1120      Q=Y*SD-D*CD
1121      ET[0]=P-AW1
1122      ET[1]=P-AW2
1123      if abs(Q)<EPS:  Q=F0
1124      if abs(ET[0])<EPS: ET[0]=F0
1125      if abs(ET[1])<EPS: ET[1]=F0
1126#--------------------------------
1127#----- REJECT SINGULAR CASE -----
1128#--------------------------------
1129#----- ON FAULT EDGE
1130      if (Q==F0 and ((XI[0]*XI[1]<F0 and ET[0]*ET[1]==F0)
1131                   or (ET[0]*ET[1]<F0 and XI[0]*XI[1]==F0) )):
1132          self.UX=F0
1133          self.UY=F0
1134          self.UZ=F0
1135          self.UXX=F0
1136          self.UYX=F0
1137          self.UZX=F0
1138          self.UXY=F0
1139          self.UYY=F0
1140          self.UZY=F0
1141          self.UXZ=F0
1142          self.UYZ=F0
1143          self.UZZ=F0
1144          self.IRET=1
1145          return
1146#----- ON NEGATIVE EXTENSION OF FAULT EDGE
1147      KXI[0]=0
1148      KXI[1]=0
1149      KET[0]=0
1150      KET[1]=0
1151      R12=sqrt(XI[0]*XI[0]+ET[1]*ET[1]+Q*Q)
1152      R21=sqrt(XI[1]*XI[1]+ET[0]*ET[0]+Q*Q)
1153      R22=sqrt(XI[1]*XI[1]+ET[1]*ET[1]+Q*Q)
1154      if XI[0]<F0 and R21+XI[1]<EPS: KXI[0]=1
1155      if XI[0]<F0 and R22+XI[1]<EPS: KXI[1]=1
1156      if ET[0]<F0 and R12+ET[1]<EPS: KET[0]=1
1157      if ET[0]<F0 and R22+ET[1]<EPS: KET[1]=1
1158#=====
1159      for K in range(0,2):
1160          for J in range(0,2):
1161              self.DCCON2(XI[J],ET[K],Q,SD,CD,KXI[K],KET[J])
1162              self.UA(XI[J],ET[K],Q,DD1,DD2,DD3)
1163              DUA=self.DUA
1164              self.UB(XI[J],ET[K],Q,DD1,DD2,DD3)
1165              DUB=self.DUB
1166              self.UC(XI[J],ET[K],Q,ZZ,DD1,DD2,DD3)
1167              DUC=self.DUC
1168#-----
1169              for I in range(0,10,3):
1170                  DU[I]=DUA[I]+DUB[I]+Z*DUC[I]
1171                  DU[I+1]=((DUA[I+1]+DUB[I+1]+Z*DUC[I+1])*CD
1172                        -(DUA[I+2]+DUB[I+2]+Z*DUC[I+2])*SD )
1173                  DU[I+2]=((DUA[I+1]+DUB[I+1]-Z*DUC[I+1])*SD
1174                        +(DUA[I+2]+DUB[I+2]-Z*DUC[I+2])*CD)
1175                  if I==9:
1176                      DU[9]=DU[9]+DUC[0]
1177                      DU[10]=DU[10]+DUC[1]*CD-DUC[2]*SD
1178                      DU[11]=DU[11]-DUC[1]*SD-DUC[2]*CD
1179              for I in range(0,12):
1180                  if (J+K)!=1: U[I]=U[I]+DU[I]
1181                  if (J+K)==1: U[I]=U[I]-DU[I]
1182#=====
1183      self.UX=U[0]
1184      self.UY=U[1]
1185      self.UZ=U[2]
1186      self.UXX=U[3]
1187      self.UYX=U[4]
1188      self.UZX=U[5]
1189      self.UXY=U[6]
1190      self.UYY=U[7]
1191      self.UZY=U[8]
1192      self.UXZ=U[9]
1193      self.UYZ=U[10]
1194      self.UZZ=U[11]
1195      self.IRET=0
1196
1197
1198    def UA(self,XI,ET,Q,DISL1,DISL2,DISL3):
1199#      IMPLICIT REAL*8 (A-H,O-Z)
1200#      DIMENSION U(12),DU(12)
1201#
1202#********************************************************************
1203#*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-A)             *****
1204#*****    DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM   *****
1205#********************************************************************
1206#
1207#***** INPUT
1208#*****   XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM
1209#*****   DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS
1210#***** OUTPUT
1211#*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES
1212#
1213 #     COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
1214 #     COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32,
1215#     *           EY,EZ,FY,FZ,GY,GZ,HY,HZ
1216#      DATA F0,F2,PI2/0.D0,2.D0,6.283185307179586D0/
1217
1218
1219      from Numeric import zeros, Float
1220
1221      U=zeros(12,Float)
1222      DU=zeros(12,Float)
1223      DUA=zeros(12,Float)
1224      F0 =0.0
1225      F2=2.0
1226      PI2=6.283185307179586
1227
1228      ALP1=self.ALP1
1229      ALP2=self.ALP2
1230      ALP3=self.ALP3
1231      ALP4=self.ALP4
1232      ALP5=self.ALP4
1233      SD=self.SD
1234      CD=self.CD
1235      SDSD=self.SDSD
1236      CDCD=self.CDCD
1237      SDCD=self.SDCD
1238      S2D=self.S2D
1239      C2D=self.C2D
1240      XI2=self.XI2
1241      ET2=self.ET2
1242      Q2=self.Q2
1243      R=self.R
1244      R2=self.R2
1245      R3=self.R3
1246      R5=self.R5
1247      Y=self.Y
1248      D=self.D
1249      TT=self.TT
1250      ALX=self.ALX
1251      ALE=self.ALE
1252      X11=self.X11
1253      Y11=self.Y11
1254      X32=self.X32
1255      Y32=self.Y32
1256      EY=self.EY
1257      EZ=self.EZ
1258      FY=self.FY
1259      FZ=self.FZ
1260      GY=self.GY
1261      GZ=self.GZ
1262      HY=self.HY
1263      HZ=self.HZ
1264
1265#-----
1266      for I in range(0,12):
1267        U[I]=F0
1268      XY=XI*Y11
1269      QX=Q *X11
1270      QY=Q *Y11
1271#======================================
1272#=====  STRIKE-SLIP CONTRIBUTION  =====
1273#======================================
1274      if DISL1 != F0 :
1275        DU[0]=    TT/F2 +ALP2*XI*QY
1276        DU[1]=           ALP2*Q/R
1277        DU[2]= ALP1*ALE -ALP2*Q*QY
1278        DU[3]=-ALP1*QY  -ALP2*XI2*Q*Y32
1279        DU[4]=          -ALP2*XI*Q/R3
1280        DU[5]= ALP1*XY  +ALP2*XI*Q2*Y32
1281        DU[6]= ALP1*XY*SD        +ALP2*XI*FY+D/F2*X11
1282        DU[7]=                    ALP2*EY
1283        DU[8]= ALP1*(CD/R+QY*SD) -ALP2*Q*FY
1284        DU[9]= ALP1*XY*CD        +ALP2*XI*FZ+Y/F2*X11
1285        DU[10]=                    ALP2*EZ
1286        DU[11]=-ALP1*(SD/R-QY*CD) -ALP2*Q*FZ
1287        for I in range(0,12):
1288            U[I]=U[I]+DISL1/PI2*DU[I]
1289#======================================
1290#=====    DIP-SLIP CONTRIBUTION   =====
1291#======================================
1292      if DISL2!=F0:
1293        DU[0]=           ALP2*Q/R
1294        DU[1]=    TT/F2 +ALP2*ET*QX
1295        DU[2]= ALP1*ALX -ALP2*Q*QX
1296        DU[3]=        -ALP2*XI*Q/R3
1297        DU[4]= -QY/F2 -ALP2*ET*Q/R3
1298        DU[5]= ALP1/R +ALP2*Q2/R3
1299        DU[6]=                      ALP2*EY
1300        DU[7]= ALP1*D*X11+XY/F2*SD +ALP2*ET*GY
1301        DU[8]= ALP1*Y*X11          -ALP2*Q*GY
1302        DU[9]=                      ALP2*EZ
1303        DU[10]= ALP1*Y*X11+XY/F2*CD +ALP2*ET*GZ
1304        DU[11]=-ALP1*D*X11          -ALP2*Q*GZ
1305        for I in range(0,12):
1306            U[I]=U[I]+DISL2/PI2*DU[I]
1307#========================================
1308#=====  TENSILE-FAULT CONTRIBUTION  =====
1309#========================================
1310      if DISL3!=F0:
1311        DU[0]=-ALP1*ALE -ALP2*Q*QY
1312        DU[1]=-ALP1*ALX -ALP2*Q*QX
1313        DU[2]=    TT/F2 -ALP2*(ET*QX+XI*QY)
1314        DU[3]=-ALP1*XY  +ALP2*XI*Q2*Y32
1315        DU[4]=-ALP1/R   +ALP2*Q2/R3
1316        DU[5]=-ALP1*QY  -ALP2*Q*Q2*Y32
1317        DU[6]=-ALP1*(CD/R+QY*SD)  -ALP2*Q*FY
1318        DU[7]=-ALP1*Y*X11         -ALP2*Q*GY
1319        DU[8]= ALP1*(D*X11+XY*SD) +ALP2*Q*HY
1320        DU[9]= ALP1*(SD/R-QY*CD)  -ALP2*Q*FZ
1321        DU[10]= ALP1*D*X11         -ALP2*Q*GZ
1322        DU[11]= ALP1*(Y*X11+XY*CD) +ALP2*Q*HZ
1323        for  I in range(0,12):
1324            U[I]=U[I]+DISL3/PI2*DU[I]
1325      #for I in range (0,12):
1326          #DUA[I]=U[I]
1327      self.DUA=U
1328
1329
1330    def UB(self,XI,ET,Q,DISL1,DISL2,DISL3):
1331      from math import sqrt, atan,log
1332 #SUBROUTINE  UB(XI,ET,Q,DISL1,DISL2,DISL3,U)
1333 #     IMPLICIT REAL*8 (A-H,O-Z)
1334 #     DIMENSION U(12),DU(12)
1335
1336      """********************************************************************
1337      C*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-B)             *****
1338      C*****    DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM   *****
1339      C********************************************************************
1340      C
1341      C***** INPUT
1342      C*****   XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM
1343      C*****   DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS
1344      C***** OUTPUT
1345      C*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES """
1346
1347#      COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
1348#      COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32,
1349#     *           EY,EZ,FY,FZ,GY,GZ,HY,HZ
1350#      DATA  F0,F1,F2,PI2/0.D0,1.D0,2.D0,6.283185307179586D0/
1351
1352      from Numeric import zeros, Float
1353
1354      DUB=zeros(12,Float)
1355      DU=zeros(12,Float)
1356      U=zeros(12,Float)
1357
1358      F0=0.0
1359      F1=1.0
1360      F2=2.0
1361      PI2=6.283185307179586
1362
1363      ALP1=self.ALP1
1364      ALP2=self.ALP2
1365      ALP3=self.ALP3
1366      ALP4=self.ALP4
1367      ALP5=self.ALP4
1368      SD=self.SD
1369      CD=self.CD
1370      SDSD=self.SDSD
1371      CDCD=self.CDCD
1372      SDCD=self.SDCD
1373      S2D=self.S2D
1374      C2D=self.C2D
1375      XI2=self.XI2
1376      ET2=self.ET2
1377      Q2=self.Q2
1378      R=self.R
1379      R2=self.R2
1380      R3=self.R3
1381      R5=self.R5
1382      Y=self.Y
1383      D=self.D
1384      TT=self.TT
1385      ALX=self.ALX
1386      ALE=self.ALE
1387      X11=self.X11
1388      Y11=self.Y11
1389      X32=self.X32
1390      Y32=self.Y32
1391      EY=self.EY
1392      EZ=self.EZ
1393      FY=self.FY
1394      FZ=self.FZ
1395      GY=self.GY
1396      GZ=self.GZ
1397      HY=self.HY
1398      HZ=self.HZ
1399
1400      RD=R+D
1401      D11=F1/(R*RD)
1402      AJ2=XI*Y/RD*D11
1403      AJ5=-(D+Y*Y/RD)*D11
1404      if CD!=F0:
1405        if XI==F0:
1406          AI4=F0
1407        else:
1408          X=sqrt(XI2+Q2)
1409          AI4=(F1/CDCD*( XI/RD*SDCD
1410            +F2*atan((ET*(X+Q*CD)+X*(R+X)*SD)/(XI*(R+X)*CD)) ))
1411        AI3=(Y*CD/RD-ALE+SD*log(RD))/CDCD
1412        AK1=XI*(D11-Y11*SD)/CD
1413        AK3=(Q*Y11-Y*D11)/CD
1414        AJ3=(AK1-AJ2*SD)/CD
1415        AJ6=(AK3-AJ5*SD)/CD
1416      else:
1417        RD2=RD*RD
1418        AI3=(ET/RD+Y*Q/RD2-ALE)/F2
1419        AI4=XI*Y/RD2/F2
1420        AK1=XI*Q/RD*D11
1421        AK3=SD/RD*(XI2*D11-F1)
1422        AJ3=-XI/RD2*(Q2*D11-F1/F2)
1423        AJ6=-Y/RD2*(XI2*D11-F1/F2)
1424#-----
1425      XY=XI*Y11
1426      AI1=-XI/RD*CD-AI4*SD
1427      AI2= log(RD)+AI3*SD
1428      AK2= F1/R+AK3*SD
1429      AK4= XY*CD-AK1*SD
1430      AJ1= AJ5*CD-AJ6*SD
1431      AJ4=-XY-AJ2*CD+AJ3*SD
1432#=====
1433      for I in range(0,12):
1434        U[I]=F0
1435        QX=Q*X11
1436        QY=Q*Y11
1437#======================================
1438#=====  STRIKE-SLIP CONTRIBUTION  =====
1439#======================================
1440      if DISL1!=F0:
1441        DU[0]=-XI*QY-TT -ALP3*AI1*SD
1442        DU[1]=-Q/R      +ALP3*Y/RD*SD
1443        DU[2]= Q*QY     -ALP3*AI2*SD
1444        DU[3]= XI2*Q*Y32 -ALP3*AJ1*SD
1445        DU[4]= XI*Q/R3   -ALP3*AJ2*SD
1446        DU[5]=-XI*Q2*Y32 -ALP3*AJ3*SD
1447        DU[6]=-XI*FY-D*X11 +ALP3*(XY+AJ4)*SD
1448        DU[7]=-EY          +ALP3*(F1/R+AJ5)*SD
1449        DU[8]= Q*FY        -ALP3*(QY-AJ6)*SD
1450        DU[9]=-XI*FZ-Y*X11 +ALP3*AK1*SD
1451        DU[10]=-EZ          +ALP3*Y*D11*SD
1452        DU[11]= Q*FZ        +ALP3*AK2*SD
1453        for I in range(0,12):
1454            U[I]=U[I]+DISL1/PI2*DU[I]
1455#======================================
1456#=====    DIP-SLIP CONTRIBUTION   =====
1457#======================================
1458      if DISL2!=F0:
1459        DU[0]=-Q/R      +ALP3*AI3*SDCD
1460        DU[1]=-ET*QX-TT -ALP3*XI/RD*SDCD
1461        DU[2]= Q*QX     +ALP3*AI4*SDCD
1462        DU[3]= XI*Q/R3     +ALP3*AJ4*SDCD
1463        DU[4]= ET*Q/R3+QY  +ALP3*AJ5*SDCD
1464        DU[5]=-Q2/R3       +ALP3*AJ6*SDCD
1465        DU[6]=-EY          +ALP3*AJ1*SDCD
1466        DU[7]=-ET*GY-XY*SD +ALP3*AJ2*SDCD
1467        DU[8]= Q*GY        +ALP3*AJ3*SDCD
1468        DU[9]=-EZ          -ALP3*AK3*SDCD
1469        DU[10]=-ET*GZ-XY*CD -ALP3*XI*D11*SDCD
1470        DU[11]= Q*GZ        -ALP3*AK4*SDCD
1471        for I in range(0,12):
1472            U[I]=U[I]+DISL2/PI2*DU[I]
1473#========================================
1474#=====  TENSILE-FAULT CONTRIBUTION  =====
1475#========================================
1476      if DISL3!=F0:
1477        DU[0]= Q*QY           -ALP3*AI3*SDSD
1478        DU[1]= Q*QX           +ALP3*XI/RD*SDSD
1479        DU[2]= ET*QX+XI*QY-TT -ALP3*AI4*SDSD
1480        DU[3]=-XI*Q2*Y32 -ALP3*AJ4*SDSD
1481        DU[4]=-Q2/R3     -ALP3*AJ5*SDSD
1482        DU[5]= Q*Q2*Y32  -ALP3*AJ6*SDSD
1483        DU[6]= Q*FY -ALP3*AJ1*SDSD
1484        DU[7]= Q*GY -ALP3*AJ2*SDSD
1485        DU[8]=-Q*HY -ALP3*AJ3*SDSD
1486        DU[9]= Q*FZ +ALP3*AK3*SDSD
1487        DU[10]= Q*GZ +ALP3*XI*D11*SDSD
1488        DU[11]=-Q*HZ +ALP3*AK4*SDSD
1489        for I in range(0,12):
1490            U[I]=U[I]+DISL3/PI2*DU[I]
1491
1492      #for I in range(0,12):
1493          #DUB[I]=U[I]
1494      self.DUB=U
1495
1496
1497
1498    def UC(self,XI,ET,Q,Z,DISL1,DISL2,DISL3):
1499 #     SUBROUTINE  UC(XI,ET,Q,Z,DISL1,DISL2,DISL3,U)
1500 #     IMPLICIT REAL*8 (A-H,O-Z)
1501  #    DIMENSION U(12),DU(12)
1502
1503
1504      """********************************************************************
1505      C*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-C)             *****
1506      C*****    DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM   *****
1507      C********************************************************************
1508      C
1509      C***** INPUT
1510      C*****   XI,ET,Q,Z   : STATION COORDINATES IN FAULT SYSTEM
1511      C*****   DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS
1512      C***** OUTPUT
1513      C*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES  """
1514#     COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
1515#      COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32,
1516 #    *           EY,EZ,FY,FZ,GY,GZ,HY,HZ
1517 #     DATA F0,F1,F2,F3,PI2/0.D0,1.D0,2.D0,3.D0,6.283185307179586D0/
1518
1519      from Numeric import zeros, Float
1520
1521      DUC=zeros(12,Float)
1522      DU=zeros(12,Float)
1523      U=zeros(12,Float)
1524
1525      F0=0.0
1526      F1=1.0
1527      F2=2.0
1528      F3=3.0
1529      PI2=6.283185307179586
1530
1531      ALP1=self.ALP1
1532      ALP2=self.ALP2
1533      ALP3=self.ALP3
1534      ALP4=self.ALP4
1535      ALP5=self.ALP4
1536      SD=self.SD
1537      CD=self.CD
1538      SDSD=self.SDSD
1539      CDCD=self.CDCD
1540      SDCD=self.SDCD
1541      S2D=self.S2D
1542      C2D=self.C2D
1543      XI2=self.XI2
1544      ET2=self.ET2
1545      Q2=self.Q2
1546      R=self.R
1547      R2=self.R2
1548      R3=self.R3
1549      R5=self.R5
1550      Y=self.Y
1551      D=self.D
1552      TT=self.TT
1553      ALX=self.ALX
1554      ALE=self.ALE
1555      X11=self.X11
1556      Y11=self.Y11
1557      X32=self.X32
1558      Y32=self.Y32
1559      EY=self.EY
1560      EZ=self.EZ
1561      FY=self.FY
1562      FZ=self.FZ
1563      GY=self.GY
1564      GZ=self.GZ
1565      HY=self.HY
1566      HZ=self.HZ
1567#-----
1568      C=D+Z
1569      X53=(8.0*R2+9.0*R*XI+F3*XI2)*X11*X11*X11/R2
1570      Y53=(8.0*R2+9.0*R*ET+F3*ET2)*Y11*Y11*Y11/R2
1571      H=Q*CD-Z
1572      Z32=SD/R3-H*Y32
1573      Z53=F3*SD/R5-H*Y53
1574      Y0=Y11-XI2*Y32
1575      Z0=Z32-XI2*Z53
1576      PPY=CD/R3+Q*Y32*SD
1577      PPZ=SD/R3-Q*Y32*CD
1578      QQ=Z*Y32+Z32+Z0
1579      QQY=F3*C*D/R5-QQ*SD
1580      QQZ=F3*C*Y/R5-QQ*CD+Q*Y32
1581      XY=XI*Y11
1582      QX=Q*X11
1583      QY=Q*Y11
1584      QR=F3*Q/R5
1585      CQX=C*Q*X53
1586      CDR=(C+D)/R3
1587      YY0=Y/R3-Y0*CD
1588#=====
1589      for  I in range(1,12):
1590          U[I]=F0
1591#======================================
1592#=====  STRIKE-SLIP CONTRIBUTION  =====
1593#======================================
1594      if DISL1!=F0:
1595        DU[0]= ALP4*XY*CD           -ALP5*XI*Q*Z32
1596        DU[1]= ALP4*(CD/R+F2*QY*SD) -ALP5*C*Q/R3
1597        DU[2]= ALP4*QY*CD           -ALP5*(C*ET/R3-Z*Y11+XI2*Z32)
1598        DU[3]= ALP4*Y0*CD                  -ALP5*Q*Z0
1599        DU[4]=-ALP4*XI*(CD/R3+F2*Q*Y32*SD) +ALP5*C*XI*QR
1600        DU[5]=-ALP4*XI*Q*Y32*CD            +ALP5*XI*(F3*C*ET/R5-QQ)
1601        DU[6]=-ALP4*XI*PPY*CD    -ALP5*XI*QQY
1602        DU[7]= (ALP4*F2*(D/R3-Y0*SD)*SD-Y/R3*CD
1603                                 -ALP5*(CDR*SD-ET/R3-C*Y*QR))
1604        DU[8]=-ALP4*Q/R3+YY0*SD  +ALP5*(CDR*CD+C*D*QR-(Y0*CD+Q*Z0)*SD)
1605        DU[9]= ALP4*XI*PPZ*CD    -ALP5*XI*QQZ
1606        DU[10]= ALP4*F2*(Y/R3-Y0*CD)*SD+D/R3*CD -ALP5*(CDR*CD+C*D*QR)
1607        DU[11]=         YY0*CD    -ALP5*(CDR*SD-C*Y*QR-Y0*SDSD+Q*Z0*CD)
1608        for I in range(0,12):
1609            U[I]=U[I]+DISL1/PI2*DU[I]
1610
1611#======================================
1612#=====    DIP-SLIP CONTRIBUTION   =====
1613#======================================
1614      if DISL2!=F0 :
1615        DU[0]= ALP4*CD/R -QY*SD -ALP5*C*Q/R3
1616        DU[1]= ALP4*Y*X11       -ALP5*C*ET*Q*X32
1617        DU[2]=     -D*X11-XY*SD -ALP5*C*(X11-Q2*X32)
1618        DU[3]=-ALP4*XI/R3*CD +ALP5*C*XI*QR +XI*Q*Y32*SD
1619        DU[4]=-ALP4*Y/R3     +ALP5*C*ET*QR
1620        DU[5]=    D/R3-Y0*SD +ALP5*C/R3*(F1-F3*Q2/R2)
1621        DU[6]=-ALP4*ET/R3+Y0*SDSD -ALP5*(CDR*SD-C*Y*QR)
1622        DU[7]= ALP4*(X11-Y*Y*X32) -ALP5*C*((D+F2*Q*CD)*X32-Y*ET*Q*X53)
1623        DU[8]=  XI*PPY*SD+Y*D*X32 +ALP5*C*((Y+F2*Q*SD)*X32-Y*Q2*X53)
1624        DU[9]=      -Q/R3+Y0*SDCD -ALP5*(CDR*CD+C*D*QR)
1625        DU[10]= ALP4*Y*D*X32       -ALP5*C*((Y-F2*Q*SD)*X32+D*ET*Q*X53)
1626        DU[11]=-XI*PPZ*SD+X11-D*D*X32-ALP5*C*((D-F2*Q*CD)*X32-D*Q2*X53)
1627        for I in range(0,12):
1628            U[I]=U[I]+DISL2/PI2*DU[I]
1629#========================================
1630#=====  TENSILE-FAULT CONTRIBUTION  =====
1631#========================================
1632      if DISL3!=F0:
1633        DU[0]=-ALP4*(SD/R+QY*CD)   -ALP5*(Z*Y11-Q2*Z32)
1634        DU[1]= ALP4*F2*XY*SD+D*X11 -ALP5*C*(X11-Q2*X32)
1635        DU[2]= ALP4*(Y*X11+XY*CD)  +ALP5*Q*(C*ET*X32+XI*Z32)
1636        DU[3]= ALP4*XI/R3*SD+XI*Q*Y32*CD+ALP5*XI*(F3*C*ET/R5-F2*Z32-Z0)
1637        DU[4]= ALP4*F2*Y0*SD-D/R3 +ALP5*C/R3*(F1-F3*Q2/R2)
1638        DU[5]=-ALP4*YY0           -ALP5*(C*ET*QR-Q*Z0)
1639        DU[6]= ALP4*(Q/R3+Y0*SDCD)   +ALP5*(Z/R3*CD+C*D*QR-Q*Z0*SD)
1640        DU[7]=(-ALP4*F2*XI*PPY*SD-Y*D*X32
1641                        +ALP5*C*((Y+F2*Q*SD)*X32-Y*Q2*X53))
1642        DU[8]=(-ALP4*(XI*PPY*CD-X11+Y*Y*X32)
1643                         +ALP5*(C*((D+F2*Q*CD)*X32-Y*ET*Q*X53)+XI*QQY))
1644        DU[9]=  -ET/R3+Y0*CDCD -ALP5*(Z/R3*SD-C*Y*QR-Y0*SDSD+Q*Z0*CD)
1645        DU[10]= (ALP4*F2*XI*PPZ*SD-X11+D*D*X32
1646                        -ALP5*C*((D-F2*Q*CD)*X32-D*Q2*X53))
1647        DU[11]= (ALP4*(XI*PPZ*CD+Y*D*X32)
1648                        +ALP5*(C*((Y-F2*Q*SD)*X32+D*ET*Q*X53)+XI*QQZ))
1649        for I in range(0,12):
1650            U[I]=U[I]+DISL3/PI2*DU[I]
1651
1652      #for I in range(0,12):
1653          #DUC[I]=U[I]
1654      self.DUC=U
1655
1656    def DCC0N0(self,ALPHA,DIP):
1657      from math import sin, cos
1658#      SUBROUTINE  DCCON0(ALPHA,DIP)
1659#      IMPLICIT REAL*8 (A-H,O-Z)
1660#
1661      """*******************************************************************
1662      C*****   CALCULATE MEDIUM CONSTANTS AND FAULT-DIP CONSTANTS    *****
1663      C*******************************************************************
1664      C
1665      C***** INPUT
1666      C*****   ALPHA : MEDIUM CONSTANT  (LAMBDA+MYU)/(LAMBDA+2*MYU)
1667      C*****   DIP   : DIP-ANGLE (DEGREE)
1668      C### CAUTION ### IF COS(DIP) IS SUFFICIENTLY SMALL, IT IS SET TO ZERO """
1669
1670#      COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D
1671 #     DATA F0,F1,F2,PI2/0.D0,1.D0,2.D0,6.283185307179586D0/
1672#      DATA EPS/1.D-6/
1673      F0=0.0
1674      F1=1.0
1675      F2=2.0
1676      PI2=6.283185307179586
1677      EPS=1.0e-6
1678
1679#-----
1680      ALP1=(F1-ALPHA)/F2
1681      ALP2= ALPHA/F2
1682      ALP3=(F1-ALPHA)/ALPHA
1683      ALP4= F1-ALPHA
1684      ALP5= ALPHA
1685      #print ALP1
1686#-----
1687      P18=PI2/360.0
1688      SD=sin(DIP*P18)
1689      CD=cos(DIP*P18)
1690      if abs(CD)<EPS:
1691        CD=F0
1692        if SD>F0: SD= F1
1693        if SD<F0: SD=-F1
1694      SDSD=SD*SD
1695      CDCD=CD*CD
1696      SDCD=SD*CD
1697      S2D=F2*SDCD
1698      C2D=CDCD-SDSD
1699
1700      self.ALP1=ALP1
1701      self.ALP2=ALP2
1702      self.ALP3=ALP3
1703      self.ALP4=ALP4
1704      self.ALP5=ALP5
1705      self.SD=SD
1706      self.CD=CD
1707      self.SDSD=SDSD
1708      self.CDCD=CDCD
1709      self.SDCD=SDCD
1710      self.S2D=S2D
1711      self.C2D=C2D
1712
1713    def DCCON1(self,X,Y,D):
1714      from math import sqrt
1715#     SUBROUTINE  DCCON1(X,Y,D)
1716 #     IMPLICIT REAL*8 (A-H,O-Z)
1717
1718      """**********************************************************************
1719      C*****   CALCULATE STATION GEOMETRY CONSTANTS FOR POINT SOURCE    *****
1720      C**********************************************************************
1721      C
1722      C***** INPUT
1723      C*****   X,Y,D : STATION COORDINATES IN FAULT SYSTEM
1724      C### CAUTION ### IF X,Y,D ARE SUFFICIENTLY SMALL, THEY ARE SET TO ZERO  """
1725
1726#      COMMON /C0/DUMMY(5),SD,CD,dumm(5)
1727#      COMMON /C1/P,Q,S,T,XY,X2,Y2,D2,R,R2,R3,R5,QR,QRX,A3,A5,B3,C3,
1728#     *           UY,VY,WY,UZ,VZ,WZ
1729      F0=0.0
1730      F1=1.0
1731      F3=3.0
1732      F5=5.0
1733      EPS=1.0e-6
1734
1735      SD=self.SD
1736      CD=self.CD
1737#-----
1738      if abs(X)<EPS: X=F0
1739      if abs(Y)<EPS: Y=F0
1740      if abs(D)<EPS: D=F0
1741      P=Y*CD+D*SD
1742      Q=Y*SD-D*CD
1743      S=P*SD+Q*CD
1744      T=P*CD-Q*SD
1745      XY=X*Y
1746      X2=X*X
1747      Y2=Y*Y
1748      D2=D*D
1749      R2=X2+Y2+D2
1750      R =sqrt(R2)
1751      if R==F0: return
1752      R3=R *R2
1753      R5=R3*R2
1754      R7=R5*R2
1755#-----
1756      A3=F1-F3*X2/R2
1757      A5=F1-F5*X2/R2
1758      B3=F1-F3*Y2/R2
1759      C3=F1-F3*D2/R2
1760#-----
1761      QR=F3*Q/R5
1762      QRX=F5*QR*X/R2
1763#-----
1764      UY=SD-F5*Y*Q/R2
1765      UZ=CD+F5*D*Q/R2
1766      VY=S -F5*Y*P*Q/R2
1767      VZ=T +F5*D*P*Q/R2
1768      WY=UY+SD
1769      WZ=UZ+CD
1770
1771
1772      self.P=P
1773      self.Q=Q
1774      self.S=S
1775      self.T=T
1776      self.XY=XY
1777      self.X2=X2
1778      self.Y2=Y2
1779      self.D2=D2
1780      self.R2=R2
1781      self.R=R
1782      self.R3=R3
1783      self.R5=R5
1784      self.R7=R7
1785      self.A3=A3
1786      self.A5=A5
1787      self.B3=B3
1788      self.C3=C3
1789      self.QR=QR
1790      self.QRX=QRX
1791      self.UY=UY
1792      self.UZ=UZ
1793      self.VZ=VZ
1794      self.VY=VY
1795      self.WY=WY
1796      self.WZ=WZ
1797
1798    def DCCON2(self,XI,ET,Q,SD,CD,KXI,KET):
1799      from math import sqrt,atan,log
1800#      SUBROUTINE  DCCON2(XI,ET,Q,SD,CD,KXI,KET)
1801#      IMPLICIT REAL*8 (A-H,O-Z)
1802
1803      """**********************************************************************
1804      C*****   CALCULATE STATION GEOMETRY CONSTANTS FOR FINITE SOURCE   *****
1805      C**********************************************************************
1806      C
1807      C***** INPUT
1808      C*****   XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM
1809      C*****   SD,CD   : SIN, COS OF DIP-ANGLE
1810      C*****   KXI,KET : KXI=1, KET=1 MEANS R+XI<EPS, R+ET<EPS, RESPECTIVELY
1811      C
1812      C### CAUTION ### IF XI,ET,Q ARE SUFFICIENTLY SMALL, THEY ARE SET TO ZER0"""
1813
1814
1815      F0=0.0
1816      F1=1.0
1817      F2=2.0
1818      EPS=1.0e-6
1819#-----
1820      if abs(XI)<EPS: XI=F0
1821      if abs(ET)<EPS: ET=F0
1822      if abs(Q)<EPS:  Q=F0
1823      XI2=XI*XI
1824      ET2=ET*ET
1825      Q2=Q*Q
1826      R2=XI2+ET2+Q2
1827      R =sqrt(R2)
1828      if R==F0: return
1829      R3=R *R2
1830      R5=R3*R2
1831      Y =ET*CD+Q*SD
1832      D =ET*SD-Q*CD
1833#-----
1834      if Q==F0:
1835        TT=F0
1836      else:
1837        TT=atan(XI*ET/(Q*R))
1838#-----
1839      if KXI==1:
1840        ALX=-log(R-XI)
1841        X11=F0
1842        X32=F0
1843      else:
1844        RXI=R+XI
1845        ALX=log(RXI)
1846        X11=F1/(R*RXI)
1847        X32=(R+RXI)*X11*X11/R
1848#-----
1849      if KET==1:
1850        ALE=-log(R-ET)
1851        Y11=F0
1852        Y32=F0
1853      else:
1854        RET=R+ET
1855        ALE=log(RET)
1856        Y11=F1/(R*RET)
1857        Y32=(R+RET)*Y11*Y11/R
1858#-----
1859      EY=SD/R-Y*Q/R3
1860      EZ=CD/R+D*Q/R3
1861      FY=D/R3+XI2*Y32*SD
1862      FZ=Y/R3+XI2*Y32*CD
1863      GY=F2*X11*SD-Y*Q*X32
1864      GZ=F2*X11*CD+D*Q*X32
1865      HY=D*Q*X32+XI*Q*Y32*SD
1866      HZ=Y*Q*X32+XI*Q*Y32*CD
1867
1868      self.XI2=XI2
1869      self.Q2=Q2
1870      self.R=R
1871      self.R2=R2
1872      self.R3=R3
1873      self.R5=R5
1874      self.Y=Y
1875      self.D=D
1876      self.TT=TT
1877      self.ALX=ALX
1878      self.ALE=ALE
1879      self.X11=X11
1880      self.Y11=Y11
1881      self.X32=X32
1882      self.Y32=Y32
1883      self.EY=EY
1884      self.EZ=EZ
1885      self.FY=FY
1886      self.FZ=FZ
1887      self.GY=GY
1888      self.GZ=GZ
1889      self.HY=HY
1890      self.HZ=HZ
1891      self.ET2=ET2
1892
1893
1894
1895
1896
1897
1898
Note: See TracBrowser for help on using the repository browser.