AlllOb 05372b

N3SIR 85-3234

Dynamic Green's Functions of an

Infinite Plate - A Computer

Program

Reference PUBUCAT| 0ns

Nelson N. Hsu

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Manufacturing Engineering

Gaithersburg, MD 20SS9

August 1 985

issued November 1985

Prepared for:

. ,,o p ari-ment OF COMMERCE

Bureau of Standards

100 aurg, MD 20899

• U56

85-3234

1985

L , v . I %

NBSIR 85-3234

DYNAMIC GREEN'S FUNCTIONS OF AN

INFINITE PLATE - A COMPUTER

PROGRAM

Nelson N. Hsu

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Manufacturing Engineering

Gaithersburg, MD 20899

August 1 985

Issued November 1985

Prepared for: U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Gaithersburg, MD 20899

U.S. DEPARTMENT OF COMMERCE, Malcolm Baldrige, Secretary

NATIONAL BUREAU OF STANDARDS. Ernest Ambler. Director

national bureau

OF STANDARDS

XJBRABY

/

J

Dynamic Green's Functions of an Infinite Plate - a Computer Program

Nelson N. Hsu

National Bureau of Standards

Gaithersburg, MD 20899 USA

Introduction

This report is a FORTRAN program to compute the Green's functions of an infinite plate. The Green's function, G. . (£,x,t), is defined as the ith component of the displacement at x due to anoint force of step-function time

dependency acting at in the jth direction initiated at t=0. The Green's

function is the fundamental solution of the transient elastic wave propagation

problem. In general, the displacement field u(^,x,t) at x due to a point force

of arbitrary time dependence acting at can be computed by a convolution

integration; i . e.

,

u. U,x,t

)

COG. . (j(,x ,t)f . ( T-t )dt. (1

)

Here, G. . is the time derivative of G. . and f.(t) is the point force component

of arbi^ary time dependence acting inJ the jtA direction (summation over repeated indices is used). Displacement due to point dipoles or couple forces

can be represented by the spatial derivatives of G. .. Displacement produced by

a dynamic force distributed over a finite area can

1 ^lso be computed by numerical

integration using the Green's function as the kernel of the integral over the

finite area.

The basic formulation of the problem and derivation of the solution for

mulas were reported in (1). The method used is called "ray theory" in the

seismological literature. Our derivation was based upon John Willis' new

Fourier Inversion method (2). Similar computation results can be found in

Reference (3), (4), and (5).

The program was originally written for the analysis of acoustic emission

signals. Acoustic emission is a nondestructive testing and monitoring technique

in which the detection of the transient stress wave produced by localized defor

mation or cracks are used to locate, and to assess the criticality of the

defects. The theoretical computation is an important link to predict how the

acoustic emission waveform evolves through the structure. How the theoretical

computations are used in the study of acoustic emission can be found in

References (1), (6), (7), and (8).

This computer program is made available mainly for its application to

calibrate acoustic emission systems and sensors. By making the present program

available, duplication of efforts can be avoided, errors in theory can be

checked, and experimental results can be reproduced and verified.

How to Use the Program

I. The easiest way to run this program is by setting up nondimensional

parameters and calling the subroutine GREENFCT. First, the x-y-z coordinate.-'
are selected by choosing an x-y plane parallel to the plate and the origin at

the center of the plate directly underneath the source, i.e., the source is

always located at £=(0,0,0. 5); the x axis is aligned in the direction from the

source pointing to the detector. The required nondimensionalized input

parameters for the subroutine GREENFCT are:

1 . ALPHA = shear wave speed/longitudinal wave speed;

2. XD = x-coordinate of the detector; actual length/thickness of the plate;

3. ZD = z-coordinate of the detector; actual length/thickness of the plate.

ZD = 0.5 if the source and the detector are on the same side of the

plate

ZD = -0.5 if they are on opposite sides of the plate.

(Note: the y-coordinate of the detector is always zero because of the way the

coordinate system is chosen.)

M. INDEX = ij ; subscript of Green's function, integer number 11, 12, 13, 21,

22, 23, 31, 32, or 33- (111, 112, 113, etc. for force dipoles.)

5. TDELTA = sampling time interval in terms of nondimensionalized time unit

= actual sampling time interval * shear wave speed/thickness of the

plate; and

6. NP0INT = Total number of sampling points to be computed; must be an integer.

The subroutine will return a double precision array DISPL of dimension

(NP0INT) which corresponds to the desired G (£,x,t) or G„ (£,x,t) sampled at

equal time intervals. However, the displacements corresponding to dipoles or couples (i.e., G ^) are for linear ramp time dependency input.

Dif ferentitatingi m?ie returned DISPL with respect to time once will give the

proper Green's function due to step time dependency input.

A simple program calls GREENFCT for given nondimensionalized input paramters together with the output results is included in Appendix A; a program that

prompts the user for input parameters of arbitrary physical units and returns

displacement in physical units is included in Appendix B. The complete listing

of the subroutines is in Appexdix C.

II. Another way to use this program is by calling GREENSUB which is basically a simple function: for given one time T, it returns a displacement DISPL.

However, there are two preparatory steps that must be taken in the main program:

(1) RAYTIME (NRAY , 3), the time of arrival table, TA(NRAY) and CN (NRAY , 3),

two working arrays, should be dimensioned.

(2) The subroutines INIT and TIMEARRI must be called first.

The subroutine GREENFCT may serve as an example of how to use the subroutine GREENSUB.

III. In addition to the two subroutines mentioned above, there are three subroutines which may be called independently and the user may find them useful in

checking experimental or theoretical results.
(1) EPIDIS computes the vertical displacement at the epicenter for given

ALPHA and T.

(2) SFWAVE computes the surface wave, which corresponds to the Green's function of a semi-infinite space.

(3) TIMEARRI computes the time of arrival for various ray paths.

Comments in these three subroutines serve as instructions about how to use these subroutines.

Some Remarks

I. There are two machine dependent constants which are used in the two numerical integration routines DGLQI and GLI5T. EPMACH is the largest relative

spacing. UFLOW is the smallest positive magnitude. See Reference (9) for

details about how to set these constants for different computers.

II. The program uses double precision complex arithmetic. In order to run on those FORTRAN compilers lacking double precision complex function libraries, the user may have to modify the program to single precision arithmetic.

III. The ray method is an exact solution for the Green's function in the time

domain produced by summing the contribution from successive arrivals of

reflected rays. The numerical computation is accurate if the source and detector are near each other (say less than ten plate thicknesses) and if the maximum

time is less than ten dimensionless units. If the sensor and detector are far

apart or the maximum time is large, the number of rays arriving at approximately

the same time may be too large to compute in a reasonable time and the accumulated error may grow.

IV. The current program is limited to the test configuration where the sensor

and detector are on the surface of the plate. The displacement computed for

force dipole (i.e., INDEX = 111, 112, etc.) are for point dipole of linear ramp

time dependency. Differentiation with respect to time once of the displacement

will produce proper Green's functions due to step time dependency (See Reference

10 ).

Acknowledgement

The author would like to thank Dr. David Kahaner of the Scientific

Computing Division, NBS, for supplying the integration subroutines DGLQ and

GLIST.

References :

1. D.G. Eitzen, et. al., "Fundamental Developments for Quantitative Acoustic

Emission Measurements," EPRI NP-2089 , Project 608-1, Palo Alto, CA:

Electric Power Research Institute, Oct. 1981.

2. J.R. Willis, "Self-Similar Problems in Elastodynamics , " Phil. Trans, of Roy.

Soc. London, 274, pp. 435-491, May 1973-
3 . Y. Pao, Ralph R. Gajewski and Ahmet N. Ceranouglu, "Acoustic Emission and

Transient Waves in an Elastic Plate," J. Acoust. Soc. Am., 65(1), Jan.

1979.

4. A.N. Cernoglu and Y.H. Pao, "Propagation of Eleastic Pulses and Acoustic

Emission in a Plate," ASME Journal of Applied Mechanics, 48, pp. 125-147.

Mar. 1981.

5. T.M. Proctor, F.R. Breckenridge and Y.H. Pao, "Transient Waves in an Elastic

Plate: Theory and Experiment Compared," J. Acoust. Soc. Am., (6) Dec. 1983.

6. N.N. Hsu, J.A. Simmons and S.C. Hardy, "An Approach to Acoustic Emission

Signal Analysis — Theory and Experiment," Materials Evaluation, 35, pp.

100-106, Oct. 1977.

7. N.N. Hsu and S.C. Hardy, "Experiments in Acoustic Emission Waveform Analysis

for Characterization of AE Sources, Sensors and Structures; In Elastic Waves

and Non-Destructive Testing of Materials," AMD-Vol. 29, Y.H. Pao, Editor,

The American Society of Mechanical Engineers, New York, pp. 85-106, Dec.

1978.

8. C. Chang and W. Sachse, "Analysis of Elastic Wave Signals from An Extended

Source in A Plate," J. Acoust. Soc. Am., 77(4), April 1985.

9. P . A . Fox, A . D . Hall and N.L. Schryer, Bell Laboratories Computer Science

Technical Report No. 47. Or: The PORT Mathematical Subroutine Library, ACM

Transactions on Mathematical Software, Vol. 4, No. 2, pp 104-106.

10. J.E. Sinclair, "Epicentre Solutions for Point Multipole Sources in An

Elastic Half-Space," J. Phys. D: Appl. Phys., Vol. 12, 1979, pp 1309-1315.
Appendix A

Sample Program and Result

^****************************************

C

Cc

c

cc

c

c

c

c

PROGRAM TSGREEN

TEST PROGRAM TO COMPUTE GREEN'S FUNCTION OF A PLATE

DOUBLE PRECISION DISPLC<256)

URITE(5,*) 'TYPE INDEX, ALPHA, XD, AND ZD' READ(5,*) INDEX, ALPHA, XD, ZD INDEX=33 ALPHA=0. 59902

XD=1 .0 ZD=-0.5

WRITE(5,*) 'TYPE DELTAT AND NUMBER OF POINTS' READ(5,*) TDELTA,NPOINT

TDELTA=0.01

NPOINT=256

CALL GREENFCTCALPHA, XD, ZD, INDEX, TDELTA,NPOINT,DISPLC)

DO 1000 1=1,256

TIME = 0.01*REAL(I

)

WRITEC6, 1001 ) TIME,DISPLC( I

)

1000 CONTINUE

1001 FORMAT(2X, 'TIME = ',F8.3,', STOP

END

DISPLACEMENT = \D21.14)
TIME = 0.010, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.020, DISPLACEMENT = 0.00000000000000D+00

TIME = 0.030, DISPLACEMENT = 0.00000000000000D+00

TIME = 0.040, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.050, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.060, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.070, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.080, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.090, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.100, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.110, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.120, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.130, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.140, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.150, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.160, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.170, DISPLACEMENT = O.OOOOOCOOOOOOOOD+OO

TIME = 0.180, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.190, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.200, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.210, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.220, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.230, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.240, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.250, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.260, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.270, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.280, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.290, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME - 0.300, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.310, DISPLACEMENT = 0.000000000000000+00

TIME = 0.320, DISPLACEMENT r O.OOOOOOOOOOOOOOD+OO

TIME = 0.330, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.340, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.350, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.360, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.370, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.380, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.390, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.400, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.410, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.420, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.430, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.440, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.450, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.460, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.470, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.480, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.490, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.500, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.510, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.520, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.530, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.540, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO
TIME = 0.550, DISPLACEMENT = 0.00000000000000D+00

TIME = 0.560, DISPLACEMENT = 0.00000000000000D+00

TIME = 0.570, DISPLACEMENT = 0.00000000000000D+00

TIME = 0.580, DISPLACEMENT = 0.00000000000000D+00

TIME = 0.590, DISPLACEMENT = 0.00000000000000D+00

TIME = 0.600, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.610, DISPLACEMENT = 0.00000000000000D+00

TIME = 0.620, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.630, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.640, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.650, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.660, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.670, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME - 0.680, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.690, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.700, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.710, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.720, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.730, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.740, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.750, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.760, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.770, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.780, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.790, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.800, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.810, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.820, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.830, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.840, DISPLACEMENT = O.OOOOOOOOOOOOOOD+OO

TIME = 0.850, DISPLACEMENT = 0.11 07506769871 2D+00

TIME = 0.860, DISPLACEMENT = 0.11 0986061 56782D+00

TIME = 0.870, DISPLACEMENT = 0.111 1 8024981 467D+00

TIME = 0.880, DISPLACEMENT = 0.11 133277774482D+00

TIME = 0.890, DISPLACEMENT = 0.11 14431 8482884D+00

TIME = 0.900, DISPLACEMENT = 0.1 1 151 1026521 21D+00

TIME = 0.910, DISPLACEMENT = 0.11 1535876944 15D+00

TIME = 0.920, DISPLACEMENT = 0.11 15173360582 1D+00

TIME = 0.930, DISPLACEMENT = 0.1 11455035416960+00

TIME = 0.940, DISPLACEMENT = 0.11 134864281 672D+00

TIME = 0.950, DISPLACEMENT = 0.1111 9787020707D+00

TIME = 0.960, DISPLACEMENT = 0.11 100247472 178D+00

TIME = 0.970, DISPLACEMENT = 0. 1 1076226552545D+00

TIME = 0.980, DISPLACEMENT = 0.1 104771 076421 9D+00

TIME = 0.990, DISPLACEMENT = 0. 1 1 014692394960D+00

TIME = 1.000, DISPLACEMENT = 0.1 09771 71 429770D+00

TIME = 1.010, DISPLACEMENT = 0. 10935 1 509141 28D+00

TIME = 1.020, DISPLACEMENT = 0.1 08886450335 15D+00

TIME = 1.030, DISPLACEMENT = 0 . 1 0837678979646D+00

TIME = 1.040, DISPLACEMENT = 0 . 1 0782270845288D+00

TIME = 1.050, DISPLACEMENT = 0.1 0722461 062043D+00

TIME = 1.060, DISPLACEMENT = 0 . 1 0658293623297D+00

TIME = 1.070, DISPLACEMENT = 0 . 1 0589820949205D+00

TIME = 1.080, DISPLACEMENT = 0. 105171 1121 083 1 D+00

TIME = 1.090, DISPLACEMENT = 0 . 1 0440220346099D+00
TIME = 1.100, DISPLACEMENT = 0 . 1 0359233809357D+00

TIME = 1.110, DISPLACEMENT = 0 . 1 0274239332029D+00

TIME - 1.120, DISPLACEMENT = 0. 10185342388500D+00

TIME = 1.130, DISPLACEMENT = 0.1 009263 1946377D+00

TIME = 1.140, DISPLACEMENT = 0. 999623 13600975D- 01 TIME = 1.150, DISPLACEMENT = 0.98962647939981D- 01 TIME = 1.160, DISPLACEMENT = 0 . 9792875478595 1 D 01 TIME = 1.170, DISPLACEMENT = 0. 96861 858478749D- 01 TIME = 1.180, DISPLACEMENT = 0. 957635660 16467D- 01 TIME = 1.190, DISPLACEMENT = 0 . 94635474207730D - 01 TIME = 1.200, DISPLACEMENT = 0 . 93479376354866D - 01 TIME = 1.210, DISPLACEMENT = 0 . 9229681 6327325D -01 TIME = 1.220, DISPLACEMENT = 0.91 089747924766D - 0

1

TIME = 1.230, DISPLACEMENT = 0.89860093569871D-01

TIME = 1.240, DISPLACEMENT = 0 .88609970379555D - 01 TIME = 1.250, DISPLACEMENT = 0. 87341 207961 754D- 01 TIME = 1.260, DISPLACEMENT = 0 . 86056060561 906D - 01 TIME = 1.270, DISPLACEMENT = 0. 847567281 27925D- 01 TIME = 1.280, DISPLACEMENT = 0. 834455961 72436D- 01 TIME = 1.290, DISPLACEMENT = 0 . 821 24730252424D - 01 TIME = 1.300, DISPLACEMENT = 0 . 80796622274542D - 01 TIME = 1.310, DISPLACEMENT = 0 . 79463686723 1 94D - 0

1

TIME = 1.320, DISPLACEMENT = 0.781 28380592089D - 01 TIME = 1.330, DISPLACEMENT = 0 . 76793325943265D - 01 TIME = 1.340, DISPLACEMENT = 0 . 75460795794622D - 01 TIME = 1.350, DISPLACEMENT = 0. 74133473644589D- 01 TIME = 1.360, DISPLACEMENT = 0. 7281 393901 21 66D- 01 TIME = 1.370, DISPLACEMENT = 0. 71 50491 4223321 D- 01 TIME = 1.380, DISPLACEMENT = 0.70208761 123235D-01

TIME = 1.390, DISPLACEMENT = 0.68928225721454D-01

TIME = 1.400, DISPLACEMENT = 0 . 63777644085523D - 0

1

TIME = 1.410, DISPLACEMENT = -0. 31205499826477D-01 TIME = 1.420, DISPLACEMENT = 0. 36764444097372D+00

TIME = 1.430, DISPLACEMENT = 0.391 058631 5 1275D+00

TIME = 1.440, DISPLACEMENT = 0. 40949206785671 D+00

TIME = 1.450, DISPLACEMENT = 0. 4258496636743 1D+00

TIME = 1.460, DISPLACEMENT = 0.441 02467848492D+00

TIME = 1.470, DISPLACEMENT = 0 . 45540362930489D+00

TIME = 1.480, DISPLACEMENT = 0. 4691 8701 3091 17D+00

TIME = 1.490, DISPLACEMENT = 0. 4824881 0354954D+00

TIME = 1.500, DISPLACEMENT = 0.4953792919241 1D+00

TIME = 1.510, DISPLACEMENT = 0 . 50790390548020D+00

TIME = 1.520, DISPLACEMENT = 0.52009143324477D+00

TIME = 1.530, DISPLACEMENT = 0.53 1961 91 904305D+00

TIME = 1.540, DISPLACEMENT = 0.543528255164190+00

TIME = 1.550, DISPLACEMENT = 0. 554802671341 78D+00

TIME = 1.560, DISPLACEMENT = 0. 565791 47374873D+00

TIME = 1.570, DISPLACEMENT = 0 . 57650038247907D+00

TIME = 1.580, DISPLACEMENT = 0. 5869330375 1993D+00

TIME = 1.590, DISPLACEMENT = 0. 597095366361 70D+00

TIME = 1.600, DISPLACEMENT = 0.6069898995 10390+00

TIME = 1.610, DISPLACEMENT = 0 .61661 997278485D+00

TIME = 1.620, DISPLACEMENT = 0.62598797275653D+00

TIME = 1.630, DISPLACEMENT = 0 . 635 098976464060+00

TIME = 1.640, DISPLACEMENT = 0. 64395545056472D+00
TIME = 1.650, DISPLACEMENT = 0. 65256088669347D+00

TIME = 1.660, DISPLACEMENT = 0.66091 81 3088975D+00

TIME = 1.670, DISPLACEMENT = 0.6690325571 71 67D+00

TIME = 1.680, DISPLACEMENT = 0. 67690728626371 D+00

TIME = 1.690, DISPLACEMENT = 0. 684546391 70679D+00

TIME = 1.700, DISPLACEMENT = 0 . 691 95338607382D+00

TIME = 1.710, DISPLACEMENT = 0.6991 3401 239547D+00

TIME = 1.720, DISPLACEMENT = 0.706091 992259960+00

TIME = 1.730, DISPLACEMENT = 0.71 2831 858301 27D+00

TIME = 1.740, DISPLACEMENT = 0.71 9358221 54785D+00

TIME = 1.750, DISPLACEMENT = 0. 725675 16353530D+00

TIME = 1.760, DISPLACEMENT = 0.731788611 32 189D+00

TIME = 1.770, DISPLACEMENT = 0.737702708981 97D+00

TIME = 1.780, DISPLACEMENT = 0 . 7434222332603 1 D+00

TIME = 1.790, DISPLACEMENT = 0. 74895 144993679D+00

TIME = 1.800, DISPLACEMENT = 0.754296 19306931 D+00

TIME = 1.810, DISPLACEMENT = 0 . 75946069053969D+00

TIME = 1.820, DISPLACEMENT = 0 . 764449680843430+00

TIME = 1.830, DISPLACEMENT = 0 . 7692674 1601 862D+00

TIME = 1.840, DISPLACEMENT = 0.77391947485854D+00

TIME = 1.850, DISPLACEMENT = 0.7784099991 1627D+00

TIME = 1.860, DISPLACEMENT = 0 . 78274354050484D+00

TIME = 1.870, DISPLACEMENT = 0. 786924 19559247D+00

TIME = 1.880, DISPLACEMENT = 0. 7909571 8088098D+00

TIME = 1.890, DISPLACEMENT = 0 . 79484643653490D+00

TIME = 1.900, DISPLACEMENT = 0.86806930229072D+00

TIME = 1.910, DISPLACEMENT = 0 . 87369508356078D+00

TIME = 1.920, DISPLACEMENT = 0 . 879233 1 4459302D+00

