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

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

more explanations on input

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