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