TIME = 1.930, DISPLACEMENT = 0 . 884687498388460+00

TIME = 1.940, DISPLACEMENT = 0.89006261 168853D+00

TIME = 1.950, DISPLACEMENT = 0 . 89536286267026D+00

TIME = 1.960, DISPLACEMENT = 0. 900592044445 14D+00

TIME = 1.970, DISPLACEMENT = 0.905755352541 14D+00

TIME = 1.980, DISPLACEMENT = 0.91 085639280247D+00

TIME = 1.990, DISPLACEMENT = 0.91 58991 8253638D+00

TIME = 2.000, DISPLACEMENT = 0.9208871 7546076D+00

TIME = 2.010, DISPLACEMENT = 0 . 925825 1 5984505D+00

TIME = 2.020, DISPLACEMENT = 0.93071 640538221 D+00

TIME = 2.030, DISPLACEMENT = 0.935564571 71 401 D+00

TIME = 2.040, DISPLACEMENT = 0. 94037277440757D+00

TIME = 2.050, DISPLACEMENT = 0.9451 454 1664089D+00

TIME = 2.060, DISPLACEMENT = 0 . 94988543954263D+00

TIME = 2.070, DISPLACEMENT = 0.9545961 61 31 288D+00

TIME = 2.080, DISPLACEMENT = 0.95928037273261 D+00

TIME = 2.090, DISPLACEMENT = 0.9639421 23260 14D+00

TIME = 2.100, DISPLACEMENT = 0.96858404264596D+00

TIME = 2.110, DISPLACEMENT = 0.97320913182266D+00

TIME = 2.120, DISPLACEMENT = 0. 97781 987758092D+00

TIME = 2.130, DISPLACEMENT = 0. 98242001 202845D+00

TIME = 2.140, DISPLACEMENT = 0.98701 187638821 D+00

TIME = 2.150, DISPLACEMENT = 0.991 5981 8288714D+00

TIME = 2.160, DISPLACEMENT = 0 . 9961 81 57630830D+00

TIME = 2.170, DISPLACEMENT = 0.1 0007641 9879 12D+01 TIME = 2.180, DISPLACEMENT = 0. 10053494399595D+01

TIME = 2.190, DISPLACEMENT = 0.1 0099393 168333D+01
TIME = 2.200, DISPLACEMENT = 0. 10145362237750D+01

TIME = 2.210, DISPLACEMENT = 0.10191420575858D+01

TIME = 2.220, DISPLACEMENT = 0.1 02375 99774067D+01

TIME = 2.230, DISPLACEMENT = 0.1 028391 771 2906D+01

TIME = 2.240, DISPLACEMENT = 0.1 03303961 40693D+01

TIME = 2.250, DISPLACEMENT = 0.1 037705 1845379D+01

TIME = 2.260, DISPLACEMENT = 0.1 042391 4489334D+01

TIME = 2.270, DISPLACEMENT = 0.1 047099993231 4D+01 TIME = 2.280, DISPLACEMENT = 0. 10518328029171D+01

TIME = 2.290, DISPLACEMENT = 0. 10565913655695D+01

TIME = 2.300, DISPLACEMENT = 0.106137849050690+01

TIME = 2.310, DISPLACEMENT = 0.1 1476631 005725D+01

TIME = 2.320, DISPLACEMENT = 0. 1 1533185695767D+01

TIME = 2.330, DISPLACEMENT = 0. 1 1590500263745D+01

TIME = 2.340, DISPLACEMENT = 0.11648631801910D+01

TIME = 2.350, DISPLACEMENT = 0.11 707621 604698D+01

TIME = 2.360, DISPLACEMENT = 0.11 7675 16974758D+01

TIME = 2.370, DISPLACEMENT = 0.1 182836581 1738D+01

TIME = 2.380, DISPLACEMENT = 0. 11890210681560D+01

TIME = 2.390, DISPLACEMENT = 0.1 19531 12520959D+01

TIME = 2.400, DISPLACEMENT = 0. 120171 15269290D+01

TIME = 2.410, DISPLACEMENT = 0. 12082269375580D+01

TIME = 2.420, DISPLACEMENT = 0. 12148619575064D+01

TIME = 2.430, DISPLACEMENT = 0.1 221 62303280 19D+01 TIME = 2.440, DISPLACEMENT = 0.1 2285 147870224D+01

TIME = 2.450, DISPLACEMENT = 0. 12355425415795D+01

TIME = 2.460, DISPLACEMENT = 0.124271099911260+01

TIME = 2.470, DISPLACEMENT = 0 . 1 2500269925297D+01

TIME = 2.480, DISPLACEMENT = 0.1 2574953844797D+01

TIME = 2.490, DISPLACEMENT = 0.12651217887352D+01

TIME = 2.500, DISPLACEMENT = 0. 127291 11430751D+01

TIME = 2.510, DISPLACEMENT = 0. 12808706969625D+01

TIME = 2.520, DISPLACEMENT = 0. 12890055553986D+01

TIME = 2.530, DISPLACEMENT = 0.1 29732 16346622D+01

TIME = 2.540, DISPLACEMENT = 0. 1305824 1075455D+01

TIME = 2.550, DISPLACEMENT = 0. 13145206649727D+01

TIME = 2.560, DISPLACEMENT = 0.1 32341 665 14682D+01 TIME = 2.080, DISPLACEMENT = 0.95928037273261 D+00

TIME = 2.090, DISPLACEMENT = 0 . 9639421 232601 4D+00

TIME = 2.100, DISPLACEMENT = 0 . 96858404264596D+00

TIME = 2.110, DISPLACEMENT = 0.973209131822660+00

TIME = 2.120, DISPLACEMENT = 0. 97781 987758092D+00

TIME = 2.130, DISPLACEMENT = 0. 98242001 202845D+00

TIME = 2.140, DISPLACEMENT = 0.98701 187638821 D+00

TIME = 2.150, DISPLACEMENT = 0.991 5981 828871 4D+00

TIME = 2.160, DISPLACEMENT = 0 . 996181 57630830D+00

TIME = 2.170, DISPLACEMENT = 0.1 0007641 98791 2D+01 TIME = 2.180, DISPLACEMENT = 0. 10053494399595D+01

TIME = 2.190, DISPLACEMENT = 0. 10099393168333D+01

TIME = 2.200, DISPLACEMENT = 0. 1014536223775 0D+01 TIME = 2.210, DISPLACEMENT = 0. 10191420575858D+01

TIME = 2.220, DISPLACEMENT = 0. 102375997740670+01

TIME = 2.230, DISPLACEMENT = 0.1 028391 771 2906D+01

TIME = 2.240, DISPLACEMENT = 0. 10330396140693D+01

TIME = 2.250, DISPLACEMENT = 0.1 037705 1845379D+01
TIME = 2.260, DISPLACEMENT = 0.1 042391 4489334D+01

TIME = 2.270, DISPLACEMENT = 0.1 04709999323 14D+01 TIME = 2.280, DISPLACEMENT = 0.1 05 183280291 71D+01

TIME = 2.290, DISPLACEMENT = 0. 10565913655695D+01

TIME = 2.300, DISPLACEMENT = 0. 10613784905069D+01

TIME = 2.310, DISPLACEMENT = 0.1 1476631005725D+01

TIME = 2.320, DISPLACEMENT = 0. 1 1533185695767D+01

TIME = 2.330, DISPLACEMENT = 0.11 590500263745D+01

TIME = 2.340, DISPLACEMENT = 0.1 1648631801910D+01

TIME = 2.350, DISPLACEMENT = 0.11 70762 16046980+01

TIME = 2.360, DISPLACEMENT = 0.11 7675 16974758D+01

TIME = 2.370, DISPLACEMENT = 0.1 182836581 17380+01

TIME = 2.380, DISPLACEMENT = 0.1 1890210681560D+01

TIME = 2.390, DISPLACEMENT = 0.1 19531 12520959D+01

TIME = 2.400, DISPLACEMENT = 0. 120171 15269290D+01

TIME = 2.410, DISPLACEMENT = 0 . 1 2082269375580D+01

TIME = 2.420, DISPLACEMENT = 0. 12148619575 064D+01

TIME r 2.430, DISPLACEMENT = 0.1 221 62303280 190+01 TIME = 2.440, DISPLACEMENT = 0. 12285147870224D+01

TIME = 2.450, DISPLACEMENT = 0. 12355425415795D+01

TIME = 2.460, DISPLACEMENT = 0.1 24271 0999 11260+01 TIME = 2.470, DISPLACEMENT = 0.12500269925297D+01

TIME = 2.480, DISPLACEMENT = 0. 12574953844797D+01

TIME = 2.490, DISPLACEMENT = 0. 12651217887352D+01

TIME = 2.500, DISPLACEMENT = 0. 127291 11430751D+01

TIME = 2.510, DISPLACEMENT = 0. 12808706969625D+01

TIME = 2.520, DISPLACEMENT = 0. 12890055553986D+01

TIME = 2.530, DISPLACEMENT = 0.1 29732 16346622D+01

TIME = 2.540, DISPLACEMENT = 0. 13058241075455D+01

TIME = 2.550, DISPLACEMENT = 0. 13145206649727D+01 TIME = 2.560, DISPLACEMENT = 0. 13234166514682D+01

Appendix B

Another sample program - use physical units rather nondimensionalized numbers.

PROGRAM GPLATE DOUBLE PRECISION DISPV(1024) DIMENSION T I ME ( 1024) INTEGER YORNO

C URITEO ,*) 'THIS PROGRAM WILL COMPUTE THE DISPLACEMENT U SUB I AT WRITEO ,*) ' A POINT ON THE SURFACE OF AN INFINITE PLATE DUE TO A URITEO ,*) 1 POINT STEP IMPULSIVE FORCE F SUB J OR A POINT LINEAR WRITEO ,*) 1 RAMP FORCE DIPOLE F SUB J,K ACTING AT ANOTHER POINT 1 WRITEO,*)' ON THE SURFACE OF THE PLATE. 1 WRITEO,*)' INPUTS REQUIRED ARE PROMPTED AND TYPED IN ON THE ' WRITEO,*)' KEYBOARD.' WRITEO,*)' '

10 WRITEO,*)' SUBSCRIPTS OF G = ? O. e. 11,13,131,322, etc.)'

READ(1,*) INDEX WRITEO,*)' SHEAR WAVE SPEED = ?' READ(1,*) SHSPD WRITEO,*)' LONGITUDINAL WAVE SPEED = ?' READ(1,*) PSPD

100 WRITEO,*) 'ARE SOURCE AND DETECTER ON THE SAME SIDE OF THE PLATE WRITEO,*)' PLEASE TYPE 1 FOR YES' READ (1,*) YORNO

IF (YORNO. EQ.1) THEN

ZD = 0.5 ELSE

ZD = -0.5 END IF WRITEO,*)' WHAT IS THE THICKNESS OF THE PLATE ?' READ(1,*) HTHICK WRITEO,*)' WHAT IS THE DISTANCE BETWEEN DETECTER AND SOURCE' WRITEO,*)' OR EPICENTER ?' READO,*) X URITEO,*)' SAMPLING TIME INTERVAL = ?' READO,*) DT WRITEO,*)' NUMBER OF SAMPLING POINTS = ? (MAX = 1024)'

READO,*) NPT WRITEO,*)' SHEAR MODULUS OF THE PLATE = ?' READO,*) SHMDUL WRITEO,*)' THE OUTPUT WILL BE TIME VERSUS DISPLACEMENT.' WRITEO,*)' '

C ALPHA = SHSPD/PSPD

XD = X/HTHICK

TSCALE = HTHICK/SHSPD

DELTAT = DT/TSCALE

DSCALE = 1 . / ( 2

. *AS I N ( 1 .0)*SHMDUL*HTHICK)

C WRITE(2,*) ' ' URITE(2,*)' THE FOLLOWINGS ARE INPUT PARAMETERS:' WRITE(2,*) ' INDEX = '.INDEX WRITE(2,*) ' SHEAR WAVE SPEED = ', SHSPD WRITE(2,*)' LONGITUDINAL WAVE SPEED = ',PSPD WRITE (2,*) 'THICKNESS OF THE PLATE = ' , HTHICK

WRITE(2,*) 'SAMPLING TIME INTERVAL = ' ,DT
WRITE(2,*) 1 SHEAR MODULUS = ' ,SHMDUL WRITE (2,*) 'NUMBER OF POINTS TO BE COMPUTED = ',NPT

IF (YORNO.EQ. 1 ) THEN WRITE(2,*) ' THE DETECTER AND THE SOURCE ARE ON THE SAME SIDE 1 ELSEWRITE(2,*)'THE DETECTER AND THE SOURCE ARE ON OPPOSITE SIDES' END IF

WRITE(2,*) ' 1 CALL GREEN FCT( ALPHA, XD, ZD, INDEX,DELTAT,NPT,DISPV)

DO 200 1=1, NPT

TIME( I )=DELTAT*REAL( I )*TSCALE

D I SPV( I )=DI SPV( I )*DSCALE

WRITE(2,*) 1 TIME = ',TIME(I),' DISPLACEMENT = ',DISPV(I)

200 CONTINUE

STOP

END
Appendix C

Listing of Subroutines for Computing Green's function of an infinite plate.

c GR000020

SUBROUTINE GREENFCT(ALPHA f XD, ZD, INDEX, TDELTA, NPOINT, DISPL) GR000030

c GR000040

q*******************-************************************************ GR000050

c GR000060

C TO COMPUT GREEN'S FUNCTION - DISPLACEMENT AS A FUNCTION OF TIME GR000070

C - FOR GIVEN TEST CONFIGURATION GR000080

C GR000090

C INPUTS: GR000100

C ALPHA = SHEAR WAVE SPEED / LONGITUDINAL WAVE SPEED GR000110

C XD = X COORDINATE OF THE DETECTOR GR000120

C ZD = Z COORDINATE OF DETECTOR GR000130

C ZD = 0.5 => DETECTOR & SOURCE ON THE SAME SIDE GR000140

C ZD =-0.5 => DETECTOR & SOURCE ON OPPOSIT SIDES GR000150

C INDEX = SUBSCRIPT OF GREEN'S FUNCTION, i.e. 33, 322 GR000160

C TDELTA = TIME INCREMENT GR000170

C NPOINT = TOTAL NUMBER OF POINTS TO BE COMPUTED GR000180

C GR000190

C OUTPUT: GR000200

C DISPL = DISPLACEMENT ** DOUBLE PRECISION, MUST BE GR000210

C DIMENSIONED (NPOINT) GR000220

C GR000230

C GR000240

C NOTE 1: GR000250

C ALL INPUTS AND OUTPUTS ARE NOND I MENS I ONAL I ZED PARAMETERS GR000260

C DISPL = NORMALIZED DISPLACEMENT = ACTUAL DISPLACEMENT * PI GR000270

C * SHEAR MODULUS * PLATE THICKNESS / FORCE GR000280

C T = NORMALIZED TIME = ACTUAL TIME *SHEAR WAVE SPEED / PLATE GR000290

C THICKNESS GR000300

C DISPL MUST BE DIMENSIONED IN THE CALLING PROGRAM (NPOINT) GR000310

C GR000320

C GR000330

C GR000340

C NOTE 2: GR000350

C THE ARRIVAL TIMES ARE STORED IN RAYTIME(N, I ) GR000360

C N = THE Nth RAY GR000370

C RAYTIME(N,3) = THE ARRIVAL TIME GR000380

C RAYTIME(N, 1 ) = Nl, NUMBER OF P TRIPS GR000390

C RAYTIME(N,2) = N2, NUMBER OF S TRIPS GR000400

C IF N2<0, IT IS A HEAD WAVE GR000410

C IF N1=- 1 , IT IS SURFACE WAVE GR000420

C IF N1=0 AND N2=0, IT IS A RAYLEIGH WAVE GR000430

C GR000440

C THE TIME OF ARRIVAL I M FORMAT I ON IS ALSO PRINTED ON LU 6 GR000450

C GR000460

C GR000470

C GR000480

DOUBLE PRECISION DISPL(NPOINT) GR000490

DIMENSION TA(101),CN(101 ,3) , RAYTIME ( 101 ,3) GR000500

C NRAY = NUMBER OF MAXIMUM RAYS EXPECTED GR000510

DATA NRAY/101/ GR000520

C GR000530

CALL INIT(ALPHA,XD, ZD, INDEX) GR000540
CALL T IMEARR I (ALPHA , XD , 2D , NRAY, RAYTIME , TA, CN

)

WRITE(6,*) 1 I N1 N2 TIME ARRIVAL' DO 500 1=1, NRAY WRITE (6,499) I , RAYTIME( I

, 1 ), RAYTIME ( I, 2) , RAYTIME ( I ,3) 499 FORMAT (2X,I5,3X,3F10.5)

500 CONTINUE

DC 1000 1=1 .NPOINT

TIME I = FLOAT ( I )*TDELTA

CALL GREENSUB (ALPHA , I NDEX , XD , ZD , NRAY , T I ME I , RAYT I ME , D I SPL ( I )

)

WRITE(6,*) 'TIME = ' ,T IMEI ,

' DISPL = ’,DISPL(I)

1000 CONTINUE

RETURN

END

SUBROUTINE INIT(ALPHA,XD, ZD, INDEX)

C£*********************************************

IMPLICIT DOUBLE COMPLEX (Y) IMPLICIT DOUBLE PRECISION (D) COMMON /BLKO/DXD,DT,DZ,DA

COMMON /BLK1/KASE ,M, L

, K, N1 ,N2,D,DTHETA

COMMON /BLK3/Y I

, YONE , YTWO, YTHREE , YFOUR , YE I GHT , YZERO

COMMON /BLK5/YR,Y1,Y2

DXD = DBLE(XD)

DZ = DBLE(ZD)

DA = DBLE(ALPHA)

IF (ZD.EQ.0.5) KASE = 1 IF (ZD . EQ. -0.5) KASE = 2

IF (INDEX.LT. 100) THEN

K = 0 M = INDEX/10

L = INDEX - M*10

ELSE IF ((INDEX. GT. 100). AND. (INDEX.LT. 334)) THEN

M = INDEX/100

L = (INDEX - M*100)/10

K = (INDEX - M*100 - L*10) ELSEWRITE(6,*) ' WRONG INDEX ' PAUSE

END IF

YI = DCMPLX(0.D0, 1 .DO) YONE = DCMPLX(1 .D0,0.D0) YTWO = DCMPLX(2.D0,0.D0)

YTHREE = DCMPLX(3.D0,0.D0)

YFOUR = DCMPLX(4.D0,0.D0)

YEIGHT = DCMPLX(8.D0, 0.D0) YZERO = DCMPLX(0.D0,0.D0)

CALL RAYRT(YR,Y1,Y2, ALPHA)

RETURN

END

C SUBROUT I NE GREENSUB( ALPHA, I NDEX , XD , ZD , NRAY , T , RAYT I ME , D I SPL

)

C

GR000550

GR000560

GR000570

GR000580

GR000590

GR000600

GR000610

GR000620

GR000630

GR000640

GR000650

GR000660

GR000670

GR000680

GR000690

GR000700

GR000710

GR000720

GR000730

GR000740

GR000750

GR000760

GR000770

GR000780

GR000790

GR000800

GR000810

GR000820

GR000830

GR000840

GR000850

GR000860

GR000870

GR000880

GR000890

GR000900

GR000910

GR000920

GR000930

GR000940

GR000950

GR000960

GR000970

GR000980

GR000990

GR001000

GR001010

GR001020

GR001030

GR001040

GR001050

GR001060

GR001070

GR001080

GR001090
£«•*»«•»*********************•****************************************

C SUBROUTINE TO COMPUTE THE GREENS FUNCTION OF A PLATE AT T C

C INPUTS:

C CONFIGURATION PARAMETERS:

C DETECTOR LOCATION:

C X = XD

C Y = 0

C Z = ZD, ZD =0.5 => TOP; ZD=-0.5 => BOTTOM

C SOURCE LOCATION:

C XO = 0

C YO = 0

C ZO = +0.5 (ON TOP OF THE PLATE)

C

C MATERIALS PARAMETERS:

C ALPHA = RATIO OF SHEAR VELOCITY TO LONGITUDINAL VELOCITY

C

C OTHER PARAMETERS:

C INDEX: SUBSCRIPTS OF GREENS FUNCTION i.e. 33, or 113

C NRAY : MAXIMUM NUMBER OF RAYS EXPECTED. -100

C RAYTIME: TIME OF ARRIVAL TABLE OF DIMENSION (NRAY, 3)

C

C T = NORMALIZED TIME = ACTUAL TIME *SHEAR WAVE SPEED / PLATE

C THICKNESS

C

C DISPLACEMENT WILL BE CORRECT ONLY IF T < RAYTIME (NRAY, 3)

C

C OUTPUT:

C DISPL = NORMALIZED DISPLACEMENT = ACTUAL DISPLACEMENT * PI

C * SHEAR MODULUS * PLATE THICKNESS / FORCE

C

C NOTE: SUBROUTINE INIT(ALPHA,XD, ZD, INDEX) AND

