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

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

Replaced 'print' statements with log.critical() calls.

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