source: anuga_core/source/anuga/shallow_water/tsunami_okada.py @ 5485

Last change on this file since 5485 was 5252, checked in by herve, 16 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#
86# Okada class
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
109class Okada_func:
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
205        #rad=1.745329252e-2
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):
225                #st = radians(strikes[ist])
226                st=radians(strikes[ist])
227                csst = cos(st)
228                ssst = sin(st)
229                cs2st = cos(2.0*st)
230                ss2st = sin(2.0*st)
231#
232                #di = radians(dips[ist])
233                di=radians(dips[ist])
234                csdi =cos(di)
235                ssdi=sin(di)
236#
237                #ra= radians(rakes[ist])
238                ra=radians(rakes[ist])
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
345        # f2py -c okada.f -m okada
346        #z1 = okada(xr,yr,depth,length,width,dip,rake,slip)
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.