C SUBROUTINE TIMEARRK ) SHOULD BE CALLED FIRST IN C THE MAIN PROGRAM.

C IMPLICIT DOUBLE COMPLEX (Y) IMPLICIT DOUBLE PRECISION (D) DOUBLE PRECISION DINTG DIMENSION RAYTIME(NRAY,3)

COMMON /BLKO/DXD,DT,DZ,DA

COMMON /BLK1/KASE,M,L,K,N1,N2,D,DTHETA

COMMON /BLK2/YD,Y, YP, YQ,YR0,YRP, YRM,YDELTA,YPSQP1 , YQSQP1

,

1 YPSQM1 , YQSQM1 ,YSQRP1 , YSQRQ1 , YY,YA,YAA, YETA1 , YETA2

2 , YDPHI COMMON /BLK3/YI ,YONE,YTWO, YTHREE,YFOUR,YEIGHT, YZERO

COMMON /BLK5/YR, Y1 ,Y2

C EXTERNAL DINTG

C

C

IF ( (INDEX. EQ. 12). OR. (INDEX. EQ. 21). OR.

+ ( INDEX. EQ. 23). OR. ( INDEX. EQ. 32). OR.

+ ( INDEX. EQ. 121 ). OR. ( INDEX. EQ. 211). OR. + ( INDEX. EQ. 231 ). OR. ( INDEX. EQ. 321). OR. + ( INDEX. EQ. 123). OR. ( INDEX. EQ. 213). OR.

GR001100

GR001110

GR001120

GR001130

GR001140

GR001150

GR001160

GR001170

GR001180

GR001190

GR001200

GR001210

GR001220

GR001230

GR001240

GR001250

GR001260

GR001270

GR001280

GR001290

GR001300

GR001310

GR001320

GR001330

GR001340

GR001350

GR001360

GR001370

GR001380

GR001390

GR001400

GR001410

GR001420

GR001430

GR001440

GR001450

GR001460

GR001470

GR001480

GR001490

GR001500

GROO 1510 GR001520

GR001530

GR001540

GR001550

GR001560

GR001570

GR001580

GR001590

GR001600

GR001610

GR001620

GR001630

GR001640
c

c

c

c

c

c

c

c***

c

c***

c***

c

C8000 8000

( INDEX. EQ. 233). OR. ( INDEX. EO. 323). OR.

( INDEX. EQ. 132). OR. ( INDEX. EQ. 312). OR. (INDEX. EQ. 112). OR. (INDEX. EQ. 222). OR. (INDEX. EQ. 332) ) THEN DISPL=0.0D0

RETURN

ELSE

DZERO=O.DO

NMAX = NRAY; Maximum number of rays. The same as DIM RAYTIME

NMAX = NRAY

IF(T.GT. (RAY TIME (NRAY , 3 ) ) ) THEN

yp I TE C 6 * ) * ************************************

*

WRITE(6,*) 'ERROR IN SUBROUTINE GREEN' WRITE(6,*)'T > TARRIVAL(NRAY) , INCREASE NRAY.' WRITE(6,*) 'PROGRAM IS PAUSED'

yp I TE C 6 * ) * ***********************************

•

PAUSE 'GREENSUB'

END IF

DT=DBLE(T)

DSURF = 0.D0

IF (KASE.EQ.1) THEN

CALL SFWAVE(ALPHA,XD,T, INDEX, DSURF)

ELSE IF (KASE.EQ.2)THEN

DSURF = 0.D0

END IF DISPL=0.D0

DO 8000 N=1 ,NMAX

IF (T.LE.RAYTIME(N,3)) GO TO 8001

N1 = RAYTIME(N , 1

)

N2 = RAYTIME(N, 2)

I F ( ( N 1 . EQ.O) . AND . (N2. EQ . 0) ) GOTO 8000

I F ( ( N 1 . LT.O) .AND . (N2.EQ.0) ) GOTO 8000

IT IS RAYLEIGH WAVE ARRIVAL WHICH HAS BEEN TAKEN

CARE OF BY SUBROUTINE SFWAVE IF(N2.LT.O) THEN

IT IS A HEAD WAVE N2=-N2 DK2=D8LE(REAL(N2)

)

DTSARR = DSQRT (DXD*DXD + DK2*DK2)

IF (DT.GE. DTSARR) GO TO 8000 DHEAD = (DT/DA) -DK2*DSQRT ( 1 .0D0-DA*DA)/DA

DRAYI=DINTG( -DXD, -DHEAD) ELSE

IT IS A REGULAR RAY DRAYI=DINTG(DZERO, -DXD) END IF

DISPL = DISPL + DRAY I Y=YZERO ;set try root eq. 0 Y=YZERO

GR001650

GR001660

GR001670

GR001680

GR001690

GR001700

GR001710

GR001720

GR001730

GR001740

GR001750

GR001760

GR001770

GR001780

GR001790

GR001800

GR001810

GR001820

GR001830

GR001840

GR001850

GR001860

GR001870

GR 00 1880 GR001890

GR001900

GR001910

GR001920

GR001930

GR001940

GR001950

GR001960

GR001970

GR001980

GR001990

GR002000

GR002010

GR002020

GR002030

GR002040

GR002050

GR002060

GR002070

GR002080

GR002090

GR002100

GR0021 10 GR002120

GR002130

GR002140

GR002150

GR002160

GR002170

GR002180

GR002190
8001 D I SPL=0 I SPL/(6. 2831853071 795900)

D I SPL=DI SPL+DSURF

RETURN

END IF END

C DOUBLE PRECISION FUNCTION DINTG(DINIT,DEND)

C£*«*****************************************************************

c IMPLICIT DOUBLE PRECISION (D)

CC DOUBLE PRECISION EPS, R, E,W( 50,6) , FMIN, FMAX,

F

INTEGER NINT,NMAX,KF, I FLAG LOGICAL RST

C

C EXTERNAL F C EPS = 1 .0-6

RST = .FALSE. NMAX = 50 NINT = 1 W(1,1) = DINIT

W(2, 1 ) = DEND

C CALL DGLQ1 ( F,DINIT, DEND, EPS, R,E, NINT, RST,U, NMAX, FMIN, FMAX, KF

1 , 1 FLAG)

C

C DINTG=R

RETURN

END Q*******************************************************************

c DOUBLE PRECISION FUNCTION F(DS)

CQ*******************************************************************

c IMPLICIT DOUBLE COMPLEX <Y) IMPLICIT DOUBLE PRECISION (D) DOUBLE PRECISION ETA1 ,ETA2SQ,DD,AA DIMENSION YQD(3,3),YRSTAR(3,3), YQS(3,3) COMMON /BLK0/DXD,DT,DZ,DA

COMMON /BLK1/KASE,M,L,K,N1 ,N2,D,DTHETA

COMMON /BLK2/YD, Y, YP, YQ, YRO, YRP, YRM, YDELTA, YPSQP1 , YQSQP1

,

1 YPSQM1 , YQSQM1 , YSQRP1 , YSQRQ1 , YY, YA , YAA , YETA1 , YETA2

2 , YDPHI COMMON /BLK3/Y I , YONE , YTWO, YTHREE , YFOUR , YE I GHT , YZERO

C DK1=DBLE(REAL(N1 )

)

DK2=DBLE(REAL(N2)

)

D=DS

GR002200

GR002210

GR002220

GR002230

GR002240

GR002250

GR002260

GR002270

GR002280

GR002290

GR002300

GR002310

GR002320

GR002330

GR002340

GR002350

GR002360

GR002370

GR002380

GR002390

GR002400

GR002410

GR002420

GR002430

GR002440

GR002450

GR002460

GR002470

GR002480

GR002490

GR002500

GR002510

GR002520

GR002530

GR002540

GR002550

GR002560

GR002570

GR002580

GR002590

GR002600

GR002610

GR002620

GR002630

GR002640

GR002650

GR002660

GR002670

GR002680

GR002690

GR002700

GR002710

GR002720

GR002730

GR002740
c CALL PHIEQ0(DK1 , DK2, DA, DT,DXD,DS,Y, DERR, NROOT)

C WRITE(6,*)DK1 ,DK2,DS,Y, DERR, NROOT

IF (NROOT. EQ.O) THEN F=0.D0

C ELSE

CALL COMMON (DS) CALL RSTAR(YRSTAR, IDUMMY) CALL QD(YQD, IDUMMY) CALL QSU(YQS, IDUMMY)

C YINTG = YZERO

DO 200 1=1,2 DO 100 J=1 ,

2

YINTG=YINTG+YQD(M, I )*YRSTAR( I

, J )*YQS(L, J

)

100 CONTINUE

200 CONTINUE

C

IF (<N1. EQ.O). AND. (N2.GT.0))

1 YINTG=YINTG+YQD(M,3)*YQS(L ,3) C YINTG=YINTG*YDPHI/CDSQRT (DCMPLX(DXD*DXD-DS*DS)

)

F=DBLE(YINTG)

C END IF RETURN

END

C SUBROUTINE PHIEQ0(DK1 ,DK2,DA,DT,DXD, D,YROOT, DERR, NROOT)

C

C

C GIVEN DK1 ,DK2,DA,DT,D TO SOLVE PHI(YROOT)=0

IMPLICIT DOUBLE PRECISION (D) IMPLICIT DOUBLE COMPLEX (Y) DOUBLE PRECISION T,R,A,C1,C2

COMMON /PARM1/T,R,A,C1 ,C2

C YOLD=YROOT

IER=0 DYREAL=O.DO

YPHI=DCMPLX(0.D0,0.D0)

T = DT

R = D A = DA

Cl = OKI

C2 = DK2

EPS=1 . E - 1

2

IEND=100

C NROOT=1

IF (DK2.EQ.0.D0) THEN

C

GR002750

GR002760

GR002770

GR002780

GR002790

GR002800

GR002810

GR002820

GR002830

GR002840

GR002850

GR002860

GR002870

GR002880

GR002890

GR002900

GR002910

GR002920

GR002930

GR002940

GR002950

GR002960

GR002970

GR002980

GR002990

GR003000

GR003010

GR003020

GR003030

GR003040

GR003050

GR003060

GR003070

GR003080

GR003090

GR003100

GR003110

GR003120

GR003130

GR003140

GR003150

GR003160

GR003170

GR003180

GR003190

GR003200

GR003210

GR003220

GR003230

GR003240

GR003250

GR003260

GR003270

GR003280

GR003290
DSTSQ = (DT/DA)*(DT/DA) -DK1*DK1 DSTAR = DSQRT (DSTSQ)

IF ( D.LT.O. DO )DSTAR=- DSTAR

IF (D.EQ.O.DO) THEN DYIMAG=-DK1/(DSTAR*DA)

YROOT = DCMPLX(O.DO,DYIMAG)

ELSE IF (D.EQ. DSTAR) THEN DYREAL=DT/(DA*DA*DSTAR)

YROOT = DCMPLX(DYREAL f O.ODO)

ELSEDYREAL = D*DT/(DA*DA*DSTSQ)

DSMD=DSTSQ-D*D

IF(DSMD.LT.O.DO) WRITE(6,*) DT,D, DSHD, DSTSQ DYIMAG = * (DK1/(DSTSQ*DA) )*DSQRT (DSMD) YROOT = DCMPLX(DYREAL, DYIMAG)

END IF YOLD = YROOT

CALL DCRTNI (YROOT, YPHI , YDPHI , YOLD, EPS, I END, IER)

LSE IF (DK1.EQ.0.D0) THEN DTHARR = DA*DABS(DXD)+DK2*DSQRT(1.0D0-DA*DA)

DTSARR = DSQRT (DXD*DXD+DK2*DK2)

DSTSQ = (DT)*(DT) -DK2*DK2 DSTAR = DSQRT (DSTSQ)

IF (DT.GT. DTSARR) THEN

IF (D.LT.O. DO) DSTAR=-DSTAR

IF (D.EQ.O.DO) THEN DYIMAG=-DK2/DABS(DSTAR)

YROOT = DCMPLX(0. DO, DYIMAG)

ELSE IF (D.EQ. DSTAR) THEN DYREAL=OT/ (DSTAR) YROOT = DCMPLX(DYREAL, O.ODO)

ELSEDYREAL = D*DT/(DSTSQ)

DYIMAG = - (DK2/DSTSQ)*DSQRT (DSTSQ-D*D

)

YROOT = DCMPLX(DYREAL, DYIMAG)

END IF ELSE IF ((DT.LT. DTSARR). AND. (DT.GT. DTHARR)) THEN DHEAD = (DT/DA)-DK2*DSQRT(1.0DO-DA*DA)/DA

IF (D. LT. 0. DO )DHEAD=- DHEAD

IF (D.LT.O.DO)DSTAR=-DSTAR

IF (D.EQ.O.DO) THEN DYIMAG=-DK2/DSTAR

YROOT = DCMPLX(0. DO, DYIMAG)

ELSE IF (D.EQ. DSTAR) THEN DYREAL=DT/(DSTAR)

YROOT = DCMPLX(DYREAL, O.ODO) ELSE IF (DABS (D ) .GT. DABS (DSTAR ) ) THEN

IF (D.LT.O. DO) THEN DYREAL = (D*DT-DK2*DSQRT (D*D - DSTSQ ))/( DSTSQ)

ELSE IF (D.GT.O.DO) THEN

GR003300

GR003310

GR003320

GR003330

GR003340

GR003350

GR003360

GR003370

GR003380

GR003390

GR003400

GR003410

GR003420

GR003430

GR003440

GR003450

GR003460

GR003470

GR003480

GR003490

GR003500

GR003510

GR003520

GR003530

GR003540

GR003550

GR003560

GR003570

GR003580

GR003590

GR003600

GR003610

GR003620

GR003630

GR003640

GR003650

GR003660

GR003670

GR003680

GR003690

GR003700

GR003710

GR003720

GR003730

GR003740

GR003750

GR003760

GR003770

GR003780

GR003790

GR003800

GR003810

GR003820

GR003830

GR003840
DYREAL = (D*DT+DK2*DSQRT(D*D-DSTSQ))/(DSTSQ)

ENDIF DYIMAG = O.DO YROOT = DCMPLX<DYREAL, DYIMAG)

END IF END IF YOLD = YROOT

CALL DCRTNI (YROOT, YPHI , YDPHI , YOLD, EPS, IEND, IER)

C ELSE IF ((DK1. NE. O.DO). AND. (DK2.NE. O.DO)) THEN IF((D.EQ.O.DO).OR.(YOLD.EQ.DCMPLX(O.DO,O.DO))) THEN

CALL DEQO(DT, DA, DK1,DK2, YROOT, NOROOT)

ENDIF YOLD = YROOT

CALL DCRTNI (YROOT, YPH I, YDPHI , YOLD, EPS, IEND, IER) END IF

I F( IER.NE.O) NROOT=0

C

C TEST ROOT TO FIND ERROR

C DERR=CDABS(YPHI

)

C RETURN

END

C SUBROUTINE DEQO (DT , DA , DK 1

, DK2 , Y , NOROOT

)

C£*******************************************************************

C

C COMPUTE THE ROOT BY SOLVING QUARTIC EQ.

C IMPLICIT DOUBLE PRECISION (D) IMPLICIT DOUBLE COMPLEX (Y) DOUBLE PRECISION B DIMENSION DB(5) ,B(3),DROOT1 (2) ,DROOT2(2)

C YF( Y)=-DCMPLX(DT , O.DO)*Y-DCMPLX(DK1 ,O.DO)*

1 CDSQRT (Y*Y*DCMPLX(DA*DA,O.DO) -DCMPLX( 1 .DO, O.DO)

)

2 -DCMPLX(DK2,0.D0)*CDSQRT (Y*Y-DCMPLX( 1 .DO, O.DO)

)

DT2 = DT*DT

DK12 =DK1*DK1 DK22 = DK2*DK2

DAA=DA*DA

C DB( 1 )=( (DT2-DK12*DAA-DK22)**2 - 4.D0*DK12*DK22*DAA)

DB(3)=2.D0*((DT2-DK12*DAA-DK22)*(DK12+DK22)

1 +2.D0*DK12*DK22*(DAA+1.D0))

DB(5 )=(DK12+DK22)**2-4.D0*DK12*DK22

C

B ( 1 ) = DB( 1

)

B(2) = DB(3)

B C 3 ) = DB(5)

T0L=1.E-20

C CALL QUADRTIC EQ. TO FIND Y**2

GR003850

GR003860

GR003870

GR003880

GR003890

GR003900

GR003910

GR003920

GR003930

GR003940

GR003950

GR003960

GR003970

GR003980

GR003990

GR004000

GR004010

GR004020

GR004030

GR004040

GR004050

GR004060

GR004070

GR004080

GR004090

GR004100

GR0041 10 GR004120

GR004130

GR004140

GR004150

GR004160

GR004170

GR004180

GR004190

GR004200

GR004210

GR004220

GR004230

GR004240

GR004250

GR004260

GR004270

GR004280

GR004290

GR004300

GR004310

GR004320

GR004330

GR004340

GR004350

GR004360

GR004370

GR004380

GR004390
CALL QUADRT ( B , DROOT 1

, DROOT2 , TOL , NOROOT

)

C Yl=CDSORT (DCMPLX(DROOTl ( 1 ), DROOT 1 (2) )

)

Y2=CDSQRT (DCMPLX(DROOT2( 1 ) ,DROOT2(2) )

)

IF (DIMAG(YI).GT.O.DO) Y1=-Y1

IF (DIMAG(Y2) .GT.O.DO) Y2=-Y2

YF1=YF(Y1

)

YF2=YF(Y2)

DERR1=CDABS(YF(Y1 )) DERR2=CDABS(YF(Y2))

IF (DERR1.LT. DERR2) Y=Y1

IF (DERR1 .GE.DERR2) Y=Y2

C

IF (DIMAG(Y) .GE.O.DO) Y=DCONJG(Y)

C RETURN

END

C DOUBLE COMPLEX FUNCTION PHIO(Z)

C

c

C FUNCTION PHIO

C DOUBLE PRECISION T,R f A,C1,C2

DOUBLE COMPLEX Z, YAA, YONE COMMON /PARM1/T,R,A,C1 ,C2

C YONE=DCMPLX( 1 .DO,O.DO)

YAA=DCMPLX(A*A, 0.D0)

C PHIO=DCMPLX(R, 0.D0) -DCMPLX(T,0.D0)*Z

1 *DCMPLX(C1 , O.DO)*CDSQRT (YAA*Z*Z- YONE)

2 -DCMPLX(C2,O.DO)*CDSQRT(Z*Z-YONE)

C RETURN

END £****************************************-**************************

c DOUBLE COMPLEX FUNCTION DPHIO(Z)

C

c

C FUNCTION DPHIO

Cc DOUBLE PRECISION T,R,A f Cl,C2 DOUBLE COMPLEX Z f YAA, YONE COMMON /PARM1/T,R,A,C1,C2

C YONE=DCMPLX(1.DO,O.DO)

YAA=DCMPLX(A*A, 0.D0)

C

GR004400

GR004410

GR004420

GR004430

GR004440

GR004450

GR004460

GR004470

GR004480

GR004490

GR004500

GR004510

GR004520

GR004530

GR004540

GR004550

GR004560

GR004570

GR004580

GR004590

GR004600

GR004610

GR004620

GR004630

GR004640

GR004650

GR004660

GR004670

GR004680

GR004690

GR004700

GR004710

GR004720

GR004730

GR 004740

GR004750

GR004760

GR004770

GR004780

GR004790

GR004800

GR004810

GR004820

GR004830

GR004840

GR004850

GR004860

GR004870

GR004880

GR004890

GR004900

GR004910

GR004920

GR004930

DPHI0=-DCMPLX(T,0.D0) GR004940
1 - DCMPLX ( Cl , 0 . DO )*YAA*Z/CDSQRT ( YAA*Z*Z - YONE

)

2 -DCMPLX(C2,O.DO)*Z/CDSQRT(Z*Z-YONE)

C RETURN

END Q*******************************************************************

c SUBROUTINE COMMON(DS)

C

C

C TO COMPUTE PARAMETERS DEPENDING ON Y AND TO PUT IN

C A COMMON BLK2

C IMPLICIT DOUBLE COMPLEX (Y) IMPLICIT DOUBLE PRECISION (D) DOUBLE PRECISION ETA1 ,ETA2SQ,DD,AA

COMMON /BLKO/DXD,DT,DZ,DA

COMMON /BLK1/KASE,M,L,K,N1 ,N2,D,DTHETA

COMMON /BLK2/YD, Y, YP, YQ,YR0, YRP, YRM,YDELTA, YPSQP1 , YQSQP1

,

1 YPSQM1 , YQSQM1 , YSQRP1 .YSQRQl , YY,YA, YAA,YETA1 , YETA2

2 , YDPHI COMMON /BLK3/YI ,YONE, YTUO, YTHREE,YFOUR, YEIGHT,YZERO

C AA=DA*DA

DD=DS

D=DS DN1=DBLE(REAL(N1 )

)

