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

Last change on this file since 5224 was 5206, checked in by herve, 16 years ago

modification of test_tsunami_okada.py to make txt file available.

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