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

Last change on this file since 6410 was 6410, checked in by rwilson, 16 years ago

numpy changes.

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