DN2=DBLE(REAL(N2)

YZERO=DCMPLX(O.DO,O.DO)

YI=DCMPLX(0.D0, 1 .ODO)

YAA=DCMPLX(AA,O.DO)

YA=DCMPLX(DA,0.D0)

YONE=DCMPLX(1 .ODO.O.ODO)

YTUO=OCMPLX( 2.00,0. DO)

YFOUR=DCMPLX(4.0D0,0.D0)

C ETA1=DS/DXD

ETA2SQ=1.D0-ETA1*ETA1

YETA1 = DCMPLX(ETA1 ,0.D0)

YETA2S= DCMPLX (ETA2SQ,0. DO)

YETA2=CDSQRT(YETA2S)

C YD = DCMPLX(DD,0.D0)

YY=Y*Y YP=CDSQRT ( YAA*YY • YONE

)

YQ=CDSQRT(YY- YONE)

YPSQP1=YAA*YY

YQSQP1=YY

YQSQM1=YQ*YQ- YONE

YPSQM1 =YP*YP - YONE

C YSQRP1=Y*YA

YSQRQ1=Y

YDELTA=YQSQM1*YQSQM1+YFCUR*YP*YQ

GR004950

GR004960

GR004970

GR004980

GR004990

GR005000

GR005010

GR005020

GR005030

GR005040

GR005050

GR005060

GR005070

GR005080

GR005090

GR005100

GR0051 10 GR005120

GR005130

GR005140

GR005150

GR005160

GR005170

GR005180

GR005190

GR005200

GR005210

GR005220

GR005230

GR005240

GR005250

GR005260

GR005270

GR005280

GR005290

GR005300

GR005310

GR005320

GR005330

GR005340

GR005350

GR005360

GR005370

GR005380

GR005390

GR005400

GR005410

GR005420

GR005430

GR005440

GR005450

GR005460

GR005470

GR005480

GR005490
YM=(YQSQM1*YQSQM1 - YFOUR*YP*YQ)

YRO=YM/YDELTA

YRP=- (YF0UR*YP*YQSQM1 )/( YA*YDELTA)

YRM=( YR0*YR0- YONE)/YRP

YDPHI=Y0NE/(YD+(DCMPLX(DN1 , 0.D0)/YP)+(DCMPLX(DN2,0.D0)/YQ)

)

RETURN

END

C

SUBRCXIT I NE RSTARCYRSTAR, IDUMMY)

C

c

C TO COMPUTE RSTAR (JS,J2) EQ.6.39

C IMPLICIT DOUBLE COMPLEX (Y) IMPLICIT DOUBLE PRECISION (D) DIMENSION YRSTAR(3,3),KSUM(20)

COMMON /BLK0/DXD,DT,DZ,DA

COMMON /BLK1/KASE,M,L,K, N1 ,N2,D,DTHETA

COMMON /BLK2/YD, Y, YP, YQ, YRO, YRP,YRM,YDELTA, YPSQP1 , YQSQP1

,

1 YPSQM1 , YQSQM1 , YSQRP1 ,YSQRQ1 ,YY,YA, YAA, YETA1 , YETA2

2 , YDPHI COMMON /BLK3/YI , YONE, YTWO, YTHREE, YFOUR, YEIGHT,YZERO

EXTERNAL KDELTA

C

C IDUMMY=0

YRSTAR(3,3) = YONE YRSTAR(1 ,3) = YZERO

YRSTAR(2,3) = YZERO

YRSTAR(3, 1 ) = YZERO

YRSTAR(3,2) = YZERO

C DO 2000 JS=1 ,

2

DO 1000 J2=1 ,2

IF (J2.E0.1) N1=N1 -

1

IF ( J2.EQ.2) N2=N2-1

ISIGN = <-1)**(J2*(N1+N2)+N2)

CALL COEFF(N1,N2, JS, J2,KD f KU,KL,KSUM)

C WRITE(7,*) N1 , N2, JS, J2,KD, KU,XL,KSUM YSUM=YZERO

DO 100 NM=<L , KU YSUM=YSUM+DCMPLX<DBLE(REAL(KSUM(NM+1 ) ) )

, 0.D0)*

1 ( (YR0*YR0- YONE)/(YRO*YRO) )**NM

100 CONTINUE

KD2=KDELTA( 1,JS)-KDELTA(1,J2)

C YSUM = ( (YR0/YRP)**(KD2) )*YSUM

RSI GN=REAL( ISIGN) DRSIGN=RSIGN

DKD=DBLE(REALOCD))

YRSTAR( JS, J2)=DCMPLX(DRSI GN, 0.D0)*(YR0**(N1+N2)

)

GR005500

GR005510

GR005520

GR005530

GR005540

GR005550

GR005560

GR005570

GR005580

GR005590

GR005600

GR005610

GR005620

GR005630

GR005640

GR005650

GR005660

GR005670

GR005680

GR005690

GR005700

GR005710

GR005720

GR005730

GR005740

GR005750

GR005760

GR005770

GR005780

GR005790

GR005800

GR005810

GR005820

GR005830

GR005840

GR005850

GR005860

GR005870

GR005880

GR005890

GR005900

GR005910

GR005920

GR005930

GR005940

GR005950

GR005960

GR005970

GR005980

GR005990

GR006000

GR006010

GR006020

GR006030

GR006040
1 *(DCMPLX(DKD,O.DO)+YSUM)

IF (J2.EQ.1) N1=N1+1

IF (J2.EQ.2) N2=N2+1

C1000 CONTINUE

2000 CONTINUE

RETURN

END c*******************************************************************

c SUBROUTINE QD(YQD, IDUMMY)

C

C

C TO COMPUTE QD(

I

( J ) BASED ON EQ. 6.21

C IMPLICIT DOUBLE COMPLEX (Y) IMPLICIT DOUBLE PRECISION (D) COMMON /BLKO/DXD,DT,DZ,DA

COMMON /BLK1/KASE,M,L,K,N1 ,N2,D,DTHETA

COMMON /BLK2/YD, Y,YP,YQ,YRO,YRP, YRM, YDELTA, YPSQP1 ,YQSQP1

,

1 YPSQM1 , YQSQM1 , YSQRP1 , YSQRQ1 , YY , YA , YAA , YETA1 , YETA2

2 , YDPHI COMMON /BLK3/Y I

, YONE , YTUO, YTHREE , YFOUR , YE I GHT , YZERO

C DIMENSION YQD(3,3)

C 1DUMMY=0

DO 100 1=1,3 DO 90 J=1 ,3 YGD(I,J) = YZERO

90 CONTINUE

100 CONTINUE

C IF(M.EQ.I) THEN

YQD< 1, 1 )=YETA1*YFOUR*YP*YQ*YQSQP1

1 / (YSQRP1*YDELTA)

YQD( 1 ,2)=- YTW0*YETA1*YQ*YQSQM1*YSQRQ1

1 /(YDELTA)

YQD( 1,3) =- YTWO*YETA2

C ELSE IF (M.EQ.2)THEN

YQD(2, 1 )=YETA2*YF0UR*YP*YQ*YQSQP1

1 /( YSQRP1*YDELTA)

YQD(2,2)=-YTW0*YETA2*YQ*YQSQM1*YSQRQ1

1 /(YDELTA)

YQD(2,3)=YTW0*YETA1

C ELSE IF (M.EQ.3) THEN YQD(3 , 1 )=- YTU0*YP*YQSQM1*YQSQP1

1 /( YSQRP1*YDELTA)

YQD(3,2)=-YF0UR*YP*YQ*YSQRQ1

1 /(YDELTA)

YQD(3,3)=YZERO

GR006050

GR006060

GR006070

GR006080

GR006090

GR006100

GR006110

GR006120

GR006130

GR006140

GR006150

GR006160

GR006170

GR006180

GR006190

GR006200

GR006210

GR006220

GR006230

GR006240

GR006250

GR006260

GR006270

GR006280

GR006290

GR006300

GR006310

GR006320

GR006330

GR006340

GR006350

GR006360

GR006370

GR006380

GR006390

GR006400

GR006410

GR006420

GR006430

GR006440

GR006450

GR006460

GR006470

GR006480

GR006490

GR006500

GR006510

GR006520

GR006530

GR006540

GR006550

GR006560

GR006570

GR006580

GR006590
END IF RETURN

END

c*******************************************************************

C INTEGER FUNCTION KDELTAd , J) CC*******************************************************************

c

IF (I.EQ.J) THEN KDELTA=1

ELSEKDELTA=0

END IF RETURN

END Q*******************************************************************

c INTEGER FUNCTION IBINO(M,N)

C

IF (M.EQ.N)GO TO 11 IF ( (M. LT .0) .OR . (N. LT.O) ) GO TO 10 K=M-N

IF (1C) 10,11,12

10 IBINO = 0 RETURN

11 IBINO = 1 RETURN

12 IF (N.EQ.O)GO TO 11 IN = 1 IP = M

IQ = 1 DO 20 1=1,

<

IN = IN* IP/ IQ

IP = IP-1

IQ = IQ +1 20 CONTINUE

IBINO = IN RETURN

END

INTEGER FUNCTION KD1 (N1 ,N2, JS, J2)

Q*******************************************************************

C EXTERNAL KDELTA KD1=KDELTA<0,N1)*KDELTA(2,J2)+KDELTA<0,N2)*KDELTA(1,J2)

KD1=ICD1*KDELTA<JS,J2) RETURN

END c*******************************************************************

c

GR006600

GR006610

GR006620

GR006630

GR006640

GR006650

GR006660

GR006670

GR006680

GR006690

GR006700

GR006710

GR006720

GR006730

GR006740

GR006750

GR006760

GR006770

GR006780

GR006790

GR006800

GR006810

GR006820

GR006830

GR006840

GR006850

GR006860

GR006870

GR006880

GR006890

GR006900

GR006910

GR006920

GR006930

GR006940

GR006950

GR006960

GR006970

GR006980

GR006990

GR007000

GR007010

GR007020

GR007030

GR007040

GR007050

GR007060

GR007070

GR007080

GR007090

GR007100

GR0071 10 GR007120

GR007130

GR007140
INTEGER FUNCTION KUPPER(N1 ,N2, JS, J2)

C

EXTERNAL KDELTA

K=N2+KDELTA(2, J2)*KDELTA( 1

, JS) -KDELTA( 1

, J2)*KDELTA(2, JS)

IF (Nl.LE.K) THEN KUPPER=N1

ELSEKUPPER=K

END IF RETURN

END

INTEGER FUNCTION KLOWER(N1 , N2, JS , J2)

EXTERNAL KDELTA KLOUER=KDELTA(2, J2)+KDELTA( 1

, J2)*KDELTA( 1, JS) RETURN

END £*************************************************

C SUBROUTINE COEFF(Nl , N2, JS, J2,KD,KU,KL,KSUM)

C

C

C TO COMPUTE THE COEFFICIENTS FOR THE RSTAR TURM

C

C INPUTS:

C N1 = NUMBER OF P TRIPS

C N2 = NUMBER OF S TRIPS

C JS = FINAL TRIP TYPE, 1 OR 2

C J2 = INITIAL TRIP TYPE, 1 OR 2 C

C OUTPUT:

C KD = FIRST COEFF

C KU = UPPER LIMIT IN SUMMATION

C KL = LOWER LIMIT IN SUMMATION

C KSUM = COEFF. IN SUMMATION, MUST BE DIMENSIONED IN CALLING

C PROGRAM TO BE MAX(N1,N2)

C

C

c

C PROGRAM TESTCOEF

C DIMENSION KSUM(20) cc***

C 1 WRITE(5, 10)

C 10 FORMAT (2X, 1 TYPE K1 ,K2, JS, J2 •

)

C READ(5,*)I,J,JS,J2

C WRITE(6,30)I,J,JS,J2

C 30 FORMAT (2X, 1 I =• , 15, 1 J= , 15 ,

• JS=',I2,' J2=',I2) CC***

GR007150

GR007160

GR007170

GR007180

GR007190

GR007200

GR007210

GR007220

GR007230

GR007240

GR007250

GR007260

GR007270

GR007280

GR007290

GR007300

GR007310

GR007320

GR007330

GR007340

GR007350

GR007360

GR007370

GR007380

GR007390

GR007400

GR007410

GR007420

GR007430

GR007440

GR007450

GR007460

GR007470

GR007480

GR007490

GR007500

GR007510

GR007520

GR007530

GR007540

GR007550

GR007560

GR007570

GR007580

GR007590

GR007600

GR007610

GR007620

GR007630

GR007640

GR007650

GR 007660 GR007670

GR007680

GR007690
END IF RETURN

END

q* a*****************************************************************

c INTEGER FUNCTION KDELTA(I.J)

C

IF (I.EQ.J) THEN KDELTA=1

ELSEKDELTA=0

END IF RETURN

END

INTEGER FUNCTION IBINO(H,N)

IF (M.EQ.N)GO TO 11 IF ((M.LT.O).OR.(N.LT.O)) GO TO 10 K=M-N

IF (K) 10,11,12

10 IBINO = 0 RETURN

11 IBINO = 1 RETURN

12 IF (N.EQ.O)GO TO 11 IN = 1 IP = M

IQ = 1 DO 20 1=1 ,

K

IN = IN* IP/ IQ IP = IP-1

IQ = IQ +1 20 CONTINUE

IBINO = IN RETURN

END £*******************************

C INTEGER FUNCTION KD1 (N1 , N2, JS, J2)

C

EXTERNAL KDELTA KD1=KDELTA(0,N1)*KDELTA(2,J2)+KDELTA(0,N2)*KDELTA(1,J2)

KD1=KD1*KDELTA( JS, J2) RETURN

END £*******************************************-************************

C

GR006600

GR006610

GR006620

GR006630

GR006640

GR006650

GR006660

GR006670

GR006680

GR006690

GR006700

GR006710

GR006720

GR006730

GR006740

GR006750

GR006760

GR006770

GR006780

GR006790

GR006800

GR006310

GR006820

GR006830

GR006840

GR006850

GR006860

GR006870

GR006880

GR006890

GR006900

GR006910

GR006920

GR006930

GR006940

GR006950

GR006960

GR006970

GR006980

GR006990

GR007000

GR007010

GR007020

GR007030

GR007040

GR007050

GR007060

GR007070

GR007080

GR007090

GR007100

GR0071 10 GR007120

GR007130

GR007140
INTEGER FUNCTION KUPPER(N1 , N2, JS, J2)

CC*******************************************************************

C EXTERNAL KDELTA K=N2+KDELTA(2, J2)*KDELTA( 1

, JS)-KDELTA( 1, J2)*KDELTA(2, JS)

IF (Nl.LE.K) THEN KUPPER=N1

ELSEKUPPER=K

END IF RETURN

END

c INTEGER FUNCTION KLOWER(N1 ,N2, JS, J2)

CQ*******************************************************************

c EXTERNAL KDELTA KLOUER=KDELTA<2, J2)+KDELTA(1,J2)*KDELTA(1, JS) RETURN

END c*******************************************************************

c

c

SUBROUTINE COEFF(N1 ,N2, JS, J2,KD,KU, KL,KSUM)

c

c TO COMPUTE THE COEFFICIENTS FOR THE RSTAR TURM

cc INPUTS:

c N1 = NUMBER OF P TRIPS

c N2 = NUMBER OF S TRIPS

c JS = FINAL TRIP TYPE, 1 OR 2 c J2 = INITIAL TRIP TYPE, 1 OR 2 c

c OUTPUT:

c KD = FIRST COEFF

c KU = UPPER LIMIT IN SUMMATION

c KL = LOWER LIMIT IN SUMMATION

c KSUM = COEFF. IN SUMMATION, MUST BE DIMENSIONED IN CALLING

c PROGRAM TO BE MAX(N1,N2)

cr...

c

C PROGRAM TESTCOEF

C DIMENSION KSUM(20)

CC***

C 1 WRITE(5, 10)

C 10 FORMAT (2X, 1 TYPE K1 ,K2, JS, J2'

)

C READ (5,*) I

, J, JS, J2

C URITE(6,30)I , J, JS, J2

C 30 FORMAT (2X, 1 I =• , 15, J=* , 15, JS=',I2,' J2=',I2) CC***

GR007150

GR007160

GR007170

GR007180

GR007190

GR007200

GR007210

GR007220

GR007230

GR007240

GR007250

GR007260

GR007270

GR007280

GR007290

GR007300

GR007310

GR007320

GR007330

GR007340

GR007350

GR007360

GR007370

GR007380

GR007390

GR007400

GR007410

GR007420

GR007430

GR007440

GR007450

GR007460

GR007470

GR007480

GR007490

GR007500

GR007510

GR007520

GR007530

GR007540

GR007550

GR007560

GR007570

GR007580

GR007590

GR007600

GR007610

GR007620

GR007630

GR007640

GR007650

GR007660

GR007670

GR007680

GR007690
IF (KASE.EQ.1) THEN

C SOURCE AT TOP SURFACE

C DO 200 1=1,3

YQS( I , 1 )=YRRP( 1

, 1 )*YUP( I

, 1 )*YP1+YRRP( 1 ,2)*YUP( I

, 2)*YP2

1

- YUM( I

, 1 )*YP1

YQS( I ,2)=YRRP(2, 1 )*YUP( I

, 1 )*YP1

1 +YRRP(2,2)*YUP(I,2)*YP2 - YUM(I,2)*YP2

YQS(I,3)=(YRRP(3,3)*YUP(I,3) - YUM( I ,3))*YP2

200 CONTINUE

C

C ELSE IF (KASE.EQ.2) THEN

C SOURCE AT BOTTOM SURFACE

C DO 300 1=1,3

YQS( I

, 1 )=YRRM( 1, 1 )*YUM( I, 1 )*YP1+YRRM< 1 ,2)*YUM( I, 2)*YP2

1 - YUP( I, 1 )*YP1

YQS( I

, 2)=YRRM(2, 1 )*YUM( I

, 1 )*YP1

1 +YRRM(2,2)*YUM( I ,2)*YP2 - YUP<I,2)*YP2

YQS( I ,3)=(YRRM(3,3)*YUM( 1,3) - YUP(I,3) )*YP2 300 CONTINUE

C END IF RETURN

END

c *****************************************

C SUBROUTINE EPIDIS(A,T,U)

C

Q **************************

C

C

C THIS PROGRAM CALCULATES THE VERTICAL DISPLACEMENT

C AT THE EPICENTER OF A LARGE PLATE DUE TO AN UNIT

C STEP FUNCTION VERTICAL FORCE.

C DOUBLE PRECISION IS USED.

C

C INPUT PARAMETERS : C

C A = SHEAR WAVE SPEED / LONGITUDINAL WAVE SPEED

C e.g. ALUMINUM : A=0.495

C T = NORMALIZED TIME = ACTUAL TIME *SHEAR WAVE SPEED / PLATE

C THICKNESS

C e.g. SHEAR WAVE SPEED ~ 3 mm /micro second

C

C OUTPUT PARAMETERS : C

C U = NORMALIZED DISPLACEMENT = ACTUAL DISPLACEMENT * 3.14159

C * SHEAR MODULUS * PLATE THICKNESS / FORCE

C

C SUBROUTINES REQUIRED:

C

