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

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

following Rajaraman advices to make the code more efficient

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