source: anuga_core/source_numpy_conversion/anuga/shallow_water/tsunami_okada.py @ 5901

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

NumPy? conversion.

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