C QPARA(T,A>, AND FUNCTION BINO

C

GR008800

GR008S10

GR008820

GR008830

GR008840

GR008850

GR008860

GR008870

GR008880

GR008890

GR008900

GR008910

GR008920

GR008930

GR008940

GR008950

GR008960

GR008970

GR008980

GR008990

GR009000

GR009010

GR009020

GR009030

GR009040

GR009050

GR009060

GR009070

GR009080

GR009090

GR009100

GR0091 10 GR009120

GR009130

GR009140

GR009150

GR009160

GR009170

GR009180

GR009190

GR009200

GR009210

GR009220

GR009230

GR009240

GR009250

GR009260

GR009270

GR009280

GR009290

GR009300

GR009310

GR009320

GR009330

GR009340
c- GR009350

c GR009360

c THE PROGRAM MAY BE SET UP AS THE FOLLOWINGS: GR009370

c GR009380

c- PROGRAM EPICT GR009390

c- DOUBLE PRECISION A,DT,DELTAT,DU(1024) GR009400

c- A = 0.584D0 GR009410

c- NPT=1024 GR009420

c- DELTAT=0.00277D0 GR009430

c- DO 100 1=1 f NPT GR009440

c- DT=DELTAT*DBLE(REAL( I )

)

GR009450

c- CALL EPIDIS(A,DT,DU( I )

)

GR009460

c- 100 CONTINUE GR009470

c- WRITE(7) DU GR009480

c- STOP GR009490

c- END GR009500

c GR009510

GR009520

c GR009530

IMPLICIT DOUBLE PRECISION (Q) GR009540

DOUBLE PRECISION U,T,A,TA,DXM GR009550

INTEGER BINO GR009560

COMMON /BLOCK/ GN1 ,QN2,QO,QP,QM,QPP,QPS,QSP,QSS,QPHI GR009570

c GR009580

c NRAY = MAX. POSSIBLE NUMBER OF P RAYS GR009590

c NRAY2 = MAX. POSSIBLE NUMBER OF S RAYS GR009600

NRAY = 20 GR009610

NRAY2 = 12 GR009620

U = 0.0 GR009630

DO 300 NA = 1 , NRAY GR009640

DO 200 NB = 1.NRAY2 GR009650

N1 = NA1

GR009660

N2 = NB- GR009670

c N1+N2 MUST BE AN ODD NUMBER GR009680

IF (

(

FLOAT (N1+N2)/2. ) .EQ. FLOAT ( (N1+N2)/2) )GO TO 200 GR009690

TA = N2 + N1*A GR009700

IF (T.LE.TA) GO TO 200 GR009710

QN1 = N1 GR009720

QN2 = N2 GR009730

c GR009740

CALL QPARA(T,A) GR009750

c GR009760

NTERM = N2 GR009770

IF (N1.LT.N2) NTERM=N1 GR009780

NTERM=NTERM+1 GR009790

QSUM = 0.0 GR009800

DO 100 MT = 1, NTERM GR009810

MM = MT1

GR009820

QFPP = BINO(N1+1,MM+2)*BINO(N2+1,MM+1) GR009830

QFPS = BIN0(N2+1,MM+1)*BIN0(N1+1,MM+1) GR009840

QFSS = BINO(N2+1 ,MM+2)*BIN0(N1+1 ,MM+1

)

GR009850

MO = N1 + N2 - 1

- 2*MM GR009860

c GR009870

QSUM = (

(

- 1 . )**N2)* GR009880

1 ( (QFPP*QPP-QFSS*QSS)*(DXM(QO, MO)*DXM( (QM*QP*( - 1 . ) )

, MM)

)

GR009890
2 +QFPS*(QPS*DXM(QO, (MO+1 ) )*DXM( (

(

- 1 . )*QP) ,MM)*DXM(QM,MM)

3 -OSP*DXM(QO, (MO+1 ) )*OXM( ( (

- 1 . )*QM) ,MM)*DXM(QP, MM) )

)

4 +QSUM

C 100 CONTINUE

U = U + QSUM/QPHI

200 CONTINUE

300 CONTINUE

RETURN

END

C

Q *****************************************************************

c SUBROUTINE QPARA(T.A)

C

Q ***********************

C IMPLICIT DOUBLE PRECISION (Q) DOUBLE PRECISION T^B.AC.BAC.XX.X.SX.SAX

COMMON /BLOCK/ QN1 , QN2,QO, QP,QM,QPP,QPS,QSP,QSS f QPHI

C

C CALCULATE X C Q1=QN1*QN1

Q2=QN2*QN2

IF(QNI.EQ.O.DO) THEN XX=(T*T-Q2)/Q2

ELSE IF (QN2.EQ.0.D0) THEN XX=(T*T-Q1*A*A)/Q1

ELSE

B = ((Q1+Q2)*T*T -(Q1*A*A-Q2)*(Q1-Q2))/((Q1-Q2)**2)

AC= ( (Q1*A*A-Q2)**2+T*T*(T*T - 2

. *Q1 *A*A- 2. *Q2)

)

1 /( (Q1 -Q2)**2)

BAC =B*B -AC

IF (BAC. LE. 0.0) 3AC = 0.0

XX = B-DSQRT (BAC)

END IF X = DSQRT(XX)

SX = DSQRT ( 1 .+XX) SAX = DSQRT (A*A+XX)

C

C CALCULATE R'S, Q'S, AND QPHI

C QA= (1.+2.*XX)**2

QB= 4.*SAX*SX*XX

QO= (QA+QB)/(QA-QB)

QP= -4.*X*A*SX*(1.+2.*XX)/(QA-QB)

QM= 4.*X*SAX*(1.+2.*XX)/((QA-Q8)*A)

QPP= 2.*SAX*QA/((QA-QB)**2)

QPS= SAX/ (QA-QB)

QSP=-SAX/(QA-QB)

QSS= -2.*SAX*QB/ ( (QA-QB )**2)

QPHI= 2.*(QN1/SAX + QN2/SX) RETURN

END

GR009900

GR009910

GR009920

GR009930

GR009940

GR009950

GR009960

GR009970

GR009980

GR009990

GR010000

GR010010

GR010020

GR010030

GR010040

GR010050

GR010060

GR010070

GR010080

GR010090

GR010100

GR0101 10 GRO 101 20 GR010130

GR010140

GR010150

GR010160

GR010170

GR010180

GR010190

GR010200

GR010210

GR010220

GR010230

GR010240

GR010250

GR010260

GR010270

GR010280

GR010290

GR010300

GR010310

GR010320

GR010330

GR010340

GR010350

GR010360

GR010370

GR010380

GR010390

GR010400

GR010410

GR010420

GR010430

GRO 10440
Q *****************************************************************

c INTEGER FUNCTION BINO(J1,J2)

C

Q *****************************************************************

c M = J1-2

N = J2-2

IF ((M.EQ.N).OR.(N.EQ.O)) THEN

BINO=1 RETURN

ELSE IF ((M.LT.O).OR.(N.LT.O).OR.CM.LT.N)) THEN BINO = 0 RETURN

ELSEK=M-N

IN = 1 IP = M

IQ = 1 DO 20 1=1,

K

IN = IN* IP/ IQ

IP = IP-1

IQ = IQ +1 20 CONTINUE

BINO = IN RETURN ENDIF

END AT**** C FUNCTION DXM(D,M)

C

C

C THE FUNCTION DXM(D,M) IS

C THE SAME AS D**M. BOTH D AND DXM MUST BE DECLARED IN C THE MAIN PROGRAM AS DOUBLE PRECISION VARIABLES.

C DOUBLE PRECISION DXM,D

IF (M.EQ.O) DXM=1.0D0

IF (M.GT.O) GO TO 11 IF (M.LT.O) GO TO 22 RETURN

C

11 DXM=1.0DO

DO 10 1=1,

M

DXM=DXM*D

10 CONTINUE

RETURN

C22 DXM=1.0D0

N=-M DO 20 1 = 1,

N

DXM=DXM/D

GR010450

GR010460

GRO 10470 GR010480

GR010490

GR010500

GR010510

GR010520

GR010530

GR010540

GR010550

GR010560

GR010570

GR010580

GR010590

GR010600

GR010610

GR010620

GR010630

GR010640

GR010650

GRO 10660

GRO 10670 GR0106S0

GRO 10690 GR010700

GR010710

GR010720

GR010730

GR010740

GR010750

GR010760

GR010770

GR010780

GRO 10790 GR010800

GR010810

GR010820

GR010830

GR010840

GR010850

GRO 10860 GR010870

GR010880

GR010890

GR010900

GR010910

GR010920

GR010930

GR010940

GR010950

GRO 10960

GRO 10970

GRO 10980 GR010990
20 CONTINUE GR011000

RETURN GR011010

C GR011020

END GR011030

GR011040

C GR011060

SUBROUTINE SFWAVE(ALPHA,XR, TTIME, INDEX, DISPL) GRO 11060

C GR011070

c** •Jrt***************************** GR011080

c SFWAVE GR011090

c SUBROUTINE TO COMPUTE THE SURFACE DISPLACEMENT GRO 111 00

c INPUTS: GR011110

c ALPHA=RATIO OF SHEAR WAVE SPEED TO LONG. WAVE SPEED GR011120

c XR =D I STANCE BETWEEN SOURCE AND DETECTOR GR011130

c =ACTUAL DIST./H GR011140

c TTIME =NOND IMENSIONALIZED TTIME GR011150

c =ACTUAL TIME*SHEAR WAVE SPEED / H GR011160

c INDEX=SU8SCRIPT OF THE GREEN'S FUNCTION GR011170

c e.g. 11,22,33,12,13 or 131,113, etc. GR011180

c GR011190

c OUTPUT

:

GR011200

c D I SPL=NOND IMENSIONALIZED DISPLACEMENT GR011210

c =ACTUAL DISPL.* PI * SHEAR MODULUS * H / FORCE GR011220

c NOTE: GR011230

c THE PARAMETERS IN COMMON BLKS 0,1, & 5 SHOULD BE ENTERED GR011240

c BEFORE THE FIRST CALL OF THIS SUBROUTINE, GR011250

c I. E. CALL INIT(ALPHA,XR, ZD, INDEX) FIRST. GR011260

c GR011270

c SUBROUITNES REQUIRED: GR011280

c GR011290

c DSFINT,DQR,QU,COMMON,RAYRT & INTEGRATION ROUTINES GR011300

c GR011310

c REMARK: THE COMPUTATION IS DONE IN Y-PLANE, 4tn QUARD. ON A GR011320

c THREE LEGGED U-PATH TO REPLACE THE INTEGRATION FROM GR011330

c -V TO V. GR011340

c GR011350

c GR011360

c GR011370

IMPLICIT DOUBLE PRECISION (D) GR011380

IMPLICIT DOUBLE COMPLEX (Y) GR011390

COMMON /BLKO/DXD,DT,DZ,DA GR011400

COMMON /BLK1/KASE,M,L,K,N1,N2,D,DTHETA GR011410

COMMON /BLK5/YR, Y1 , Y2 GR011420

c GR011430

N1=0 GR011440

N2=0 GR011450

DYR=OBLE(YR) GR011460

DT=DBLE(TTIME) GR011470 DZERO=O.ODO GR011480

DMONEH=-0.5D0 GR011490

DSQRYR=DSQRT (DYR) GR011500

RWSPD=DSQRYR GR011510

c GR011520

DPF=1.2D0/DBLE(ALPHA) GR011530

DV=DXD/DT GR011540
IF (TTIME .LE . (ALPHA*XR) ) THEN DISPL=O.ODO

RETURN

ELSE DISPL=O.DO

D I SPL=D I SPL+DSF I NG(DZERO, DMONEH , 1

)

2 +DSFING(DV,DPF,2)

3 +DSFING(DMONEH,DZERO,3)

END IF END

c DOUBLE PRECISION FUNCTION DSFINGCDINIT,DEND, I PATH)

C£************1Ht******************************************************

C IMPLICIT DOUBLE PRECISION (D)

C DOUBLE PRECISION EPS,R, E ,U(50,6) , FMIN, FMAX,

F

INTEGER NINT,NMAX,KF, I FLAG LOGICAL RST

C

C EXTERNAL DQR1 ,DQR2,DQR3

C EPS = 1.D-8

RST = .FALSE. NMAX =50

NINT = 1 W(1,1) = DINIT W(2,1) = DEND

IF ( IPATH.EQ. 1 ) THEN

CALL DGLQ1 (DQR1,DINIT , DEND, EPS, R,E, NINT, RST, W,NMAX, FMIN, FMAX,

+ KF, I FLAG) ELSE IF (IPATH.EQ. 2)THEN

CALL DGLQ1(DQR2,DINIT,DEND,EPS,R,E, NINT, RST,W,NMAX, FMIN, FMAX,

+ KF, I FLAG) ELSE IF (IPATH.EQ. 3) THEN

CALL DGLQ1 (DQR3, DINIT, DEND, EPS, R,E, NINT, RST, W, NMAX, FMIN, FMAX,

+ KF, I FLAG) ELSE PAUSE 'DSFING: IPATH .NE. 1,2, OR 3' END IF C D2PEI=8.D0*DASIN(1 .DO)

DSFING=-R/D2PEI

RETURN

END C*******************************************************************

c DOUBLE PRECISION FUNCTION DQRI(DL)

CC  ** ***** ************** ******************************************

IMPLICIT DOUBLE PRECISION (D) EXTERNAL DQR

GR011550

GR011560

GR011570

GR011580

GR011590

GR011600

GR011610

GRO 11620 GR011630

GR011640

GR011650

GR011660

GR011670

GRO 11680 GR011690

GR011700

GR011710

GR011720

GR011730

GRO 11740 GR011750

GR011760

GR011770

GRO 11780

GRO 11790 GR011800

GR011810

GR011820

GRO 11830 GR011840

GR011850

GR011860

GR011870

GRO 11880 GR011890

GR011900

GRO 1 1910 GR011920

GR011930

GRO 11940 GR011950

GR011960

GR011970

GRO 11980

GRO 11990 GR012000

GR012010

GR012020

GR012030

GR012040

GR012050

GR012060

GR012070

GR012080

GR012090
DQR1 =DQR(DL, 1

)

RETURN

END

(2 *******************************************************************

C DOUBLE PRECISION FUNCTION DQR2(DL)

C£*******************************************************************

IMPLICIT DOUBLE PRECISION (D) EXTERNAL DQR

DQR2=DQR(DL, 2) RETURN

END

C DOUBLE PRECISION FUNCTION DQR3CDL)

C

IMPLICIT DOUBLE PRECISION (D) EXTERNAL DQR

DQR3=DQR(DL,3)

•RETURN END £***********************

C DOUBLE PRECISION FUNCTION DQR(DL, IPATH)

C

C IMPLICIT DOUBLE COMPLEX (Y) IMPLICIT DOUBLE PRECISION (D) DIMENSION YQD(3,3) , YUP(3,3) ,YQK(3) , YFEE(3) COMMON /BLKO/DXD,DT,DZ,DA

COMMON /BLKl/KASE,M,L,K f N1,N2,D,DTHETA

COMMON /BLK2/YD, Y, YP, YQ, YR0,YRP, YRM, YDELTA,YPSQP1 , YQSQP1

,

1 YPSQM1 , YQSQM1 , YSQRP1 , YSQRQ1 , YY, YA, YAA, YETA1 , YETA2

2 , YDPHI COMMON /BLK3/YI , YONE, YTWO, YTHREE, YFOUR, YEIGHT, YZERO

Cc DV=DXD/DT

IF (IPATH. EQ.1) THEN

Y= DCMPLX( DV,DL)

ELSE IF (IPATH. EQ. 2) THEN Y= DCMPLX(DL, -0.5D0) ELSE IF (IPATH. EQ. 3) THEN

DPF=1 .2D0/DA Y= DCMPLX(DPF,DL)

ELSEPAUSE 'DRQ: IPATH ?' END IF CALL SFCOM(DL)

CALL QD(YQD, IDUMMY) CALL QU(YUP, IDUMMY)

C

GR012100

GR012110

GR012120

GR012130

GR012140

GR012150

GR012160

GR012170

GR012180

GR012190

GR012200

GR012210

GR012220

GR012230

GR012240

GR012250

GR012260

GR012270

GR012280

GR012290

GR012300

GR012310

GR012320

GR012330

GR012340

GR012350

GR012360

GRO 12370 GR012380

GR012390

GR012400

GR012410

GR012420

GR012430

GR012440

GR012450

GRO 12460 GR012470

GR012480

GR012490

GR012500

GRO 1 25 1

0

GR012520

GR012530

GR012540

GR012550

GR012560

GR012570

GR012580

GR012590

GR012600

GR012610

GR012620

GR012630

GR012640
YFEEC 1 )=YAA/ (YP*YD) GR012650

YFEE(2)=Y0NE/(YQ*YD) GRO 12660 YFEE(3)=Y0NE/(YQ*YD) GR012670

c GRO 12680

YQR = YZERO GR012690

c GR012700

IF (K.EQ.O) THEN GR012710

YQK(1 )=YONE GR012720

YQK(2)=YONE GR012730

YQK(3)=YONE GRO 12740

ELSE IF (K.EQ.1) THEN GR012750

YQK(1 )=YETA1/Y GRO 12760 YQK(2)=YETA1/Y GR012770

YQK(3)=YETA1/Y GR012780

ELSE IF OC.EQ.2) THEN GRO 12790

YQK( 1 )=YETA2/Y GR012800

YQK(2)=YETA2/Y GR012810

YQK(3)=YETA2/Y GR012820

ELSE IF (K.EQ.3) THEN GR012830

YQK( 1 )=- YP/Y GRO 12840 YQK(2)=- YQ/Y GR012850

YQK(3)=- YQ/Y GR012860

ELSE GRO 12870 PAUSE ' FUUCTION DOR: WRONG INDEX, < .NE. 0,1,2, OR 3' GRO 12880

END IF GR012890

DO 100 1=1,3 GR012900

YQR=YQR+YQD(M, I )*YUP(L, I )*YQK( I ) *Y FEE C I

)

GR012910

100 CONTINUE GRO 12920

c GR012930

IF ( ( I PAT H.EQ.1).OR.(I PATH . EQ . 3 ) ) THEN GR012940

YQR=YQR*YI/CDSQRT (DCMPLX(DV*DV,0.D0) - YY) GR012950

ELSE GRO 12960 YQR=YQR/CDSQRT (DCMPLX(DV*DV, 0.D0) - YY) GRO 12970

END IF GRO 12980 DQR=DBLE(YQR) GR012990

c GR013000

RETURN GR013010

END GR013020

^*****************************Hr****** GR013030

c GR013040 SUBROUTINE SFCOM(DL) GR013050

c GR013060

0***r***********************1***************************************** GR013070

c GR013080

c TO COMPUTE PARAMETERS DEPENDING ON Y AND TO PUT IN GR013090

c A COMMON BLK2 GR013100

c GR0131 10 IMPLICIT DOUBLE COMPLEX (Y) GR013120

IMPLICIT DOUBLE PRECISION (D) GR013130

DOUBLE PRECISION ETA1 ,ETA2SQ,DD,AA GR013140

COMMON /BLK0/DXD,DT,DZ,DA GR013150

COMMON /BLK1/KASE,M,L,K,N1 ,N2,D,DTHETA GR013160

COMMON /BLK2/YD, Y, YP, YQ, YRO, YRP, YRM, YDELTA, YPSQP1 , YQSQP1

,

GR013170

1 YPSQM1 , YQSQM1 , YSQRP1 , YSQRQ1 , YY , YA , YAA, YETA1 , YETA2 GR013180

2 , YDPHI GR013190
COMMON /BLK3/YI , YONE, YTWO, YTHREE, Y FOUR, YE I GHT, YZERO

C AA=DA*DA

DV=DXD/DT

DN1=0BLE(REAL(N1 )) DN2=DBLE(REAL(N2)>

YZERO=OCMPLX( 0.00,0. DO)

YI=DCMPLX(0.D0, 1.0D0)

YAA=DCMPLX(AA, 0.D0)

YA=DCMPLX(DA,0.D0)

YONE=DCMPLX(1. 000,0.000)

YTWO=DCMPLX(2.DO,O.DO)

YFOUR=DCMPLX(4. 000,0. DO)

C YD=- Y*DCMPLX(DT,0.D0)

YETA1 = -Y/OCMPLX(DV,O.DO)

YETA2=CDSQRT(YONE-YETA1*YETA1

)

IF (OBLE(Y) .LT.O.DO) YETA2=-YETA2

YY=Y*Y YP=CDSQRT ( YAA*YY - YONE

)

YQ=CDSQRT<YY-YONE)

YPSQP1=YAA*YY

YQSQP1=YY

YQSQM1=YQ*YQ- YONE YPSQM1=YP*YP- YONE

C YSQRP1=Y*YA

YSQRQ1=Y

YDELTA=YQSQM1*YQSQM1+YF0UR*YP*YQ

YM=(YQSQM1*YQSQM1 - YFOUR*YP*YQ)

YRO=YM/YDELTA

YRP=* (YF0UR*YP*YQSQM1 )/(YA*YDELTA)

YRM=(YRO*YRO-YONE)/YRP

C RETURN

END C*******************************************************************

SUBROUTINE QU(YUP, IDUMMY)

£*******************************

C

C TO COMPUTE QU(I,J) BASED ON EQ. 6.4

C IMPLICIT DOUBLE COMPLEX (Y) COMMON /BLKO/DXD , DT , DZ , DA COMMON /BLK1/KASE,M,L,K,N1,N2,D,DTHETA

COMMON /BLK2/YD,Y, YP,YQ,YR0, YRP, YRM, YDELTA, YPSQP1 , YQSQP1

,

1 YPSQM1 , YQSQM1 , YSQRP1 , YSQRQ1 , YY , YA , YAA , YETA1 , YETA2

2 , YDPHI COMMON /BLK3/Y I

, YONE , YTWO , YTHREE , YFOUR , YE I GHT , YZERO

C DIMENSION YUP(3,3)

C IDUMMY=0

GR013200

GR013210

GR013220

GRO 13230 GR013240

GR013250

GR013260

GR013270

GR013280

GR013290

GR013300

GR013310

GR013320

GR013330

GR013340

GR013350

GR013360

GR013370

GR013330

GR013390

GR013400

GR013410

GR013420

GR013430

GR013440

GR013450

GR013460

GRO 13470 GR013480

GR013490

GR013500

GR013510

GR013520

GR013530

GR013540

GR013550

GR013560

GR013570

GR013580

GR013590

GR013600

GR013610

GR013620

GR013630

GR013640

GR013650

GR013660

GR013670

GRO 13680 GR013690

GR013700

GR013710

GR013720

GR013730

GR013740
YUP(1 , 1 )=YETA1/YSQRP1

YUP(1,2)=-YETA1*YQ/YSQRQ1

YUP< 1

, 3)=- YETA2

