source: branches/numpy/anuga/shallow_water/tsunami_okada.py @ 7208

Last change on this file since 7208 was 7208, checked in by rwilson, 15 years ago

Fixed problems with back-merge.

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