YUP(2, 1 )=YETA2/YSQRP1

YUP(2,2)=- YETA2*YQ/YSQRQ1

YUP(2,3)=YETA1

YUP(3 f 1 )=-YP/YSQRP1

YUP(3 f 2)=- Y0NE/YSQRQ1

YUP(3,3)=YZERO

C RETURN

END C*******************************************************************

C SUBROUTINE RAYRT(YR,YI 1 ,YI2, ALPHA)

CQ*******************************************************************

c

C TO SOLVE RAYLEIGH EQUITION

CC INPUT: ALPHA = SHEAR WAVE SPEED / LONG. WAVE SPEED

C

C OUTPUT: YR, YI 1

, YI2 ARE THE ROOTS TO

C

C Y*Y*Y-8*Y*Y+8*(3-2*AA)*Y- 16*( 1 - AA) = 0 C DOUBLE COMPLEX YR,YI1,YI2

DOUBLE PRECISION A,AA,RT1 ,RT2,RT3 DIMENSION A(4),RT1(2),RT2(2),RT3(2)

C

I F( (ALPHA. LE.O.). OR. (ALPHA. GE. 1.0)) THEN WRITE(6,*) 'RAYRT : WRONG VALUE OF ALPHA' WRITE(6,*) ALPHA

WRITE (6,*) 'RAYRT : COMPUTATION PAUSED' PAUSE

END IF C AA=ALPHA*ALPHA

TOL= 1.E-37

A( 1 ) = 1.0 A(2) = -8.0DO

A(3) = 8.*(3. -2.*AA)

A(4) = -16.*(1 . -AA)

C CALL CUB I C( A , RT 1

, RT2 , RT3 , TOL , NOROOT

)

C

C IF (ALPHA.LT. 0.567) THEN

C

C ONLY ONE REAL ROOT

C

C ELSE

C

C ALL THREE ROOTS ARE REAL

C

IF ( (RT1 ( 1 ) . LT . RT2( 1 ) ) .AND . (RT1 ( 1 ) . LT.RT3( 1 ) ) )THEN

GR013750

GRO 13760 GR013770

GR013780

GRO 13790 GR013800

GR013810

GR013820

GR013830

GR013840

GR013850

GR013860

GRO 13870 GR013880

GR013890

GR013900

GR013910

GR013920

GR013930

GR013940

GR013950

GR013960

GR013970

GR013980

GR013990

GR014000

GR014010

GR014020

GR014030

GR014040

GR014050

GR014060

GR014070

GR014080

GR014090

GR014100

GR0141 10 GRO 141 20 GR014130

GR014140

GR014150

GR014160

GR014170

GR014180

GR014190

GR014200

GR014210

GR014220

GR014230

GR014240

GR014250

GR014260

GR014270

GR014280

GR014290
YR=DCMPLX(RT1 (1),RT1(2))

YI1=DCMPLX(RT2(1 ),RT2(2))

YI2=DCMPLX(RT3(1),RT3(2))

ELSEIF ((RT2(1).LT.RT1(1)) .AND . (RT2( 1 ) . LT .RT3( 1 ) ) ) THEN YI1=DCMPLX(RT1(1),RT1(2))

YR =DCMPLX(RT2(1 ),RT2(2))

YI2=DCMPLX(RT3(1 ),RT3(2)) ELSEIF ( (RT3( l).LT.RTI(D) .AND . (RT3( 1 ) . LT .RT2( 1 ) ) )THEN

YI2=DCMPLX(RT1 ( 1),RT1(2))

YI 1=0CMPLX(RT2( 1 ) ,RT2(2)

)

YR =DCMPLX(RT3( 1 ), RT3(2)

)

END IF RETURN

END

C SUBROUT I NE T I MEARR I ( A , R , ZD , NRAY , RCN , TA , CN

)

Cc*******************************************************************

c

C 9/27/82

C

C DIMENSION Z1(4),Z2(4), TA(NRAY),CN(NRAY,3),RCN(NRAY,3)

DOUBLE COMPLEX YR,YI1,YI2

I TER= 1 00 NRAY=NRAY

1=0 LM=1

C Z = ZD+0.5 ; ZD refers to origin at the center of the plate

Z = ZD+0.5

Z1(1)=0.0

Z2(1 )=Z

Z1 (2)=Z

Z2(2)=0.0

Z1(3)=-Z

Z2(3)=0.0

Z1(4)=0.0

Z2(4)=-Z

C

IF (Z.EQ.1.0) Z2(1)=0.0

OREVEN=FLOAT(INT(Z))*0.5

DO 100 J=1 ,21 DO 200 K=1 ,12 JJ=J-1

KK=K1

IF ( (K.EQ.I).AND.(J.EQ.I) ) THEN IF(Z.EQ.I.O) THEN

C THREE SURFACE WAVE ARRIVALS:

1 = 1+1 CN( I

, 1 )=- 1

.

CN( I ,2)=0.

CN( I ,3)=R*A

TA( I )=CN( I ,3)

1 = 1+1

GR014300

GR014310

GR014320

GR014330

GR014340

GR014350

GR014360

GR014370

GR014380

GRO 14390 GR014400

GR014410

GR014420

GRO 14430 GR014440

GR014450

GR014460

GR014470

GR014480

GR014490

GR014500

GR014510

GR014520

GR014530

GR014540

GR014550

GR014560

GR014570

GR014580

GR014590

GR014600

GR014610

GR014620

GR014630

GR014640

GR014650

GR014660

GR014670

GR014680

GR014690

GR014700

GR014710

GR014720

GR014730

GR014740

GR014750

GRO 14760 GR014770

GR014730

GR014790

GR014800

GR014810

GR014820

GR014830

GR014840
CN( I

, 1 )=- 1

.

CN(I,2)=-1.

CN( I ,3)=R

TA( I )=CN( 1,3) CALL RAYRT (YR,YI1,YI2,A)

1 = 1 + 1 CN( I

, 1 )=0.

CN( I ,2)=0.

CN < 1 ,3)=R/DSQRT(DBLE(YR))

TA( I )=CN( I ,3) END IF GOTO 200

ENDIF

IF ( (Z.EQ.O.O).OR.(Z.EQ.I.O) ) LM=4

IF ( ((FLOAT( JJ+KK))/2.0-OREVEN) .EQ. FLOAT ( ( J J+KK)/2) )

1 GO TO 200 DO 300 LL=1,4,LM

C1=JJ+Z1(LL)

C2=KK+Z2(LL)

IF ( (C1.LT.0.0) .OR. (C2.LT. 0.0) ) GO TO 300

1 = 1+1

IF (I.GT.NRAY) GO TO 99

CALL TARRV (A,C1 ,C2,R, I TER , ERROR, T)

IF (ERROR. GT.1.E-6) WRITE (6,603) I ,C1 ,C2, ERROR

CN( I

, 1 )=C1

CN( I ,2)=C2

CN( I ,3)=T

TA( I )=T

IF (CI.GT.O.)GO TO 300 APJ=A*C2/SQRT ( 1

.

- A*A)

IF(R.LE.APJ) GO TO 300 C2=-C2

CALL TARRV (A,C1 ,C2,R, I TER, ERROR, T)

1 = 1+1

IF (I.GT.NRAY) GO TO 99

CN( 1 , 1 )=C1

CN( I ,2)=C2

CN( I ,3)=T

TA( I )=T 300 CONTINUE

200 CONTINUE

100 CONTINUE

C 99 CONTINUE

CALL RSORT(CN,TA,RCN,NRAY,3,0)

C603 FORMAT(2X, 'ERROR IS LARGER THAN 1.E-6 • , I3,3F10.5) RETURN

END

C SUBROUTINE TARRV(A,C1 ,C2,R, ITER,ERR,TARR)

CC*******************************************************************

C

GR014850

GRO 14860

GRO 14870 GR014880

GRO 14890 GR014900

GR014910

GR014920

GR014930

GR014940

GR014950

GR014960

GR014970

GR014980

GR014990

GR015000

GR015010

GR015020

GR015030

GR015040

GR015050

GR015060

GR015070

GR015080

GR015090

GRO 1 5 1 00 GR0151 10 GR015120

GR015130

GRO 151 40 GR015150

GR015160

GR015170

GR015180

GR015190

GR015200

GR015210

GR015220

GR015230

GR015240

GR015250

GR015260

GR015270

GR015280

GR015290

GR015300

GR015310

GR015320

GR015330

GR015340

GR015350

GR015360

GR015370

GRO 15380 GR015390
IMPLICIT DOUBLE PRECISION (D) GR015400

c GR015410

NUM = 0 GR015420

DTI =0.1570796326794901 GR015430

DT2=0.0D0 GR015440

DA=D8LE(A) GR015450 DAA=DA*DA GR015460

DN1=DBLE(C1

)

GR015470

DN2=DBLE(C2) GR015480

DR=DBLE(R) GR015490

IF (C1.EQ.0.0) GO TO 2 GR015500

IF (C2.EQ.0.0) GO TO 3 GR015510

c GR015520

1 CONTINUE GR015530

NUM = NUM + 1 GR015540

DT=(DT1+DT2)/2.0D0 GR015550

DAASS=DAA*DSIN(DT)*DSIN(DT) GR015560

DAASS=1 .ODO-DAASS GR015570

DCS2=DSQRT (DAASS) GR015580

DRT = DN1*DTAN(DT) + DN2*DA*DSIN(DT)/DCS2 GR015590

c URITE (6,20) NUM,DT,DRT,C1 ,C2 GR015600

IF (NUM. GE. ITER) GO TO 10 GR015610

IF (DABS(DRT-DR) .LE.5.D- 11) GO TO 10 GR015620

IF (DRT.LE.DR) DT2=DT GR015630

IF (DRT.GT.DR) DT1=DT GR015640

GO TO 1 GR015650

c GR015660

c 20 FORMAT ('NUM= , 13, ' T= * ,D17. 10, 1 RT= ' ,017. 10,2(2X, F4. 0)

)

GR015670

10 ITER=NUM GR015680

DTARR=DA*DN1/DCOS(DT) + DN2/DCS2 GR015690

DERR=DABS(DR-DRT) GR015700

ERR=DERR GR015710

TARR=DTARR GR015720

RETURN GR015730

c GR015740

2 IF (C2.LT.0.0) GO TO 4 GR015750

DTARR=DSQRT (DR*DR + DN2*DN2) GR015760

TARR=DTARR GR015770

ERR=0.0 GR015780

RETURN GR015790

c GR015800

3 DTARR=DSQRT (DR*DR+DN1*DN1 )*DA GR015810

TARR=DTARR GR015820

ERR=0.0 GR015830

RETURN GR015840

c GR015850

4 DTARR=DR*DA - DN2*DSQRT ( 1 .DO-DAA) GR015860

c NOTE HERE DN2 < 0 GR015870

TARR=DTARR GR015880

ERR=0.0 GR015890

RETURN GR015900

c GR015910

END GR015920

Q******************************************************************* GR015930

c GR015940
SUBROUT I NE CUB I C( A , RT1 , RT2 , RT3 , TOL , NOROOT

)

C£**********************************

C

C THIS SUBROUTINE FINDS THE ROOTS OF A CUBIC EQUITION

C

C INPUTS:

C A COEFFICIENTS, A(4)

C

C OUTPUTS:

C RT1 ,RT2,RT3 THREE ROOTS

C RT1 (1 )=REAL PART OF RT1 , ETC.

C

C THE ROUTINE CALLS LINCNG,AND QUADRT

C DOUBLE PRECISION A,B,C,RT1 ,RT2,RT3,ZERO f X,Y,ZZ,SQT3,RC27f ANG DIMENSION A(4) ,RT1 (2) ,RT2(2) , RT3(2) , B(2) ,C(4)

C ZERO=TOL/10.

SQT3=DSQRT (3.0DO) RC27=1. 0/27.0

C

IF (DABS(A( 1 )

) -ZERO) 7,7,12

7 CALL QUADRT (A(2),RT1 ,RT2, TOL, NOROOT) RETURN

C

12 NOROOT =3 IT=0

C DO 1 1=2,4

A( I )=A( I )/A( 1

)

1

CONTINUE

A( 1 )=1 .0 C

IF (DABS(A(2) ) .LE .ZERO) GOTO 2 NDA=3

B ( 1 )=1 .0

B ( 2 )= - (A(2)/3.0)

CALL LINCNG(A,NDA,B,C)

I T= 1 GO TO 4

2 DO 3 1=1,4

C( I )=A( I

)

3 CONTINUE

C

4 X = C(4)*C(4)*0.25+(C(3)**3)*RC27

IF (X.LT.0.0) GO TO 9

X = DSQRT(X)

Y = -(C(4)*0.5)

1=1 RTI(I) =Y+X

C

5 N = 0

I F C RT 1 ( I ) .LT.0.0) N=1

IF (DABS(RTI(D).LE.ZERO) GO TO 6

GR015950

GR015960

GR015970

GR015980

GR015990

GR016000

GR016010

GR016020

GR016030

GR016040

GR016050

GRO 16060 GR016070

GRO 16080 GR016090

GR016100

GR016110

GR016120

GR016130

GR016140

GR016150

GR016160

GR016170

GR016180

GR016190

GR016200

GR016210

GR016220

GR016230

GR016240

GR016250

GR016260

GR016270

GR016280

GR016290

GR016300

GR016310

GR016320

GR016330

GR016340

GR016350

GRO 16360 GR016370

GR016380

GR016390

GR016400

GR016410

GR016420

GR016430

GR016440

GR016450

GRO 16460 GR016470

GR016480

GR016490
RT 1 ( I ) = DABS(RT1(I))**(1./3.) GR016500

IF (N.EQ.1) RT1(I) = -RTI(I) GR016510

6 IF (I.EQ.2) GOTO 8 GR016520

I = 2 GR016530

RT 1 C I ) = Y-X GR016540

GO TO 5 GR016550

c GR016560

8 RT2(2) = ( (RT1 ( 1 )

- RT1(2))*0.5)*SQT3 GR016570

RT1 ( 1 ) = RT1 (1 ) + RT1 (2) GR016580

RT 1 < 2 > = 0.0 GR016590

RT2(1 ) = - RT 1 ( 1 )*0.5 GR016600

RT3( 1 ) = RT2(1

)

GR016610

RT3C2) = -RT2(2) GR016620

GO TO 11 GR016630

c GRO 16640 9 ZZ = DABS(C(3>) GR016650

X = -(C(4)*0.5)/DSQRT((ZZ**3)*RC27 ) GR016660

ANG = DACOS(X) GR016670

Y = 2.0*(DSQRT (ZZ/3.0)

)

GR016680

ANG=ANG/3.0 GRO 16690

RT1 ( 1 ) = Y*DCOS(ANG) GR016700

RT2C1) = Y*DCOS(ANG + 2. 094395 1024D0) GR016710

RT3(1 ) = Y*DCOS(ANG + 4. 1887902047D0) GR016720

RT1 (2) = 0.0 GR016730

RT2(2) = 0.0 GR016740

RT3(2) = 0.0 GR016750

c GRO 16760

11 IF (IT.EQ.1) THEN GRO 16770

RT1 ( 1 ) = RT 1(1) + B(2) GRO 16780

RT2( 1 ) = RT2(1 ) + B(2) GR016790

RT3(1 ) = RT3( 1 ) + B(2) GR016800

ELSE GR016810

END IF GR016820

RETURN GR016830

END GR016840

£***•****************************

GR016850

c GR016860 SUBROUTINE LINCNGCA, NDA,B,C) GR016870

c GR016880

£****************************************************************** GR016890

c GR016900

c TO MAKE A LINEAR CHANGE OF VARIABLES GR016910

c IN A GIVEN POLYNOMIAL GR016920

c GR016930

DOUBLE PRECISION A,B,C,BIPWR GR016940

DIMENSION A(1 ), B(2), C(1

)

GR016950

c GR016960

BIPWR = 1.0D0 GR016970

NZ = NDA +1 GR016980

c GR016990

DO 4 1=1 ,NZ GR017000

C(I) = A(I) GR017010

4 CONTINUE GR017020

c GR017030

6 DO 8 1=2, NZ GR017040
C(I) = C(I) + CC 1-1 )*B(2) GR017050

CONTINUE GROI 7060

GROI 7070

C(NZ) = C(NZ)*BIPWR GROI 7080

NZ = NZ - 1 GROI 7090

BIPWR = BIPWR*B( 1

)

GR017100

IF (NZ.GT.1) GO TO 6 GR0171 10 C(1) = C( 1 )*BIPWR GR017120

RETURN GR017130

END GR017140

£************-**********+******************************************** GR01 71 50

C GR017160 SUBROUTINE DGLQ1(F,A,B,EPS f RjE^INTjRST.W, NMAX,FMIN,FMAX,KF, I FLAG)GR017170

C GR017180

q ** ************************************************

grOI 71 90

C GROI 7200

C GR017210

C *** ALL REAL VARIABLES ARE TYPED DOUBLE PRECISION GR017220

C 9/20/82 NH GR017230

C GROI 7240

C GROI 7250

C DGLQ1 IS A SUBROUTINE FOR THE AUTOMATIC EVALUATION GROI 7260

C OF DEFINITE INTEGRALS OF A USER DEFINED FUNCTION GROI 7270

C OF ONE VARIABLE PROVIDING FLEXIBLE USAGE. GR017280

C GROI 7290

C FOR AN EASY TO USE VERSION SEE SUBROUTINE DGLQ. GR017300

C GR017310

C CAPABILITIES OF DGLQ1 (IN ADDITION TO THOSE OF DGLQ) GR017320

C INCLUDE: GR017330

C ABILITY TO RESTART A CALCULATION TO GREATER GR017340

C ACCURACY WITHOUT PENALTY... GR017350

c ABILITY TO SPECIFY AN INITIAL PARTITION OF GR017360

c THE INTEGRATION INTERVAL... GROI 7370

c ABILITY TO INCREASE THE WORK SPACE TO HANDLE GROI 7380

c MORE DIFFICULT PROBLEMS... GR017390

c OUTPUT OF LARGEST/SMALLEST INTEGRAND VALUE FOR GROI 7400

c APPLICATIONS SUCH AS GRAPHING... GR017410

c GROI 7420

c A R G UMENTS IN THE CALL SEQUENCE GROI 7430

c GROI 7440

c F (INPUT) THE NAME OF YOUR INTEGRAND FUNCTION. GR017450

c THIS NAME MUST APPEAR IN AN EXTERNAL STATEMENT GROI 7460

c IN ANY PROGRAM WHICH CALLS DGLQ1. GROI 7470

c YOU MUST WRITE F IN THE FORM GROI 7480

c FUNCTION F(X) GROI 7490

c F=(EVALUATE INTEGRAND AT THE POINT X) GROI 7500

c RETURN GR017510

c END GROI 7520

c A GR017530

c B (INPUT) ENDPOINTS OF INTEGRATION INTERVAL GR017540

c EPS (INPUT) ACCURACY TO WHICH THE INTEGRAL IS TO BE CALCULAGR017550

c DGLQ1 WILL TRY TO ACHIEVE RELATIVE ACCURACY , GR017560

c E.G. SET EPS=.01 FOR 2 DIGITS, .001 FOR 3, ETGR017570

c R (OUTPUT) THE ESTIMATE OF THE INTEGRAL GROI 7580

c E (OUTPUT) THE ESTIMATE OF THE ABSOLUTE ERROR IN R. GR017590
NIMT (INPUT GR017600

OUTPUT) GR017610

AS AN OUTPUT QUANTITY, NINT GIVES THE GR017620

NUMBER OF SUBINTERVALS IN THE FINAL GR017630

PARTITION OF [A,B] . GR017640

AS AN INPUT QUANTITY, NINT MUST 8E SET TO GR017650

THE NUMBER OF SUBINTERVALS IN THE INITIAL GR017660

PARTITION OF [A,B] . FOR MOST PROBLEMS GR017670

THIS IS JUST 1, THE INTERVAL [A,B] ITSELF. GR017680

NINT IS USEFUL IF YOU WOULD LIKE TO HELP GRO 17690

DGLQ1 LOCATE A DIFFICULT SPOT ON [A,B] . GR017700

IN THIS REGARD NINT IS USED ALONG GR017710

WITH THE ARRAY W (SEE BELOW). IF YOU SET GR017720

N I NT= 1 IT IS NOT NECESSARY TO BE CONCERNED GR017730

WITH W, EXCEPT THAT IT MUST BE DIMENSIONED. . .GR01 7740

AS AN EXAMPLE OF MORE GENERAL APPLICATIONS, GR017750

IF [A,B] = [0, 1] BUT THE INTEGRAND JUMPS AT 0.3GR017760

IT WOULD BE WISE TO SET NINT=2 AND THEN SET GR017770

W(1,1)=0.0 (LEFT ENDPOINT) GR017780

W(2, 1 )=0.3 (SINGULAR POINT) GR017790

W(3, 1 )=1 .0 (RIGHT ENDPOINT) GR017800

IF YOU SET NINT GREATER THAN 1, BE SURE TO GR017810

CHECK THAT YOU HAVE ALSO SET GR017820

W( 1, 1 )=A AND W(NINT+1 , 1 )=B GR017830

RST (INPUT) A LOGICAL VARIABLE (E.G. TRUE OR FALSE) GR017840

SET RST=. FALSE. FOR INITIAL CALL TO DGLQ1 GR017850

SET RST=.TRUE. FOR A SUBSEQUENT CALL, GR017860

E.G. ONE FOR WHICH MORE ACCURACY IS GR017870

DESIRED (SMALLER EPS). A RESTART ONLY GR017880

MAKES SENSE IF THE PRECEDING CALL RETURNED GR017890

WITH A VALUE OF I FLAG (SEE BELOW) LESS THAN GR017900

ON A RESTART YOU MAY NOT CHANGE THE VALUES 0GR017910

ANY ARGUMENTS IN THE CALL SEQUENCE, EXCEPT EGR017920

W(NMAX,6) W IS AN ARRAY USED FOR SCRATCH STORAGE BY DGLQ1 .GR017930

YOU MUST INCLUDE A DIMENSION STATEMENT IN GR017940

YOUR CALLING PROGRAM TO ALLOCATE THIS STORAGE. GR017950

THIS SHOULD BE OF THE FORM GR017960

DIMENSION W(NMAX,6) GR017970

WHERE NMAX IS AN INTEGER. AN ADEQUATE VALUE OF GR017980

NMAX IS 50. GR017990

NMAX (INPUT) AN INTEGER EQUAL TO THE FIRST SUBSCRIPT IN THE GR018000

DIMENSION STATEMENT FOR THE ARRAY W. THIS IS GR018010 ALSO EQUAL TO THE MAXIMUM NUMBER OF SUB INTERVAL GR0 18020 PERMITTED IN THE INTERNAL PARTITION OF [A,B]. GR018030

A VALUE OF 50 IS AMPLE FOR MOST PROBLEMS. GR018040

FMIN GR018050

FMAX (OUTPUT) THE SMALLEST AND LARGEST VALUES OF THE I NTEGRANGR018060

WHICH OCCURRED DURING THE CALCULATION. THE GR018070

ACTUAL INTEGRAND RANGE ON [A,B] MAY, OF COURSGRO 18080

BE GREATER BUT PROBABLY NOT BY MORE THAN 10X.GR018090

KF (OUTPUT) THE ACTUAL NUMBER OF INTEGRAND EVALUATIONS USE0GR018100

BY DGLQ1 TO APPROXIMATE THIS INTEGRAL. KF GR018110

WILL ALWAYS BE AT LEAST 30. GR018120

I FLAG (OUTPUT) TERMINATION FLAG. . .POSSIBLE VALUES ARE GR018130

0 NORMAL COMPLETION, E SATISFIES GR018140
c E<EPS AND E<EPS*DABS(R) GR018150

c 1 NORMAL COMPLETION, E SATISFIES GR018160

c E<EPS, BUT NOT RELATIVE ERROR REQUEST GR018170

c 2 NORMAL COMPLETION, E SATISFIES GR018180

c E<EPS*DA8S(R) GR018190

c 3 NORMAL COMPLETION BUT EPS WAS TOO SMALL TO GR018200

c SATISFY ABSOLUTE OR RELATIVE ERROR REQUEST. GR018210

c 4 ABORTED CALCULATION BECAUSE OF SERIOUS ROUNDING GR018220

c ERROR. PROBABLY E AND R ARE CONSISTENT. GR018230

c 5 ABORTED CALCULATION BECAUSE OF INSUFFICIENT STORAGR018240

c R AND E ARE CONSISTENT. PERHAPS INCREASING NMAGR018250

c WILL PRODUCE BETTER RESULTS. GRO 18260

c 6 ABORTED CALCULATION BECAUSE OF SERIOUS DIFFICULTIGR018270

c MEETING YOUR ERROR REQUEST. GR018280

c 7 ABORTED CALCULATION BECAUSE EITHER EPS, NINT OR NGR018290

c HAS BEEN SET TO AN ILLEGAL VALUE. GR018300

c 8 ABORTED CALCULATION BECAUSE YOU SET NINT>1 BUT FOGR018310

c TO SET W( 1, 1 )=A AND W(NINT+1 , 1 )=B GR018320

c GR018330

c TYPICAL PROBLEM SET UP GR018340

c GR018350

c DIMENSION U(50,6) GR018360

c LOGICAL RST GRO 18370

c EXTERNAL F GR018380

c A=0.0 GR018390

c B=1 .0 GR018400

c W( 1, 1 )=A GR018410

c W(2,1 )=.3 [SET INTERNAL PARTITION POINT AT .3] GR018420

c W(3, 1 )=B GR018430

c N I NT =2 [INITIAL PARTITION HAS 2 INTERVALS] GR018440

c RST=. FALSE. GR018450

c EPS=.001 GR018460

c NMAX=50 GR018470

c GR018480

c CALL DGLQ1(F,A, B,EPS,R, E,NINT,RST,W,NMAX, FMIN, FMAX,KF, I FLAG) GR018490

c GR018500

c 1 A,8,EPS ,R,E,NINT,FMIN, FMAX,KF, I FLAG GR018510

c IF(EPS.EQ. .0001 .OR. I FLAG.GE.3)STOP GR018520

c RST=.TRUE GR018530

c EPS=.0001 [ASK FOR ANOTHER DIGIT] GR018540

c GO TO 1 GR018550

c END GRO 18560

c FUNCTION F(X) GR018570

c IF(X.LT. .3) GR018580

c 1 THEN GR018590

c F=X**(0.2D0)*ALOG(X) GR018600

c ELSE GR018610

c F=DSIN(X) GR018620

c END IF GR018630

c RETURN GR018640

c END GR018650

c GR018660

c END OF DO CUMENTATION GRO 18670

c GR018680

INTEGER C GR018690
DOUBLE PRECISION A , B , E , EB , EPMACH , EPS , FMAX , FMAXL , FMAXR , FM I N , FM I NL

1 , FMINR, FMN, FMX,R f RAB, RABS,RAV,SIGN, T,TE

2 ,TE1 , TE2,TR, TR1 , TR2, UFLOW, W,XM,

F

DIMENSION W(NMAX,6)

EXTERNAL F

GR018700

GR018710

GRO 18720 GR018730

GR018740

LOGICAL RST.DEBUG GR018750

DATA SIGN /-I. DO/ GRO 18760

C THE FOLLOWING DATA ARE FOR PERKIN-ELMER COMPUTERS GRO 18770 DATA EPMACH , UFL0W/Z341 0000000000000 , Z001 000000000000/

C THE FOLLOWING DATA ARE FOR VAX COMPUTERS

GR018780

GRO 18790

C DATA EPMACH f UFLOW/Z0000250000000000 , Z0000008000000000/

C

GR018800

GR018810

C EPMACH = Y

1 341 0000000000000' GR018820

C UFLOW = Y'00100000' GR018830

C EPMACH = DBLE(16**(- 13))

C UFLOW = DBLE(16**(-65))

I F(A.EQ.B) THEN

R=0.

GRO 18840 GR018850

GR018860

GRO 18870

E=0. GR018880

NINT=0 GRO 18890

I FLAG=0 GR018900

KF=1 GR018910

FMIN=F(A) GR018920

FMAX=FMIN GRO 18930 GO TO 20 GR018940

END IF GR018950

IF(RST) THEN IFdFLAG.LT. 3) GO TO 15 GO TO 20

GRO 18960

GRO 18970 GR018980

END IF GR018990

KF=0 GR019000

IFCEPS .LE. 0. .OR. NMAX .LE. 1 .OR. NINT .LE. 0) THEN

I FLAG=7

GR019010

GR019020

GO TO 20 GR019030

END IF GR019040

IF(NINT.EQ.I) THEN

C 1 THEN

GR019050

GR019060

W( 1, 1 )=A

W(2,2)=B

W(1 ,5)=A

W(1 ,6)=B

W(2,5)=A

W<2,6)=8

W(1 ,2)=A+(B-A)/2.D0

W(2,1)=A+(B-A)/2.D0

NINT=2

GR019070

GR019080

GR019090

GR019100

GR0191 10 GR019120

GR019130

GR019140

GR019150

ELSE GR019160

I F(W( 1

, 1 ) .NE .A .AND. WCNINT+1 , 1 ) .NE.B) THEN IFLAG=8

GR019170

GR019180

GO TO 20 GR019190

END IF GR019200

W(1 ,5)=A DO 89 1=1 , NINT

W< I ,2)=W< 1+1 , 1

)

W( I ,5)=W( 1,1)

GR019210

GR019220

GR019230

GR019240
W( I ,6)=W(I ,2) 89 CONTINUE

GR019250

GR019260

END IF GR 01 9270

C GR019280

DEBUG=. FALSE. GR019290

I FLAG = 0 GR019300

IROFF=0 GR019310

RABS=0.0 GR019320

DO 3 1=1 ,NINT GR019330

CALL GL15T(F,W(I f 1),W(I f 2) f W(I,5),W(I,6),

1 W(I ,4),W(I ,3),RAB,RAV, FMN, FMX)

I F(DEBUG) WRITE (7 ,*) 1 INITIALIZE (W( I

, J) ,J=1 ,6)

KF=KF+15

GRO 19340 GR019350

GR019360

GRO 19370 IF(I.EQ.I)

1 THEN

GR019380

GRO 19390

R=W( 1,4)

E=W( I ,3) RABS=RABS+RAB

GR019400

GR019410

GR019420

FMIN=FMN GR019430

FMAX=FMX GR019440

ELSE GR019450

R=R+W( I ,4) E=E+U( I ,3) RABS=RABS+RAB

GR019460

GR019470

GRO 19480 FMAX=DMAX1 ( FMAX , FMX

)

FMIN=DMIN1 ( FMIN, FMN) END IF

GR019490

GR019500

GR019510

3 CONTINUE GR019520

DO 10 I=NINT+1,NMAX GR019530

W( I ,3) = 0. 10 CONTINUE

GR019540

GR019550

15 CONTINUE GR019560

C GR019570

C MAIN SUBPROGRAM LOOP GR019580

C GR019590

I F ( 1 00 . DO*EPMACH*RABS . GE .DABS( R ) .AND. E.LT.EPSJGO TO 20 EB=DMAX1(100.D0*UFLOW,DMAXl (EPS,50.D0*EPMACH)*DABS(R)

)

IF(E.LE.EB) GO TO 20

GR019600

GR019610

GRO 19620

IF (NINT.LT.NMAX)

1 THEN

GR019630

GR019640

NINT = NINT+1 GR019650

C = NINT GR019660

ELSE GRO 19670 C=0 GR019680

16 I F(C.EQ.NMAX) THEN

I FLAG=5

GRO 19690 GR019700

GO TO 20 GR019710

END IF GRO 19720

C=C+1 GR019730

I F(W(C,3) .GT.O.DO) GO TO 16 END IF

GR019740

GR019750

C LOC=ISAMAX(NINT,W(1 , 3) , 1

)

CALL LSE(W,NMAX,NINT,LOC)

I F(DEBUG)WRITE(7,200)LOC, (W(LOC, I),I=1,6),R,E

XM = W(LOC, 1 )+(W(LOC,2)-U(LOC, 1 ))/2.

GR019760

GR019770

GR019780

GRO 19790
IF ( (DMAX1 (DA3S(U(L0C , 1 ) ) ,DABS(W(LOC, 2) ) ) ) .GT. GR019800

1 ( (

1

.D0+100.D0*EPMACH)*(DABS(XM)+0. 1D+04*UFLOW))) GR019810

2 THEN GR019820

CALL GL15T ( F,W(LOC, 1 ) ,XM, W(L0C,5) ,W(L0C,6)

,

GR019830

1 TR1 ,TE1 , RAB,RAV, FMINL, FMAXL) GR019840

KF=KF+15 GR019850

IF (TE1.LT. (EB*(XM-U(LOC, 1))/(B-A))) GR019860

A TE1 = TE1*SIGN GRO 19870 CALL GL15T(F,XM l W(LOC,2),W(LCC t 5),U(LOC,6), GR019880

1 TR2,TE2,RA8,RAV, FMINR, FMAXR) GR019890

KF=KF+15 GR019900

FMIN=DMIN1(FMIN, FMINL, FMINR) GR019910

FMAX=DMAX1 (FMAX, FMAXL, FMAXR) GR019920

IF (TE2.LT. (EB*(W(L0C,2)-XM)/(B-A))) TE2=TE2*SIGN GR019930

TE = DABS(W(L0C,3)) GR019940

TR = U(LOC, 4) GR019950

U(C, 3) = TE2 GR019960

W(C,4) = TR2 GR019970

W(C,1) = XM GR019980

W(C,2) = W(L0C,2) GR019990

W(C,5) = W(L0C,5) GR020000

W(C,6) = W(LOC, 6) GR020010

W(L0C,3) = TE1 GR020020

W(L0C,4) = TR1 GR020030

W(L0C,2) = XM GR020040

I F(DEBUG)WRITE(7, 200)C, (U(C,K),K=1,6) GR020050

I F( DEBUG) WRITE (7, 200) LOC, (W(L0C,K) ,K=1 ,6) GR020060

E = E-TE+(DABS(TE1)+DABS(TE2)) GR020070

R = R-TR+(TR1+TR2) GR020080

I F (DEBUG)WRITE(7 ,*)NINT,R,E GR020090

I F(DABS(DABS(TE1 )+DABS(TE2)-TE) .LT. .001D0*TE) THEN GR020100

IROFF=IROFF+1 GR0201 10

I F( IROFF.GE .10) THEN GR020120

I FLAG=4 GR020130

GO TO 20 GR020140

END IF GR020150

END IF GR020160

ELSE GR020170

IF (EB.GT.U(LOC,3)) GR020180

1 THEN GR020190

W(LOC,3) = 0. GR020200

ELSE GR020210

IFLAG=6 GR020220

GO TO 20 GR020230

END IF GR020240

END IF GR020250

GO TO 15 GR020260

GR020270

ALL EXITS FROM HERE GR020280

GR020290

20 CONTINUE GR020300

I F ( I FLAG . GE . 4 )RETURN GR020310

I FLAG=3 GR020320

T=EPS*DABS(R) GR020330

I F(E.GT.EPS .AND. E.GT.T)RETURN GR020340
I FLAG=2

I F(E.GT.EPS .AND. E.LT.T)RETURN

IFLAG=1 IF(E.LT.EPS .AND. E.GT.T)RETURN

IFLAG=0 RETURN

200 FORMAT ( I4,8(E1 1 .3)) END

SUBROUTINE LSE (WORK,NMAX,NI , LOC)

C

C THIS SUBPROGRAM FINDS THE CELL IN THE WORK AREA

THAT HAS THE LARGEST ABSOLUTE ERROR AND RETURNS THE LOCATION OF THAT CELL.

DOUBLE PRECISION ERROR, WORK DIMENSION WORK(NMAX,6)

INITIALIZE VARIABLES

ERROR = DABS(WORK( 1 ,3)

)

LOC = 1

MAIN SUBPROGRAM LOOP

DO 20 1=1, NI

IF (DABS(WORK(I,3)).GT.ERROR) THEN ERROR = DABS(WORK( 1,3)) LOC = I END IF 20 CONTINUE

RETURN

END

SUBROUTINE GL15T(F,A,B,XL,XR,R,AE,RA,

1 RASC, FMIN, FMAX)

C***AUTHORS

C

C

ROBERT PIESSENS AND ELISE DE DONCKER

APPL. MATH. AND PROGR. DIV.- K.U. LEUVEN DAVID KAHANER, NBS WASHINGTON

PURPOSE

TO COMPUTE I = INTEGRAL OF G(X) OVER (A, B)

,

WITH ERROR ESTIMATE

J = INTEGRAL OF DABS(G) OVER (A,B)

PARAMETERS

ON ENTRY

F - DOUBLE PRECISION

GR020350

GR020360

GR020370

GR020380

GR020390

GR020400

GR020410

GR020420

GR020430

GR020440

GR020450

GR020460

GR020470

GR020480

GR020490

GR020500

GR020510

GR020520

GR020530

GR020540

GR020550

GR020560

GR020570

GR020580

GR020590

GR020600

GR020610

GR020620

GR020630

GR020640

GR020650

GR020660

GR020670

GR020680

GR020690

GR020700

GR020710

GR020720

GR020730

GR020740

GR020750

GR020760

GR020770

GR020780

GR020790

GR020800

GR020810

GR020820

GR020830

GR020840

GR020850

GR020860

GR020870

GR020880

GR020890
c FUNCTION SUBPROGRAM DEFINING THE INTEGRAND GR020900

c FUNCTION F(X) . THE ACTUAL NAME FOR F NEEDS GR02C910

c TO BE DECLARED E X T E R N A L IN THE GR020920

c CALLING PROGRAM. GR020930

c THE FUNCTION G(X) IS DEFINED TO BE GR020940

c G(X)=F(PHI(X))*PHIP(X) GR020950

c WHERE PHI (X) IS THE CUBIC GIVEN BY GR020960

c THE ARITHMETIC STATEMENT FUNCTION BELOW. GR020970

c PHIP(X) IS ITS DERIVATIVE. THE VARIABLES GR020980

c XL AND XR ARE THE LEFT AND RIGHT ENDPOINTS GR020990

c OF A PARENT INTERVAL OF WHICH (A,B) IS A PART. GR021000

c GR021010

c A - DOUBLE PRECISION GR021020

c LOWER LIMIT OF INTEGRATION GR021030

c GR021040

c B - DOUBLE PRECISION GR021050

c UPINTG LIMIT OF INTEGRATION GR021060

c GR021070

c XL - DOUBLE PRECISION GR021080

c LEFT ENDPOINT OF PARENT INTERVAL GR021090

c GR021100

c XR - DOUBLE PRECISION GR021110

c RIGHT ENDPOINT OF PARENT INTERVAL GR021120

c GR021130

c ON RETURN GR021140

c R - DOUBLE PRECISION GR021150

c APPROXIMATION TO THE INTEGRAL I GR021160

c R IS COMPUTED BY APPLYING THE 15-POINT GR021170

c KRONROD RULE (RESK) OBTAINED 3Y OPTIMAL GR021 180 c ADDITION OF ABSCISSAE TO THE 7-POINT GAUSS GR021190

c RULE (RESG)

.

GR021200

c GR021210

c AE - DOUBLE PRECISION GR021220

c ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR, GR021230

c WHICH SHOULD NOT EXCEED DABS(I-R) GR021240

c GR021250

c RA - DOUBLE PRECISION GR021260

c APPROXIMATION TO THE INTEGRAL J GR021270

c GR021280

c RASC - DOUBLE PRECISION GR021290

c APPROXIMATION TO THE INTEGRAL OF DA3S(G- I/(B-A) )GR021300

c OVER (A, B) GR021310

c GR021320

c FMAX, FMIN - DOUBLE PRECISION GR021330

c MAX AND MIN VALUES OF THE FUNCTION F ON (A,B) GR021340

c SUBROUTINES OR FUNCTIONS NEEDED GR021350

c - F (USER-PROVIDED FUNCTION) GR021360

c - R1MACH GR021370

c - FORTRAN DABS, DMAX1 , DMIN1 GR021380

c GR021390

c . GR021400

C***END PROLOGUE GR021410

c GR021420

DOUBLE PRECISION A,AE,B,DHLGTH,EPMACH, F, FC, FSUM, FVAL1 , FVAL2, GR021430

* FV1 , FV2,HLGTH,RA,RASC,RESG,RESK,RESKH,R,R1MACH,UFLCW, GR021440
* WG,WGK,XGK

DOUBLE PRECISION XL, XR, CENTR, ABSC,U, FMAX, FMIN,PHI , PHIP

INTEGER J , JTW, JTWM1

C

CC

C

C

C

C

C

C

C

C

C

C

C

C

C

DIMENSION FV1(7),FV2(7),WG(4),WGK(8),XGK(8)

THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1) BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR CORRESPONDING WEIGHTS ARE GIVEN.

XGK - ABSCISSAE OF THE 15-POINT KRONROD RULE XGK(2) , XGK(4) , ... ABSCISSAE OF THE 7-POINT

GAUSS RULE

XGK<1), XGK(3) , ... ABSCISSAE WHICH ARE OPTIMALLY

ADDED TO THE 7-POINT GAUSS RULE

WGK - WEIGHTS OF THE 15-POINT KRONROD RULE

WG - WEIGHTS OF THE 7-POINT GAUSS RULE

DATA XGK( 1 ) ,XGK(2) ,XGK(3) ,XGK(4) ,XGK(5) ,XGK(6) ,XGK(7) ,XGK(8)/

0 .991455371 1 2081 26D+00 , 0 . 9491 0791 23427585D+00

,

0 . 8648644233597691 D+00 , 0.7415311 855993944D+00

,

0 . 586087235467691 1 D+00 , 0 . 405845 1 5 1 3773972D+00

,

0.2077849550078985D+00, 0.0D+00 / DATA WGK( 1 ) ,WGK(2) ,WGK(3) ,WGK(4) ,WGK(5) , WGK(6) ,WGK(7) ,WGK(8)/

0 . 2293532201 052922D - 01 , 0 . 6309209262997855D - 01

,

0. 1047900103222502D+00, 0. 14065325971552590+00,

0 . 1 690047266392679D+00 , 0.1 903505780647854D+00

,

* 0.2044329400752989D+00,

DATA WG(1),WG(2),WG(3),WG(4)/

* 0.1 294849661 688697D+00,

* 0.3818300505051 1890+00,

0.20948214108472780+00/

0 . 279705391 4892767D+00

,

0.41 79591 836734694D+00/

DATA EPMACH , UFL0W/Z341 0000000000000 , Z001 0000000000000/

PHI (U)=XR- (XR-XL)*U*U*(2.D0*U+3.D0)

PHIP(U)=-6.D0*U*(U+1 .DO)

C

c

c

LIST OF MAJOR VARIABLES

c

c CENTR - MID POINT OF THE INTERVAL

c HLGTH - HALF-LENGTH OF THE INTERVAL

c ABSC - ABSCISSA

c FVAL* - FUNCTION VALUE

c RESG - R OF THE 7-POINT GAUSS FORMULA

c RESK - R OF THE 15-POINT KRONROD FORMULA

c

c

RESKH - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B)

I.E. TO I/CB-A)

c

c

MACHINE DEPENDENT CONSTANTS

C

C EPMACH IS THE LARGEST RELATIVE SPACING.

C UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.

GR021450

GR021460

GR021470

GR021480

GR021490

GR021500

GR021510

GR021520

GR021530

GR021540

GR021550

GR021560

GR021570

GR021580

GR021590

GR021600

GR021610

GR021620

GR021630

GR021640

GR021650

GR021660

GR021670

GR021680

GR021690

GR021700

GR021710

GR021720

GR021730

GR021740

GR021750

GR021760

GR021770

GR021780

GR021790

GR021800

GR021810

GR021820

GR021830

GR021840

GR021850

GR021860

GR021870

GR021880

GR021890

GR021900

GR021910

GR021920

GR021930

GR021940

GR021950

GR021960

GR021970

GR021980

GR021990
c GR022000

C***F IRST EXECUTABLE STATEMENT GR022010

C EPMACH = DBLE( 16**( - 13)

)

GR022020

C UFLOW = DBLE(16**( -65)) GR022030

c GR022040

HLGTH = 0.5D+00*(B-A) GR022050

CENTR = A+HLGTH GR022060

DHLGTH = OABS(HLGTH) GR022070

c GR022080

c COMPUTE THE 15-POINT KRONROD APPROXIMATION TO GR022090

c THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR. GR022100

c GR0221 10 U= ( CENTR -XR)/(XR- XL) GR022120

FMIN=F(PHI (U)

)

GR022130

FMAX=FMI

N

GR022140

FC=FMIN*PHIP(U) GR022150

RESG = FC*UG(4) GR022160

RESK = FC*WGK(8) GR022170

RA = DABS(RESK) GR022180

DO 10 J=1 ,3 GR022190

JTU = J*2 GR022200

ABSC = HLGTH*XGK( JTU) GR022210

U=( CENTR -ABSC-XR)/(XR- XL) GR022220

FVAL1 = F(PHI(U)) GR022230

FMAX=DMAX 1 ( FMAX , FVAL1

)

GR022240

FMIN=DMIN1 (FMIN, FVAL1

)

GR022250

FVAL1=FVAL1*PHIP(U) GR022260

U=(CENTR+ABSC-XR)/(XR-XL) GR022270

FVAL2=F(PHI (U)

)

GR022280

FMAX=DMAX 1 ( FMAX , FVAL2) GR022290

FMIN=DMIN1 (FMIN, FVAL2) GR022300

FVAL2=FVAL2*PHIP(U) GR022310

FVI(JTW) = FVAL1 GR022320

FV2CJTW) = FVAL2 GR022330

FSUM = FVAL1+FVAL2 GR022340

RESG = RESG+WGC J)*FSUM GR022350

RESK = RESK+WGK( JTU)* FSUM GR022360

RA = RA+WGK( JTU)*(DABS( FVAL1 )+DABS( FVAL2)

)

GR022370

10 CONTINUE GR022380

DO 15 J = 1,4 GR022390

JTWM1 = J*2-

1

GR022400

ABSC = HLGTH*XGK( JTWM1

)

GR022410

U= ( CENTR - ABSC -XR)/ (XR-XL) GR022420

FVAL1=F(PHI <U)

)

GR022430

FMAX=DMAX1 ( FMAX, FVAL1

)

GR022440

FMIN=DMIN1 (FMIN, FVAL1

)

GR022450

FVAL1=FVAL1*PHIP(U) GR022460

U=(CENTR+ABSC-XR)/( XR-XL) GR022470

FVAL2=F(PHI(U)) GR022480

FMAX=DMAX1 ( FMAX, FVAL2) GR022490

FMI N=DMIN1 ( FMIN , FVAL2) GR022500

FVAL2=FVAL2*PHIP(U) GR022510

FVI(JTWMI) = FVAL1 GR022520

FV2( JTWM1 ) = FVAL2 GR022530

FSUM = FVAL1+FVAL2 GR022540
RESK = RESK+WGK( JTWM1 )*FSUM

RA = RA+WGK( JTWM1 )*(DABS( FVAL1 )+OABS( FVAL2)

)

15 CONTINUE

RESKH = RESK*0.5D+00

RASC = WGK(8)*DABS(FC-RESKH)

DO 20 J=1 ,7 RASC = RASC+WGK( J )*(DABS( FV1 ( J

)

- RESKH )+DABS( FV2( J

)

- RESKH )

)

20 CONTINUE

R = RESK*HLGTH

RA = RA*DHLGTH

RASC = RASC*DHLGTH

AE = DABS( (RESK-RESG)*HLGTH)

I F ( RASC . NE . 0 . 0D+00 . AND . AE . NE . 0 . 0D+00

)

* AE = RASC’DMINI (0. 1D+01

,

* (0.2D+03*AE/RASC)**1 .5D+00)

I F(RA.GT.UFLOW/(0.5D+02*EPMACH) ) AE = DMAX1

* ( (EPMACH*0.5D+02)*RA,AE)

RETURN

END

C SUBROUTINE QUADRT ( A , ROOT 1

, ROOT2 , TOL , NOROOT

)

C

C

C TO SOLVE ANY QUADRT I C EQUATION OF THE FORM

C

C A(1 )*X**2+A(2)*X+A(3) = 0 C DOUBLE PRECISION A, ROOT 1 , ROOT2,X,Y,Z,U,ZERO

DIMENSION A(3) ,ROOT1 (2) , ROOT2L2)

C ZERO = TOL/IO.O

C

IF (DABS(A(1 )) -ZERO) 2, 2,

1

2 IF (DABS(A(2) )-ZERO) 3,3,4

3 NOROOT =0 RETURN

C

1 NOROOT = 2 X = A(2)*A(2) - 4.0*A( 1 )*A(3)

Y = A( 1 ) + A( 1

)

Z = DSQRT (DA8S(X) )/Y U = -A(2)/Y

C

IF (X.LT.0.0) GO TO 7

ROOT1 ( 1 )=W+Z

ROOT1 (2)=0.0

ROOT2( 1 )=W-Z

ROOT2(2)=0.0

RETURN

C 7 ROOT1 ( 1 )=W ROOT1 (2)=Z

ROOT2( 1 )=W

GR022550

GR022560

GR022570

GR022580

GR022590

GR022600

GR022610

GR022620

GR022630

GR022640

GR022650

GR022660

GR022670

GR022680

GR022690

GR022700

GR022710

GR022720

GR022730

GR022740

GR022750

GR022760

GR022770

GR022780

GR022790

GR022S00

GR022810

GR022820

GR022830

GR022840

GR022850

GR022860

GR022870

GR022880

GR022890

GR022900

GR022910

GR022920

GR022930

GR022940

GR022950

GR022960

GR022970

GR022980

GR022990

GR023000

GR023010

GR023020

GR023030

GR023040

GR023050

GR023060

GR023070

GR023080

GR023090
ROOT2(2)=-Z

RETURN

C

4 N0R00T=1 ROOTI(I) = -A(3)/A(2)

ROOT1 (2) = 0.0 RETURN

END £•»***************•***************************************************

cc**

c

c

cc

c

c

c

ccc

c

cc

c

c

c

c

c

cc

cc

c

c

c

c

c

c

cc

cc

cc

cc

cc

c

c

c

SUBROUTINE DCRTNI (X, F,DERF,XST,EPS, IEND, IER)

****************************

PURPOSE

TO SOLVE PHI = 0 BY MEANS OF NEUTON-S ITERATION METHOD.

USAGE

CALL DCRTNI (X, F, DERF, XST, EPS, IEND, IER)

DESCRIPTION OF PARAMETERS

X - DOUBLE COMPLEX RESULTANT ROOT OF EQUATION F(X)=0.

F - DOUBLE COMPLEX RESULTANT FUNCTION VALUE AT ROOT X. DERF - DOUBLE COMPLEX RESULTANT VALUE OF DERIVATIVE

AT ROOT X. XST - DOUBLE COMPLEX INPUT VALUE WHICH SPECIFIES THE INITIAL GUESS OF THE ROOT X. EPS - SINGLE PRECISION INPUT VALUE WHICH SPECIFIES THE UPPER BOUND OF THE ERROR OF RESULT X. IEND - MAXIMUM NUMBER OF ITERATION STEPS SPECIFIED.

IER - RESULTANT ERROR PARAMETER CODED AS FOLLOWS

IER=0 - NO ERROR,

IER=1 • NO CONVERGENCE AFTER IEND ITERATION STEPS

IER=2 - AT ANY ITERATION STEP DERIVATIVE DERF WAS EQUAL TO ZERO.

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED PHIO.AND DPHIO

METHOD SOLUTION OF EQUATION F(X)=0 IS DONE BY MEANS OF NEWTON-S

ITERATION METHOO, WHICH STARTS AT THE INITIAL GUESS XST OF

A ROOT X. CONVERGENCE IS QUADRATIC IF THE DERIVATIVE OF

F(X) AT ROOT X IS NOT EQUAL TO ZERO. ONE ITERATION STEP REQUIRES ONE EVALUATION OF F(X) AND ONE EVALUATION OF THE DERIVATIVE OF F<X).

FOR REFERENCE, SEE R. ZURMUEHL, PRAKTISCHE MATHEMATIK FUER INGENIEURE UND PHYSIKER, SPRINGER, BERLIN/GCETTINGEN/

HEIDELBERG, 1963, PP.12-1 7.

DOUBLE COMPLEX X,F, DERF, XST, TOL,DX,PHIO, DPHIO

GR023100

GR0231 10 GR023120

GR023130

GR023140

GR023150

GR023160

GR023170

GR023180

GR023190

GR023200

GR023210

GR023220

GR023230

GR023240

GR023250

GR023260

GR023270

GR023280

GR023290

GR023300

GR023310

GR023320

GR023330

GR023340

GR023350

GR023360

GR023370

GR023380

GR023390

GR023400

GR023410

GR023420

GR023430

, GR023440

GR023450

GR 023460

GR023470

GR023480

GR023490

GR023500

GR023510

GR023520

GR023530

GR023540

GR023550

GR023560

GR023570

GR023580

GR023590

GR023600

GR 0236 10 GR023620

GR023630

GR023640
DOUBLE PRECISION A,DTOL,DTOLF

EXTERNAL PHI0,DPHI0

C PREPARE ITERATION

IER=0 X=XST TOL=X

F = PHIO(TOL)

DERF = DPHIO(TOL)

DTOLF=100.D0*EPS

C

C

C START ITERATION LOOP

DO 6 1=1 , IEND

I F(CDABS(F))1 , 7,

1

C

C EQUATION IS NOT SATISFIED BY X

1 I F(CDABS(DERF))2,8,2

C

C ITERATION IS POSSIBLE

2 DX=F/DERF

X=X-DX

TOL=X

F = PHIO(TOL)

DERF = DPHIO(TOL)

C

C TEST ON SATISFACTORY ACCURACY DTOL=DBLE(EPS)

A=CDABS(X)

I F(A- 1 .D0)4, 4,3

3 DTOL=DTOL*A

4 I F(CDABS(DX)-DTOL)5,5,6

5 IF(CDABS(F)-DT0LF)7,7,6

6 CONTINUE

C END OF ITERATION LOOP

C

C

C NO CONVERGENCE AFTER IEND ITERATION STEPS. ERROR RETURN.

IER=1 7 RETURN

C

C 8

SUBROUTINE RSORT(A,B,R,N,M,MS)

C

DIMENSION A(1),B(1),R(1)

C

C MOVE SORTING KEY VECTOR TO FIRST COLUMN OF OUTPUT MATRIX

C AND 3UILD ORIGINAL SEQUENCE LIST IN SECOND COLUMN

C DO 10 1=1,

N

ERROR RETURN IN CASE OF ZERO DIVISOR

IER=2 RETURN

END

GR023650

GR023660

GR023670

GR023680

GR023690

GR023700

GR023710

GR023720

GR023730

GR023740

GR023750

GR023760

GR023770

GR023780

GR023790

GR023800

GR023810

GR023820

GR023830

GR023840

GR023850

GR023860

GR023870

GR023880

GR023890

GR023900

GR023910

GR023920

GR023930

GR023940

GR023950

GR023960

GR023970

GR023980

GR023990

GR024000

GR024010

GR024020

GR024030

GR024040

GR024050

GR024060

GR024070

GR024080

GR024090

GR024100

GR024110

GR024120

GR024130

GR024140

GR024150

GR024160

GR 024 170 GR024180

GR024190
R( I )=B( I

)

GR024200

I2=I>M GR024210

10 R(I2)=I GR024220

c GR024230

c SORT ELEMENTS IN SORTING KEY VECTOR (ORIGINAL SEQUENCE LIST GR024240

r IS RESEQUENCED ACCORDINGLY) GR024250

C GR024260

L=N+1 GR024270

20 I S0RT=0 GR024280

L=L1

GR024290

I F(L .LT. 2) GO TO 50 GR024300

DO 40 1=2,

L

GR024310

I F(R( I)-R(I-I)) 30,40,40 GR024320

30 ISORT=1 GR024330

RSAVE=R< I

)

GR024340

R( I )=R( 1*1) GR024350

R( I - 1 )=RSAVE GR024360

I2=I+N GR024370

SAVER=R( 12) GR024380

R< 1 2)=R( 12-1

)

GR024390

R ( 1 2 - 1 )=SAVER GR024400

40 CONTINUE GR024410

I F( I SORT) 20,50,20 GR024420

c GR024430

c MOVE ROWS FROM MATRIX A TO MATRIX R (NUMBER IN SECOND COLUMN GR024440

c OF R REPRESENTS ROW NUMBER OF MATRIX A TO BE MOVED) GR024450

c GR024460

50 DO 80 1=1,

N

GR024470

c GR024480

c GET ROW NUMBER IN MATRIX A GR024490

c GR024500

I2=I+N GR024510

I N=R ( I 2

)

GR024520

c GR024530

IR=I -N GR024540

DO 80 J=1,M GR024550

c GR024560

c LOCATE ELEMENT IN OUTPUT MATRIX GR024570

c GR024580

IR=IR+N GR024590

c GR024600

c LOCATE ELEMENT IN INPUT MATRIX GR024610

c GR024620

CALL LOC( IN , J , IA, N,M,MS) GR024630

c GR024640

c TEST FOR ZERO ELEMENT IN DIAGONAL MATRIX GR024650

c GR024660

I FC IA) 60,70,60 GR024670

c GR024680

c MOVE ELEMENT TO OUTPUT MATRIX GR024690

c GR024700

60 R( IR)=A( IA) GR024710

GO TO 80 GR024720

70 R( I R

) =0 GR024730

80 CONTINUE GR024740
RETURN

END

Q-*******************************************************************

C SUBROUTINE LOC( I

, J, IR,N,M,MS)

C

PURPOSE COMPUTE A VECTOR SUBSCRIPT FOR AN ELEMENT IN A MATRIX OF SPECIFIED STORAGE MODE

USAGE

CALL LOC (I,J, IR,N,M,MS)

DESCRIPTION OF PARAMETERS

I

- ROW NUMBER OF ELEMENT

J - COLUMN NUMBER OF ELEMENT

IR - RESULTANT VECTOR SUBSCRIPT

N - NUMBER OF ROWS IN MATRIX

M * NUMBER OF COLUMNS IN MATRIX

MS - ONE DIGIT NUMBER FOR STORAGE MODE OF MATRIX

0 - GENERAL

1

- SYMMETRIC

2 - DIAGONAL

REMARKS

NONE

M IS UNUSED BUT IS REQUIRED TO BE PRESENT

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED

NONE

METHOD

MS=0

MS=1

MS=2

SUBSCRIPT IS COMPUTED FOR A MATRIX WITH N*M ELEMENTS

IN STORAGE (GENERAL MATRIX) SUBSCRIPT IS COMPUTED FOR A MATRIX WITH N*(N+1)/2 IN STORAGE (UPPER TRIANGLE OF SYMMETRIC MATRIX). IF ELEMENT IS IN LOWER TRIANGULAR PORTION, SUBSCRIPT IS CORRESPONDING ELEMENT IN UPPER TRIANGLE. SUBSCRIPT IS COMPUTED FOR A MATRIX WITH N ELEMENTS

IN STORAGE (DIAGONAL ELEMENTS OF DIAGONAL MATRIX).

IF ELEMENT IS NOT ON DIAGONAL (AND THEREFORE NOT IN STORAGE), IR IS SET TO ZERO.

IX=I JX=J IF(MS-I) 10,20,30

10 IRX=N*( JX- 1 )+IX GO TO 36

20 IF(IX-JX) 22,24,24

22 IRX=IX+( JX*JX- JX)/2

GO TO 36

GR024750

GR024760

GR024770

GR024780

GR024790

GR024800

GR024810

GR024820

GR024830

GR024840

GR024850

GR024860

GR024870

GR024880

GR024890

GR024900

GR024910

GR024920

GR024930

GR024940

GR024950

GR024960

GR024970

GR024980

GR024990

GR025000

GR025010

GR025020

GR025030

GR025040

GR025050

GR025060

GR025070

GR025080

GR025090

GR025100

GR0251 10 GR025120

GR025130

GR025140

GR025150

GR025160

GR025170

GR025180

GR025190

, GR025200

GR025210

GR025220

GR025230

GR025240

GR025250

GR025260

GR025270

GR025280

GR025290
24 IRX=JX+( I X* I X - 1 X )/2 GR025300

GO TO 36 GR025310

30 IRX=0 GR025320

IF(IX-JX) 36,32,36 GR025330

32 IRX=IX GR025340

36 IR=IRX GR025350

RETURN GR025360

END GR025370

NBS-114A PEv. 2-«Cl U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO.

NBSIR/ 85—3234

2. Performing Organ, Report No. 3. Publication Date

NOVEMBER 1985

4 . TITLE AND SUBTITLE

Dynamic Green's Functions of An Infinite Plate - A Computer Program

5. AUTHOR(S)

Nelson N. Hsu

6. PERFORMING ORGANIZATION (If joint or other than N BS, see instruction s)

national bureau of standards

DEPARTMENT OF COMMERCE

WASHINGTON. D.C. 20234

7. Contract/Grant No.

3. Type of Report & Period Covered

9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City. State, ZIP)

10. SUPPLEMENTARY NOTES

Document describes a computer program; SF-185, FlPS Software Summary, is attached.

11. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant

bi bl iography or literature survey, mention it here)

This report is a FORTRAN program to compute the Green's functions of an infinite plate. The Green's function, G. . (_£,x,t), is defined as the ith com- ponent of the displacement at x due to anoint force of step-function time

dependency acting at £

in the jth direction initiated at t=0. The Green's

function is the fundamental solution of the transient elastic wave propagation

problem. In general, the displacement field u(^,x,t) at x due to a point force

of arbitrary time dependence acting at can be computed by a convolution

integration; i.e.

,

ivU.X’O = G. . U,x,t)f (T-t)dx.

vJ vJ (1)

Here, G. . is the time derivative of G. . and f.(t) is the point force component

of arbitrary time dependence acting inJ the j t A direction (summation over

repeated indices is used) . Displacement due to point dipoles or couple forces

can be represented by the spatial derivatives of G. .

.

Displacement produced by

a dynamic force distributed over a finite area can

1 llso be computed by nur.er ; :

- integration using the Green's function as the kernel of the integral over the

finite area.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words .-~ Acoustic emission; dynamic displacement due to point impact; Green's function; ray

theory; transient wave propagation; waveform analysis

13. AVAILABILITY

X

,

Unlimited For Official Distribution. Do Not Release to NTIS

Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402.

X' Order From National Technical Information Service (NTIS), Springfield, VA. 22161

14. NO. OF

PRINTED PAGJ

15. Price

C MM*