Program ANRAY - Version 4.50 **************************** Program ANRAY is a modified version of packages ANRAY86 and ANRAY89 by Gajewski and Psencik [1],[2]. It is designed for ray, travel time, ray amplitude, Green's function and ray synthetic seismograms computations. Rays can be computed in two modes. In the first one, rays are specified by the point source location and the initial orientation of the slowness vector at the source. This mode is called initial-value ray tracing. In the second mode, rays are specified by the point source location and a system of regularly or irregularly distributed receivers situated along profiles on surface or an interface, or along vertical line profiles. This mode is called two-point ray tracing. Polarization vectors, geometrical spreading and reflection, transmission and conversion coefficients may be evaluated along the rays. Basic program of the package ANRAY also called ANRAY can produce three files: the first for plotting ray diagrams, travel times and ray amplitudes, the second for computation of synthetic seismograms and the third for plotting plane sections of slowness, phase velocity or group velocity surfaces. Since the program is based on the ray method it does not work properly in the so-called singular regions like vicinity of caustics, of critical points, of transitions from shadow to illuminated region, in anisotropic media in directions in which the two qS waves propagate with nearly the same phase velocity. The existence of the latter singularities does not affect only the amplitude computations but may also cause problems with convergence of iterations to two-point rays. Most important characteristics of the package ANRAY: - the model is 3-D, laterally inhomogeneous, with curved layers of nonzero thickness; - density-normalized elastic parameters inside layers can be determined either by a vertical linear interpolation between interfaces, on which constant values of elastic parameters are specified or by a 3-D interpolation by B-splines with tension from the values of density-normalized elastic parameters specified in grid points of a 3-D grid; - in the B-spline mode, a water layer can be considered; - a point source can be situated anywhere in an isotropic or anisotropic layer; - receivers can be distributed not only along surface or interface profiles, but also along vertical profiles; - receivers can record either displacement vectors or pressure (if pressure is recorded, it is stored in the x-component of the amplitude vector AMPX); - conversion coefficients can be considered at receivers situated on the surface as well as at receivers situated at internal interfaces; - surface and interface profiles do not need to pass through the vertical projection of the point source on the surface or the corresponding interface; - 3-D two-point ray tracing allowing to search for rays arriving from a source to a prescribed vicinity of a receiver specified on a surface, interface or vertical profile, based on paraxial ray approximation; - the output of the program ANRAY can be processed by modified programs from the modified 2-D package SEIS83 [3] for time-domain or frequency-domain computations of ray synthetics. - the package allows to perform calculations of the qS wave synthetics in the quasi-isotropic approximation. - the results of the computations can be displayed in the form of 2-D plots generated in the PostScript format or in the form of 3-D plots using the graphics system of the CRT package. Some useful information - when using ANRAY in the two-point ray-tracing mode, it is desirable that the box containing the model is constructed so that the source is always situated inside the model, away from the vertical boundaries of the model; this makes possible an effective use of the shooting procedure. It is also recommended to avoid receivers situated in the corners of the model box; - the bottom boundary of the model cannot be used as a reflector when ray amplitudes are to be evaluated; calculation of reflection coefficients requires knowledge of parameters below the reflector. Coordinate systems used in the program -------------------------------------- Several coordinate systems are used in the program. The most important ones from the user's point of view are the following: A) "Model" Cartesian coordinate system This coordinate system is used for the description of the model, i.e., for the delineation of interfaces and the density-normalized elastic parameters distribution. It is defined as follows. The x and y coordinate axes are situated in a horizontal plane, with the vertical (z) axis chosen positive downwards. The positive x and y axes are chosen such that the coordinate system is a right handed one. B) "Receiver" cylindrical coordinate system In the two-point ray-tracing mode, i.e. in the mode, in which rays arriving to a priori specified receivers are sought, receivers must be specified. Receivers must be always situated on profiles. These profiles are always situated in a vertical plane. In case of receivers situated along a surface profile, receivers are located along the intersection of the vertical plane with the surface. The same holds for a profile of receivers distributed along an interface. In case of receivers distributed along a vertical profile, the profile is straight vertical line specified by its x- and y-coordinates. The coordinates of receivers are specified in a cylindrical coordinate system with its vertical axis passing through a point specified by the user. In this way, receivers along a surface or interface profile are fully specified by the azimuth of the profile and distances of receivers from the vertical axis. In case of a vertical profile, receivers are fully specified by x- and y-coordinates of the vertical profile and by "model" z-coordinates of receivers. C) "Elastic tensor" Cartesian coordinate system Tensor of density-normalized elastic parameters may be specified in a coordinate system whose axes coincide with the principal axes of the elasticity tensor. In such a case, in anisotropic media with higher symmetry many of the 21 density-normalized elastic parameters are zero. Since the axes of the "elastic tensor" coordinate system may differ from the "model" coordinate system an additional information specifying how the "elasticity tensor" coordinate system is to be rotated within the model must be supplied, see description of input data for isosurface interpolation sub 7a. The x and y coordinate axes are situated in a horizontal plane, with the vertical (z) axis chosen positive downwards. The positive x and y axes are chosen such that the coordinate system is right handed. Description of the model and the source --------------------------------------- The model is situated in a "box" bounded by plane vertical boundaries parallel to the (x,z) and (y,z) planes, and by, possibly curved, top and bottom boundaries formed by interfaces. The model consists of a system of layers separated by curved non-intersecting interfaces. The layers as well as interfaces are numbered sequentially from the top to the bottom. All interfaces must be defined on the whole (x,y) cross-section of the model. The interfaces are specified by points on a rectangular grid completely covering the (x,y) model cross-section. Thus, the first and last lines of this grid must be situated on vertical boundaries of the model. Bicubic spline interpolation is used to approximate the interfaces between the grid points. The interfaces are therefore smooth without corner points. The medium inside layers may be either isotropic or anisotropic. All 21 density-normalized elastic parameters and density may vary spatially within the anisotropic layers and the P and S wave velocities and density may vary spatially in the isotropic layers. The density- normalized elastic parameters have the dimensions of velocity squared. The P and S wave velocities in an isotropic layer are, for simplicity, also referred to as density-normalized elastic parameters. Two methods of the determination of the density-normalized elastic parameters in individual layers may be used. (a) Linear interpolation of density-normalized elastic parameters along vertical lines between interfaces representing isosurfaces of density-normalized elastic parameters. In this approach, the distribution of density-normalized elastic parameters is laterally varying within the layers if their boundaries are curved or inclined. The vertical gradient of the density-normalized elastic parameters within a layer is large when the interfaces bounding the layer are close to one another and vice versa. (b) B-spline interpolation of density-normalized elastic parameters from their values specified in grid points of a 3-D grid. B-splines with tension by A.K.Cline [4]. The length and velocity units km and km/s, m and m/s or m and m/ms must be used consistently. The density should be always specified in g/cm**3 (or 1000 kg/m**3). A point source may be situated at any point of the model (for the effective functioning of the two-point ray-tracing procedure it is recommended that it is situated inside the model box, away from its vertical boundaries). If it is situated at an interface, it is automatically considered to be situated in the underlying layer. A unit isotropic radiation function is incorporated in the program ANRAY. Radiation functions of single force, double couple and explosive (implosive) point sources may be introduced during further processing of the data obtained from the program ANRAY in the programs ANRAYPL, SYNTAN or FRESAN. The polarization vectors at the point source situated in an isotropic layer are specified in the following way. The initial polarization vector of a P wave is specified by a unit vector tangent to the ray at the source. The initial polarization vectors of an S wave are specified by unit vectors situated in the plane perpendicular to the ray. The first of these vectors is situated along the intersection of a vertical plane with the plane perpendicular to the ray. It is orientated so that it has always a negative component into the z axis in the "model" coordinates The other unit vector is chosen such that the triplet formed by the first and second unit vectors in the plane perpendicular to the ray and the vector tangent to the ray at the source form a right-handed system. For this choice, the second unit vector in the plane perpendicular to the ray is horizontal. Description of elementary waves and their codes ----------------------------------------------- Various refracted, multiply reflected and converted elementary waves may be considered. Each of these waves is described by a unique code which can be manually generated. The codes are defined as follows. The whole ray is divided into segments, each of which lies between two successive points at which the ray strikes an interface. If the endpoints of a segment lie on different interfaces, the segmnent is called a simple segment. If the end points lie on the same interface, the segment is called compound segment. Any compound segment is regarded as to be formed of two simple segments. The part of the ray between the source and the point at which the ray first strikes an interface is considered as the first segment of the ray. If the source is situated on an interface and the ray strikes first this interface from below, the first segment is regarded as a compound segment. The code of a wave consists of chain of doublets corresponding to individual simple segments of its rays, successively, from the source towards the termination points. The first number in the doublet gives the number of the layer, in which the corresponding segment is situated. The second number specifies the type of the wave along the segment. The numbers 1 and 2 correspond to quasishear waves (in an isotropic layer it is sufficient to specify the number 1 for the shear wave, the number 2 gives the same wave) the number 3 corresponds to the quasicompressional wave (or compressional wave in an isotropic layer). Description of the ray computations ----------------------------------- Individual rays are specified by the source position and by two angles specifying the orientation of the slowness vector at the source. For simplicity, one of the angles is called here azimuth, and the other declination. Both angles are specified in "model" coordinate system. Let us note that in anisotropic layers, the above angles generally differ from the take-off angles of rays. In isotropic layers, the angles specifying the slowness vector coincide with the take-off angles of corresponding rays. The azimuth is measured in the (x,y) plane, positively clockwise from the positive x axis. The declination is defined as an angle between the slowness vector and its projection into the horizontal plane. It is measured positively clockwise, zero declination corresponding to a ray leaving the source with horizontal slowness vector. The program ANRAY may be run in two different modes: (a) Initial-value ray-tracing mode, (b) Two-point ray-tracing mode. In the initial-value ray-tracing mode, the rays are computed for a system of specified initial angles, without a priori knowledge of where they will terminate. In the two-point ray-tracing mode, the rays which terminate in a prescribed vicinity of receivers regularly or irregularly distributed along a surface, interface or vertical profile are computed. Two-point ray tracing is based on the shooting method. Modified versions of the routine TRAMP from the program package SEIS83 [2] are used for this purpose. In both modes, a system of initial angles must be specified, both azimuths and declinations. In the initial-value ray tracing mode, ray paths and related quantities corresponding directly to these initial angles are computed. In the two-point ray tracing mode, the system of initial angles serves only as a basic system, within which an iterative search for rays terminating at the specified receivers is performed. The procedure for finding rays terminating in a vicinity of prescribed receivers works as follows. First, search is made within an iterative procedure (routine PROFIL) for a ray terminating on the profile containing the receiver. During this search declination is held constant. If a ray terminating on the profile is found, it enters into the second iterative procedure (routine RECEIV), which deals only with rays terminating on the profile. The second iterative procedure searches for a ray terminating at the receiver. The iterative procedure is based on the paraxial ray approximation. Description of input and output files ------------------------------------- Main input data are read from the standard input by list-directed input (free format) and consist of a single line containing following data: 'LIN' 'LOU' 'LU1' 'LU2' 'LU3' 'LUBRD' 'LUGRD' 'LUIND' 'LURAY'/ Here: 'LIN' is the name of the input data file LIN. 'LOU' is the name of the output log file LOU. 'LU1' is the name of the output data file LU1. 'LU2' is the name of the output data file LU2. 'LU3' is the name of the output data file LU3. 'LUBRD' is the name of the output data file used in generating a VRML file. 'LUGRD' is the name of the output data file used in generating a VRML file. 'LUIND' is the name of the output data file used in generating a VRML file. 'LURAY' is the name of the output data file used in generating a VRML file. / is a slash recomended in batch and script files to enable future extensions. Defaults: 'LIN'='anray.dat' 'LOU'='anray.out' 'LU1'='lu1.dat' 'LU2'='lu2.dat' Example of the main input data: 'anray.an' / (B-spline interpolation; modbs.for) or 'anray.qi' / (B-spline interpolation; modbs.for) or 'anray.chs' / (linear interpolation; modis.for) or 'anray.vel' / (B-spline interpolation; modbs.for) or 'anray.rfr' / (linear interpolation; modis.for) The data for the program ANRAY are read in from the file LIN, specified by the user. Output data describing the computations are stored in the file LOU. If the file LU1 is specified, output data for plotting ray diagrams, travel time curves and amplitude-distance curves in the program ANRAYPL are stored in it. If the file LU2 is specified, output data for computing ray synthetic seismograms in the program SYNTAN are stored in it. If the file LU3 is specified, output data for plotting sections of velocity surfaces are stored in it. If the files LUBRD, LUGRD, LUIND and LURAY are specified, output data for plotting model box, interfaces, medium parameter distribution and rays in 3-D graphics mode are stored in them for the use in the program ANRAYWRM. Input data in the file LIN -------------------------- 1) MTEXT MTEXT... description of the problem under consideration, (alphanumeric text). Default 'ANRAY'. 2) INULL,ISURF INULL... in some tests (e.g. on zero elements of matrices), the numbers less than or equal RNULL=10**(-INULL) are taken as zero. Default INULL=4. ISURF... dummy parameter 3) MPRINT,NINT,N1,N2 MPRINT... controls storage of the description of the model in the file LOU. MPRINT=0: Only a copy of input data, no other description of the model is stored. Default value. MPRINT=1: Writes the most important input data plus printer plots of the form of interfaces in the model (see ZMIN, ZMAX in line 5). MPRINT=2: The same as for MPRINT=1 plus description of the unrotated and rotated compressed matrices of the density- normalized elastic parameters. NINT... number of interfaces in the model (including first interface representing the model surface and the last interface representing the model bottom). Default NINT=2. N1,N2... number of points along the first and second horizontal axes for triangulation of interfaces for 3-D images, see routine FACETS. Default N1=10, N2=10. 4) The lines 4a-4d repeat NINT times, each for one interface, from the top to the bottom of the model. 4a)MX,MY MX,MY... number of grid lines x=const and y=const for approximation of the interfaces. 4b)SX(1),...,SX(MX) SX(I)... x-coordinates of grid lines for approximation of the interfaces. 4c)SY(1),...SY(MY) SY(I)... y-coordinates of grid lines for approximation of the interfaces. 4d)Z(1),...,Z(MX*MY) Z(I)... z-coordinates (depths) of the interface at the grid points, successively for x=SX(1) from y=SY(1) to SY(MY), then on a new line for x=SX(2) from y=SY(1) to SY(MY), a.s.o. 5) This line repeats NINT times, each one for one interface, from the top to the bottom of the model. ZMIN,ZMAX ZMIN,ZMAX... the depth range at which the interfaces are displayed in the printer plots. The printer plots consist of systems of digits from 0 to 9 and possibly asterisks covering the whole horizontal extent of the model with 101 points in the x direction and 51 points in the y direction. The digits are determined from the computed z coordinates of interfaces on the 101x51 grid using the relation I=IFIX(10.*(Z-ZMIN)/(ZMAX-ZMIN). Outside the range (ZMIN,ZMAX), the asterisks are printed. These lines should be included (possibly as blank lines) even when MPRINT=0, i.e. when the printer plots are not required. ZMIN of the shallowest interface and ZMAX of the deepest interface are used as the z-coordinates of the upper and the bottom horizontal planes limiting the model box in 3-D images. 6) Cards controlling the approximation of density-normalized elastic parameters and density distribution 6a) ISQRT,IRHO ISQRT... controls the approximation of density-normalized elastic parameters. ISQRT=0: Approximation in density-normalized elastic parameters. Default value. ISQRT=1: Approximation in square roots of density-normalized elastic parameters (velocity approximation). Note: In anisotropic layers (IANI.NE.0), ISQRT=0 automatically. IRHO... controls the density distribution in the model. IRHO=0: The density RHO at any point of the model is determined from the relation RHO=1.7+0.2*VP, where VP is the square root of the density-normalized elastic parameter A11 (which in an isotropic layer corresponds to the P-wave velocity). Default value. IRHO=1: The density is constant throughout each layer with value specified in line 6b. 6b) This line is read only if IRHO.NE.0. RHO(1),RHO(2),...,RHO(NINT-1) RHO(I)... density in the I-th layer. It has meaning only if IRHO=1. RHO(I)=1 by default. 7) Description of distribution density-normalized elastic parameters in individual layers. The input data differ for isosurface interpolation and B-spline approximation. ********************** B-spline approximation ********************** The density-normalized elastic parameters are specified at grid points of a rectangular network. The network may be different in different layers and it must always cover the whole layer. The input data 7a-7d repeat (NINT-1)times, one set for each layer, from the top to the bottom of the model. 7a)IANI,SIGMA IANI... specifies properties of the layer. IANI=0: The layer is isotropic. IANI=1: The layer is anisotropic. Default value. SIGMA... the tension factor, see [4]. This value indicates the curviness desired. If ABS(SIGMA) is nearly zero (e. g. .001) the resulting surface is approximately the tensor product of cubic splines. If ABS(SIGMA) is large (e. g. 50.) the resulting surface is approximately tri- linear. If SIGMA equals zero tensor products of cubic splines result. A standard value for SIGMA is approximately 1. In absolute value. Default SIGMA=0. 7b)NX,NY,NZ NX... number of grid lines in the x-direction (NX.GE.1). For NX=1, the medium is homogeneous in the direction of the x-axis. Default NX=1. NY... number of grid lines in the y-direction (NY.GE.1). For NY=1, the medium is homogeneous in the direction of the y-axis. Default NY=1. NZ... number of grid lines in the z-direction (NZ.GE.1). For NZ=1, the medium is homogeneous in the direction of the z-axis. Default NZ=1. 7c)X1(1),X1(2),...,X1(NX) X1(I)... coordinate of the I-th grid line in the x-direction. Default X1(1)=0. 7c)Y1(1),Y1(2),...,Y1(NY) Y1(I)... coordinate of the I-th grid line in the y-direction. Default Y1(1)=0. 7c)Z1(1),Z1(2),...,Z1(NZ) Z1(I)... coordinate of the I-th grid line in the z-direction. Default Z1(1)=0. 7d) Values of density-normalized elastic parameters at the grid points of the 3-D grid. In an isotropic layer, when IANI=0, values of the squares of the P and S wave velocities are read in successively. In case of an anisotropic layer, individual density-normalized elastic parameters are read in successively. The density-normalized elastic parameters are ordered in the following way: A11,A12,A13,A14,A15,A16, A22,A23,A24,A25,A26,A33,A34,A35,A36,A44,A45,A46,A55,A56,A66. Each density-normalized elastic parameter is read in in the following way DO I=1,NX DO J=1,NY W(I,J,K),K=1,NZ W(I,J,K) contains the value of a density-normalized elastic parameter at the grid point (X(I),Y(J),Z(K)). Note1: It is possible to specify a water layer by specifying S-wave velocity equal zero. This option is applicable only in the B-spline approximation (modbs.for). Note2: The maximum number of grid lines to specify the distribution of density-normalized elastic parameters in a single layer cannot exceed, in the present version, 10. In order to change this number, the dimensions of arrays X1,Y1,Z1,W,WG1,WG2,C,VX,VY,VZ and TEMP must be changed correspondingly. The dimension of the array C must be of at least NX*NY*NZ and of the array TEMP of at least 3*MAX(NX, NY, NZ). The values of variables NW1 and NW2 in the beginning of the routine MODEL (program block MODBS.FOR) must be also changed to the new value (NW1 and NW2 are the first two dimensions of the array W in MODBS.FOR; they should be chosen so that NW1.ge.NX and NW2.ge.NY. The total number of grid points in the whole model cannot exceed, in the present version, 100. In order to change this number, the dimensions 100 of the arrays in the common block /EPAR/ must be changed correspondingly. ************************ Isosurface interpolation ************************ The density-normalized elastic parameters are specified at interfaces representing their isosurfaces. The input data 7a-7c repeat (NINT-1)times, one set for each layer, from the top to the bottom of the model. 7a)IANI,ANGU(1),ANGU(2),ANGU(3),ANGL(1),ANGL(2),ANGL(3) IANI... specifies properties of the layer. IANI=0: The layer is isotropic. IANI=1: The layer is anisotropic. Default value. ANGU(1),...,ANGU(3)... the rotation angles, by which the "elasticity tensor" coordinate system is to be rotated at the upper (U) interface bounding the layer. The angles are specified in degrees. ANGU(1) is a rotation angle around the z-axis of the "elasticity tensor" coordinate system. ANGU(1) rotates the x- and y-axes of the "elasticity tensor" system into new positions x' and y'. ANGU(1) is positive if the rotation is anticlockwise. ANGU(2) is a rotation angle around the y'-axis, which rotates the "elasticity tensor" coordinates x' and z into new positions x'' and z'. ANGU(2) is positive if the rotation is anticlockwise. ANGU(3) is a rotation angle around the z'-axis. Default ANGU(1)=0., ANGU(2)=0., ANGU(3)=0. ANGL(1),...,ANGL(3)... the same as for ANGU, but at the lower (L) interface bounding the layer. Default ANGL(1)=0., ANGL(2)=0., ANGL(3)=0. Note 1: The values of rotation angles inside a layer are determined in the same manner as the density-normalized elastic parameters. They are obtained by a linear interpolation along vertical lines between the corresponding values of ANGU and ANGL. Note2: In the case of an isotropic layer the parameters ANGU and ANGL have no meaning. 7b)This line differs for isotropic and anisotropic layer. Anisotropic layer - IANI.NE.0. (A66U(J,K),J=1,6),K=1,6) A66U... density-normalized elastic parameters in a compressed 6x6 form at the upper (U) interface bounding the layer. It is sufficient to specify only right upper corner of the matrix A66U. Isotropic layer - IANI.EQ.0. A66U(1,1),A66U(4,4) A66U(1,1),A66U(4,4)... squares of P- and S-wave velocities at the upper (U) interface bounding the layer. 7c)This line differs for isotropic and anisotropic layer. Anisotropic layer - IANI.NE.0. (A66L(J,K),J=1,6),K=1,6) A66L... density-normalized elastic parameters in a compressed 6x6 form at the lower (L) interface bounding the layer. It is sufficient to specify only right upper corner of the matrix A66U. Isotropic layer - IANI.EQ.0. A66L(1,1),A66L(4,4) A66L(1,1),A66L(4,4)... squares of P- and S-wave velocities at the lower (L) interface bounding the layer. 8) This line is read only if the file LU3 is specified (generation of data for plots of sections of velocity surfaces). NPAR,LAY NPAR... switch indicating which wave surface is to be plotted. NPAR=1: slowness surface. Default value. NPAR=2: phase velocity surface. NPAR=3: group velocity surface. LAY... number of the layer, in which the point specified on line 9 is situated. Default LAY=1. 9) This line is read only if the file LU3 is specified (generation of data for plots of velocity surfaces). X0,Y0,Z0,DDELTA,AZIM X0,Y0,Z0... coordinates of the point for which a wave surface is to be plotted. The point should not be situated closer to the border of the model or interface than two time steps along the ray, DT, see input line 13. Deafault X0=1., Y0=1., Z0=1. DDELTA... step in the angle of the phase normal (in radians) for plotting the section of a selected surface. In case of the slowness surface and phase velocity surface, it is the step with which points of the corresponding section are plotted. In case of group velocity, the density of plotted points changes according to the intensity of the energy flux. Where energy flux is denser, the points will be denser and vice versa. Default DDELTA=0.05. AZIM... azimuth (in radians) of the vertical plane intersecting the studied surface. For AZIM=0., a section by the (x,z) plane is generated. For AZIM=1.570796, a section by the (y,z) plane is generated. Default AZIM=0. 10) ICONT,MEP,MOUT,MDIM,METHOD,MREG,ITMAX,IPOL,IPREC,IRAYPL, IPRINT,IAMP,MTRNS,ICOEF,IRT,ILOC,MCOD,MORI ICONT... controls the continuation of the computations. ICONT=0: Indicates the termination of the computations. The line 10 is the last line in the input data set. ICONT=1: The computations continue, the next line to be read in is the line 10. Default value. Note: For ICONT=0, the rest of the line 10 can be blank. MEP... specifies whether initial-value ray tracing or two-point ray tracing is to be performed. MEP.EQ.0: Initial-value ray tracing is performed. Default value. MEP.NE.0: Two-point ray tracing is performed. ABS(MEP): The number of receivers, at which the rays should terminate. MEP=1: Computation of a ray terminating at a single receiver on a surface of the model or an internal interrface. MEP.LT.0 or MEP.GT.1: The receivers can be distributed along a surface profile, a vertical profile or a profile along an interface, see parameter ILOC in this line; the orientation of the profile is specified by parameter PROF in line 11. MEP.GT.1: The receivers are distributed regularly along the profile (see line 11). MEP.LT.0: The receivers are distributed irregularly along the profile (see line 11). MOUT... controls printing of the results to the file LOU. MOUT=0: Only the input data are reproduced and a list of codes is printed. Default value. MOUT=1: Only the results for the rays terminating at specified receivers, and boundary rays are printed. MOUT=2: The results for all rays terminating on a specified profile are printed. MOUT=3: The results for all rays computed. Note: In case of initial-value ray tracing, MOUT=1 already gives sufficient information about the results of computation. MDIM... controls the level of the computations. MDIM=0: Only rays and travel times along them are computed. Default value. MDIM=1: Spreading-free amplitudes, i.e. products of plane wave reflection/transmission/conversion coefficients are evaluated along the rays in addition to the quantities obtained for MDIM=0. MDIM=2: Ray amplitudes including geometrical spreading METHOD... specifies the method used for the iterative determination of rays in the two-point ray tracing procedure. METHOD=0: Paraxial ray approximation from the rays specified by the basic system of initial angles, see lines 14 and 15 or 16. Default value. METHOD=1: Paraxial ray approximation from the previously found two-point rays (fast but not always reliable). Note 1: In case of initial-value ray tracing (MEP=0), the parameter METHOD has no meaning. MREG... controls the evaluation of amplitudes at the termination points of the rays. MREG=0: Elastic case: displacement vector calculated. The conversion coefficients are applied. Default value. MREG=1: Elastic case: displacement vector calculated. The conversion coefficients are not applied. MREG=2: Acoustic case: pressure calculated. The conversion coefficients are applied. MREG=3: Acoustic case: pressure calculated. The conversion coefficients are not applied. Note: MREG is set to be equal to 1 automatically for VSP computations (when ILOC.EQ.1). In this case elastic case is considered and conversion coefficients are not used. ITMAX... number of iterations permitted in the two-point ray tracing procedure. The default value is ITMAX=10. Note: In case of initial-value ray tracing (MEP=0), the parameter ITMAX has no meaning. IPOL... controls printing of polarization vectors in the routine POLAR. IPOL=0: Polarization vectors are not printed at all. Default value. IPOL=1: Polarization vectors are printed only at the points of incidence of a ray at an interface. IPOL=2: Polarization vectors are printed at all computed points along the ray (this option makes the computations more time consuming). IPREC... switch controlling precission of the ray tracing and dynamic ray tracing calculations. IPREC=0: Precission is not tested. Default value. IPREC=1: Precission is tested and where inaccuracy exceeds the value RNULL (see description of the parameter INULL in input data 2), ray tracing and dynamic ray tracing are modified so that they satisfy the precission tests. IPREC=2: Precission is tested but no modification of results is performed. When inaccuracy exceeds the value RNULL, the values of unsatisfactory tests are printed at any point of any ray where the inaccuracy occurs. Note: As a test of precission the following relations are tested. (a) Scalar product of the slowness vector and the group velocity vector; its value should be p_i*v_i=1. (b) Derivative of the eikonal equation with respect to the two ray parameters r_1 and r_2 (take-off angles of the considered ray}; both derivatives should be zero. (c) Scalar products of the slowness vector and the vectors Q_{i1}=dx_i/dr_1 and Q_{i2}=dx_i/dr_2 at constant travel time; their values should be p_i*Q_{iJ}=0. IRAYPL... dummy parameter. IPRINT,IAMP,MTRNS,ICOEF,IRT... switches controlling storage of auxiliary quantities in the output file LOU. IPRINT controls output of quantities computed along the rays in routine and routine INV; IAMP controls the printing of the auxiliary quantities concerning amplitudes in the routine AMPL; MTRNS - transformation of the slowness vector at the interfaces in the routine TRANSL; ICOEF - reflection/ transmission/conversion coefficient in the routine COEF; IRT - iterative determination of the point of intersection of a ray with an interface in the routine OUT. It is recommended to specify all these parameters zero, which means that no auxiliary quantities are printed. Default values IPRINT=0, IAMP=0, MTRNS=0, ICOEF=0, IRT=0. ILOC... controls the orientation and position of the profile, on which receivers are situated. ILOC=0: Receivers along a surface profile. Default value. ILOC=1: Receivers along a vertical profile. ILOC.GT.1: Receivers along the ILOC-th interface. MCOD... controls specification of initial angles of rays of individual elementary waves. MCOD=0: The system of initial angles of rays is determined once in the beginning of computations. It is the same for all elementary waves. Default value. MCOD=1: The system of initial angles of rays is determined for each elementary wave separately. MORI... controls specification of ray parameters: angles specifying a ray at the source. The two angles are chosen as polar angles related to the "model" axis z and y. The reason for the two options is to avoid singular situations as that which occurs for vertical rays if polar angles are related to the z-axis. See specification of the initial slowness vector by ray angles specified in lines 14 and 15 or 16. MORI=0: The ray angles are related to the "model" z-axis. Default value. MORI=1: The ray angles are related to the "model" y-axis. 11) This line is read only if MEP.NE.0, i.e., in the two-point ray-tracing mode. Specification of a profile with receivers. Specification of receiver positions along this profile. Receivers are specified in the "receiver" cylindrical coordinate system with axis through the point [XREF,YREF]. For MEP.LT.0, NDST=-MEP. PROF,DST(1),DST(2),...,DST(NDST),XPRF,YPRF PROF... azimuthal angle specifying a surface or interface profile, along which the receivers are situated - an angle between the vertical plane containing the profile and the positive "model" x axis. The profile passes through the vertical axis specified by the coordinates XPRF and YPRF. PROF has no meaning in the VSP mode, when receivers are distributed along a vertical profile (ILOC=1). PROF is measured in radians from the positive x axis of the "model" coordinates towards the positive y axis. For example, PROF=0.0 corresponds to the profile situated in the vertical plane containing positive "model" x axis, PROF=1.5707 to the profile situated in the vertical plane containing positive "model" y axis. PROF has no meaning when receivers are distributed along a vertical profile. DST(I)... distances of receivers from the origin of the profile (XPRF,YPRF), see below, measured along the profile. In case of a surface or interface profile, DST(I) is the distance of the I-th receiver from the vertical axis specified by coordinates XPRF and YPRF. In case of a vertical profile, DST(I) is the "model" z-coordinate of the receiver. XPRF,YPRF... x- and y- "model" coordinates of the vertical axis, which contains an origin with respect to which receivers situated along surface or interface profiles are specified. The origin should be chosen so that the receivers are always to one side from it. Default values XPRF=XSOUR, YPRF=YSOUR. For XSOUR and YSOUR, see input data line 13. For MEP=1, NDST=1. XREC,YREC XREC,YREC... x- and y-coordinates of the receiver situated on the surface of the model or an internal interface (the interface is specified by the wave code). The azimuth and the epicentral distance of the receiver from the source are determined automatically and appear in the file LOU. If BMIN, BSTEP and BMAX in line 15 or 17 are not specified, they are specified automatically from the values of the azimuth and the epicentral distance. For MEP.GT.1, NDST=MEP. PROF,RMIN,RSTEP,XPRF,YPRF PROF... azimuthal angle (in radians) specifying a surface or interface profile, along which the receivers are situated - an angle between the vertical plane containing the profile and the positive "model" x axis. The profile passes through the vertical axis specified by the coordinates XPRF and YPRF. PROF has no meaning in the VSP mode, when receivers are distributed along a vertical profile (ILOC=1). PROF is measured in radians from the positive x axis of the "model" coordinates towards the positive y axis. For example, PROF=0.0 corresponds to the profile situated in the vertical plane containing positive "model" x axis, PROF=1.5707 to the profile situated in the vertical plane containing positive "model" y axis. PROF has no meaning when receivers are distributed along a vertical profile. RMIN... distance of the first receiver from the origin of the profile (XPRF,YPRF), see below. In case of a surface or interface profile, RMIN is the distance of the first receiver from the vertical axis specified by coordinates XPRF and YPRF. In case of a vertical profile, RMIN is "model" z-coordinate of the first receiver. RSTEP... the step in distances of neighbouring receivers from the origin (XPRF,YPRF). The distance of the I-th receiver is determined from the relation DST(I)=RMIN+RSTEP*(I-1). XPRF,YPRF... x- and y- "model" coordinates of the vertical axis, which contains an origin with respect to which receivers situated along surface or interface profiles are specified. The origin should be chosen so that the receivers are always to one side from it. Default values XPRF=XSOUR, YPRF=YSOUR. For XSOUR and YSOUR, see input data line 13. 12) This line is read only when ILOC.eq.1 and MEP.ne.0. XVSP,YVSP XVSP,YVSP... x and y coordinates of the vertical profile in case of VSP computations. 13) XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,REPS,PREPS XSOUR,YSOUR,ZSOUR... (x,y,z) coordinates of the source. TSOUR... initial time. Default TSOUR=0. DT... time step for the P wave in the integration of the ray tracing system. For S waves the time step is automatically multiplied by 1.7 to make the speed of computations comparable with P waves. DT should be positive. The default value is DT=1. AC... the required accuracy of the integration of the ray tracing system. Recommended values of this parameter are in the range (0.0001, 0.001). The default value is AC=0.0001. REPS... the vicinity of a receiver on the profile in the length units used for specification of coordinates. Iterative procedure terminates if a ray is found, which has its termination point within this vicinity. Note: In case of initial-value ray tracing (MEP=0), the parameter REPS has no meaning. The default value is REPS=0.05. PREPS... specifies maximum permitted deviation of the termination point of a ray from the profile with receivers in the length units used for specification of coordinates. The deviation is measured as the length of the arc of the circle passing through the termination point, with the centre at epicentre, between the profile and the termination point. Note: In case of initial-value ray tracing (MEP=0), the parameter PREPS has no meaning. The default value is PREPS=0.05. 14) This line is read only when MCOD.EQ.0 AMIN,ASTEP,AMAX 15) This line is read only when MCOD.EQ.0 BMIN,BSTEP,BMAX AMIN,ASTEP,AMAX,BMIN,BSTEP,BMAX... define the basic system of initial angles (in radians) for ray tracing. Default values BMIN=PROF(1)-.3, BSTEP=.6, BMAX=PROF(1)+.4 The initial angles have different meaning in the modes MORI.EQ.0 and MORI.NE.0. For MORI.EQ.0, the angles A and B specify a slowness vector p at the source in the "model" coordinates in the following way: p=c**(-1)*(cosBcosA, sinBcosA, sinA). For MORI.NE.0, the angles A and B specify a slowness vector p at the source in the "model" coordinates in the following way: p=c**(-1)*(cosAcosB, sinB, sinAcosB). The same values of AMIN,ASTEP,AMAX and BMIN,BSTEP,BMAX are used for all considered waves. 16) KC,KREF,CODE(1,1),CODE(1,2),CODE(2,1),CODE(2,2),... ...CODE(KREF,1),CODE(KREF,2) KC... KC.EQ.0: Ray of direct unconverted wave with no reflections and no conversions at interfaces. The ray terminates on boundaries of the model. KC.NE.0: Ray of multiply reflected, possibly converted wave. KC=1: After leaving the source, the ray is to strike first the interface below the source. KC=-1: after leaving the source, the ray is to strike first the interface above the source. KREF... the number of simple segments of the ray. In case of the direct ray (KC=0), specify KREF=1. KREF=0 indicates that no elementary wave is to be computed. CODE(I,1)... the number of the layer in which the I-th simple segment of the ray is situated. CODE(I,2)... the wave type along the I-th simple segment of the ray. CODE(I,2)=1 or 2: quasishear waves (qS1,qS2). CODE(I,2)=3: quasicompressional wave (qP). Note 1: In an isotropic layer, CODE(I,2)=1 must be used for the shear wave. Note 2: CODE(I,J) is an integer. 17) The following two lines are read only when MCOD.NE.0. They are read in separately for each considered wave. AMIN,ASTEP,AMAX BMIN,BSTEP,BMAX The meaning of the above parameters is the same as in lines 14 and 15. Example of data LIN for anisotropic model 'an' Example of data LIN for background isotropic model 'qi' for QI computations Example of data LIN for anisotropic model 'chs' (isosurface approximation) Example of data LIN for generating file for plotting phase velocity sections 'vel' Example of data LIN for generating VRML file 'rfr' Termination of computations --------------------------- The sequence of lines 16-17 can be repeated any number times, separately for each considered elementary wave. To terminate the computations of individual elementary waves put KREF=0 on the line 16 (or simply insert a blank line in the position of the line 16). The sequence of lines 10-17 can be repeated any number times to perform different computations for the same model. Each system of lines 10-17 produces file LU1 for plotting in program ANRAYPL (if ISPACE=0) or for plotting rays in a graphical mode (ISPACE=1), and file LU2 (if LU2.NE.0) for program SYNTAN, in which ray synthetic seismograms are computed. To terminate completely the computations, specify ICONT=0 (or simply insert a blank line in the position of the line 10). Output to the file LU1 ---------------------- For LU1.NE.0 and MEP.NE.0, i.e. in the two-point ray-tracing mode, certain important data concerning rays, travel times and amplitudes are stored in the file LU1 for possible further processing. The file LU1 can be used in the program ANRAYPL, included in this package, for plotting ray diagrams, time-distance and amplitude-distance curves. The data are stored in LU1 in the following formatted form: (A)Information about the source and receivers. 1) ICONT,NDST,ILOC FORMAT(26I3) ICONT... indicative index. ICONT=0: Indication of the end of the data stored in the file LU1. The line 1) is the last data line in the file LU1. ICONT=1: Indication that a data set corresponding to one elementary wave follows. NDST... number of receivers along the considered profile. ILOC... controls the orientation and position of the profile, on which receivers are situated. ILOC=0: Receivers along a surface profile. ILOC=1: Receivers along a vertical profile. ILOC.GT.1: Receivers along the ILOC-th interface. 2) RO FORMAT(8F10.5) RO... density at the source. 3) NPN,NPN,NPN FORMAT(26I3) Dummy data, NPN=2 by default. 4) APN,APN,APN,APN,APN FORMAT(5E15.5) Dummy data, APN=0.0 by default. There are two lines of these dummy data. 5) XPRF,YPRF,0.,PROF FORMAT(8F10.5) XPRF,YPRF... x- and y- "model" coordinates of the vertical axis, with respect to which receivers situated along surface or interface profiles are specified, see input data line 11. PROF... azimuthal angle specifying a surface or interface profile, along which the receivers are situated - an angle between the vertical plane containing the profile and the positive "model" x axis. PROF has no meaning in the VSP mode, when receivers are distributed along a vertical profile (ILOC=1). PROF is measured in radians from the positive x axis of the "model" coordinates towards the positive y axis. For example, PROF=0.0 corresponds to the profile situated in the vertical plane containing positive x axis, PROF=1.5707 to the profile situated in the vertical plane containing positive y axis. 6) (DST(I),I=1,NDST) FORMAT(8F10.5) DST(I)... coordinates of receivers along the profile. In case of a surface or interface profile, DST(I) is the radius of the receiver measured from the vertical axis passing through the source, see the description of the "receiver" coordinate system. In case of a vertical profile DST(I) is z-coordinate of the receiver in the "model" coordinate system. (B)Information about rays. For each ray with a termination point at the specified receiver, the quantities given in the following lines 7 and 8 are stored. The data given in 7 and 8 are stored successively for all rays terminating at specified receivers. After the last ray, the parameter N (see 7) is set to zero. Then, the quantities in section C are stored. 7) N,II FORMAT(26I3) N.NE.0... the number of points along the ray. N.EQ.0... indication of the end of the information about rays. II... specifies the number of the receiver at which the ray terminates. Note: Maximum value of N is 2000. 8) (T(J),X(J),Y(J),Z(J),P1(J),P2(J),P3(J),A11(J),A44(J),RHO(J), (E(1,K,J),K=1,3),(E(2,L,J),L=1,3),J = 1,N) FORMAT(16E15.5) Values of important quantities at the points along rays, starting from the source (J=1) and ending with the termination point (J=N). T - travel time, X,Y,Z coordinates of a point on a ray, P1,P2,P3 - components of the slowness vector, A11,A44 - density-normalized elastic parameters (squares of P and S wave velocities in isotropic media). RHO - density, E(1),E(2) - basis vectors of the ray-centred coordinate system. (C)Information about travel times and amplitudes. 9) INDEX FORMAT(26I3) ABS(INDEX)... the total number of rays of the considered elementary wave, for which travel times and ray amplitudes are stored. INDEX.LT.0 indicates a shear wave in an isotropic medium. Note: Maximum value of ABS(INDEX) is 500. 10)(INDI(I),XCOOR(I),ZCOOR(I),TIME(I),AMPX(1,I),AMPY(1,I),AMPZ(1,I) FORMAT(I5,9E15.5) INDI(I)... specifies the successive number of the receiver at which the I-th ray terminates. XCOOR(I)... radius of the termination point of the ray in the "receiver" cylindrical coordinate system, measured from the vertical axis through the source. ZCOOR(I)... depth of the termination point in the "model" coordinates. TIME(I)... travel time at the termination point of the I-th ray. AMPX(1,I),AMPY(1,I),AMPZ(1,I)... complex amplitudes of the x, y and z components of the polarization vector in the "model" coordinates; for an S wave in an isotropic medium, complex amplitudes corresponding to the first polarization vector in the plane perpendicular to the ray, see the description of the model and the source in the introduction. If pressure is recorded, it is stored in AMPX. 11) This line is stored only if INDEX.LT.0. AMPX(2,I),AMPY(2,I),AMPZ(2,I) FORMAT(6E15.5) AMPX(2,I),AMPY(2,I),AMPZ(2,I)... complex amplitudes of the components of the second polarization vector of an S wave, see the description of the model and the source in the introduction. 12)(P(I),I=1,3) FORMAT(6E15.5) P(I)... components of the slowness vector of the considered wave at the source. 13)(POL(I),I=1,3) FORMAT(6E15.5) POL(I)... components of the polarization vector of the considered wave at the source. In case of an S wave in an isotropic medium, components of its first polarization vector in the plane perpendicular to the ray, see the description of the model and the source in the introduction. 14)This line is stored only if INDEX.LT.0. (POL1(I),I=1,3) FORMAT(6E15.5) POL1(I)... components of the second polarization vector of an S wave in an isotropic medium, see the description of the model and the source in the introduction. The lines 10-14 are repeated ABS(INDEX)-times. The sequence of data 1-14 repeats for each considered elementary wave. Indication of the end of information in LU1 is ICONT=0 in 1). Output to the file LU2 ---------------------- For LU2.NE.0 and MEP.NE.0, i.e., in the two-point ray-tracing mode, data necessary for construction of synthetic seismograms are stored in the file LU2. Program SYNTAN, included in this package, can be used for this purpose. The data are stored in LU2 in the following formatted form: 1) MTEXT FORMAT(17A4) MTEXT... alphanumeric test describing the computations. 2) NDST,KSH,ILOC... FORMAT(26I3) NDST... number of receivers along the profile. KSH... dummy parameter, automatically KSH=2. ILOC... controls the orientation and position of the profile, on which receivers are situated. ILOC=0: Receivers along a surface profile. ILOC=1: Receivers along a vertical profile. ILOC.GT.1: Receivers along the ILOC-th interface. 3) XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,ROS FORMAT(8F10.5) XSOUR,YSOUR,ZSOUR... coordinates of the source in the "model" coordinate system. TSOUR... initial time. RSTEP... an average distance between neighbouring receivers. ROS... density at the source. 4) DST(I),I=1,NDST FORMAT(8F10.5) DST(I)... coordinates of receivers along the profile. In case of a surface or interface profile, DST(I) is the radius of the receiver measured from the vertical axis passing through the source, see the description of the "receiver" coordinate system. In case of a vertical profile, DST(I) is z-coordinate of the receiver in the "model" coordinate system. 7) NC,II,T,AMPX1,AMPY1,AMPZ1,TAST FORMAT(2I3,F10.5,12E12.6,F10.6) IABS(NC)... number of the considered elementary wave. NC.LT.0 indicates that the wave starts as an S wave from the source situated in an isotropic layer. II... successive number of the receiver position. T... arrival time. AMPX1,AMPY1,AMPZ1... complex amplitudes of the x, y and z components of the polarization vector in the "model" coordinates; for an S wave in an isotropic medium, complex amplitudes corresponding to the first polarization vector in the plane perpendicular to the ray, see the description of the model and the source in the introduction. If pressure is recorded, it is stored in AMPX. TAST... the global absorption factor t*; in this version, automatically t*=0. 8) This line is stored only if NC.LT.0. AMPX2,AMPY2,AMPZ2 FORMAT(6E12.6) AMPX2,AMPY2,AMPZ2... complex amplitudes of the components of the second polarization vector of an S wave, see the description of the model and the source in the introduction. 9) (P(I),I=1,3) FORMAT(3F10.5) P(I)... components of the slowness vector of the considered wave at the source. 10)(POL(I),I=1,3) FORMAT(3F10.5) POL(I)... components of the polarization vector of the considered wave at the source. In case of an S wave in an isotropic medium, components of its first polarization vector in the plane perpendicular to the ray, see the description of the model and the source in the introduction. 11)This line is stored only if NC.LT.0. (POL1(I),I=1,3) FORMAT(3F10.5) POL1(I)... components of the second polarization vector of an S wave in an isotropic medium, see the description of the model and the source in the introduction. The lines 7-11 are repeated ABS(NC)-times. Output to the file LU3 ---------------------- For LU3.NE.0, data necessary for construction of sections of slowness surfaces, phase velocity and group velocity surfaces are stored in the file LU3. Program VELPL, included in this package, can be used for plotting the sections. The data are stored in LU3 in the following formatted form: 1) ITYPE,NPAR,INUM FORMAT(16I5) ITYPE... ITYPE=1: qS1 wave. ITYPE=2: qS2 wave. ITYPE=3: qP wave. NPAR... NPAR=1: plot of section of slowness surface. NPAR=2: plot of section of phase velocity surface. NPAR=1: plot of section of group velocity surface. INUM... number of samples in the plot. 2) DDELTA,XX(I),YY(I),ZZ(I) FORMAT(4E15.5) DDELTA... angle (in radians), for which the vector (XX,YY,ZZ) is calculated. XX(I),YY(I),ZZ(I)... x,y and z components of the slowness vector (if NPAR=1), phase velocity vector (if NPAR=2) or group velocity vector (if NPAR=3) for the angle DDELTA. Maximum points: 300. The line 2 is repeated for all the angles of all the three waves. Output to the file LOU ---------------------- Storage of results in the file LOU is controlled by the switches MPRINT and MOUT. The switch MPRINT controls storage of parameters describing the model. The switch MOUT controls the storage of quantities during the initial-value or two-point ray tracing. MPRINT: MPRINT=0: Only a copy of input data, no other description of the model is presented. MPRINT=1: Writes the most important input data plus printer plots of the form of interfaces in the model (see ZMIN,ZMAX in line 5). MPRINT=2: The same as for MPRINT=1 plus description of the unrotated and rotated compressed matrices of the density- normalized elastic parameters. MOUT: MOUT=0: Only input data are reproduced and list of the codes of the computed waves. MOUT=1: For each code, results describing all rays in the case of initial-value ray tracing and the rays arriving at specified receivers in the case of two-point ray tracing are stored. In the latter case, in addition, all boundary and critical rays are stored. All the results are stored in the routine RECEIV. The results for a ray terminating at a prescribed receiver are stored in two lines as follows: DD,R,XX,YY,ZZ,T,DEC,AZIM,IND,IND1,ITER,II,INDEX UU,(APX(I),APY(I),APZ(I),I=1,2),SPR,KMAH DD: Coordinate of the receiver (in the "receiver" coordinate system), at which the ray should arrive. For initial- value ray tracing, DD has no meaning, therefore, DD=0.0 by default. R: In the two-point ray tracing to receivers situated along the surface or an interface, R is the radius of the termination point in the "receiver" coordinates. In the two-point ray tracing to receivers distributed along a vertical profile, R is z-coordinate in "model" coordinates. R should differ from DD by less then REPS, see line 13 in input data file LIN. XX,YY,ZZ: "Model" coordinates of the termination point of the ray. T: Corresponding arrival time. AZIM,DEC: Azimuth and declination specifying the slowness vector of the considered ray at the source (in radians). IND: Parameter describing the history of the ray. The most important values of IND are: (1 or 2) - the ray intersects the vertical (y,z) boundary of the model with lower (1) or larger (2) value of the coordinate x. (3) - termination of the ray on the model surface. (4) - termination of the ray on the lower boundary of the model. (5) - in the Runge-Kutta method, the basic time step DT is halved more than 10 times, without reaching the required accuracy AC, see line 13. There is probably something wrong with your velocity distribution in the model. (6) - in line 13, DT=0 (not allowed). The computation of the ray diagram does not start. (7) - in line 13, DT.lt.0 (not allowed). The computation of the ray diagram does not start. (8) - termination of the ray at an interface because its ray path does not correspond to the prescribed code. (9) - overcritical incidence of the ray. Termination of calculation of the ray. (10) - coinciding values of the phase velocities of quasishear waves for the given direction. Termination of calculation of ray. (11) - matrix for the evaluation of R/T coefficients singular. Termination of calculation of ray. (12) - travel time along the ray exceeded the value TMAX. This value is set TMAX=10000. Termination of calculation of the ray. (13) - the ray is specified by more than 2000 points. Its computation is terminated. (14) - specified code does not correspond to the position of the source in the model. (15) - prescribed reflection from the bottom boundary of the model while MDIM.NE.0. Termination of calculation of the ray. (20) - two neighbourhood interfaces in the model intersect each other. The computation terminates. (30 or 31) - the ray intersects the vertical (x,z) boundary of the model with lower (30) or larger (31) value of the coordinate y. (39) - number of iterations necessary for the determination of the point of intersection of a ray with an interface exceeded 10. Termination of calculation of the ray. Proable cause: too large integration step along a ray, DT, see input data line 13. (43) - the ray terminates at the vertical plane perpendicular to the profile source-borehole in the VSP mode. (50) - the source is situated outside the model, the computation of the ray diagram does not start. (IND.GT.100) - termination of the ray on the (IND-100)th interface. IND1: ABS(IND1) - number of the interface last hit by the ray. Negative sign of IND1 indicates at least one compound segment of the ray, see description of waves and wave codes. ITER: Number of iterations necessary to find the ray arriving to the prescribed receiver. II: Dummy parameter. INDEX: Current number of the ray in the system of rays for the given code. UU: Square root of the product of the ratios of impedances of the reflected/transmitted to the incident waves at individual points of incidence along the ray path, including the ratios of cosines of angles of reflection/ transmission to the cosines of angles of incidence at the same points. APX(J,I),APY(J,I),APZ(J,I)... for termination point in an anisotropic medium or for P wave arrival to the termination point situated in an isotropic medium, APX(1,I),APY(1,I),APZ(1,I) are complex amplitudes of the x,y and z components of the polarization vector of the considered wave in the "model" coordinates; for an S wave arrival to termination point situated in an isotropic medium, APX(J,I),APY(J,I),APZ(J,I) are complex amplitudes of the x,y and z components corresponding to the two (J=1,2) unit vectors in the plane perpendicular to the ray, describing the polarization of S waves, see the description of the model and the source in the introduction. SPR... geometrical spreading. KMAH... index specifying how many times the considered ray touched a caustic. The information about boundary or critical rays determined in iterations along the profile is stored in one line of the form: R,ZZ,T,DEC,IND,IND1,IRES 'BOUNDARY (CRITICAL) RAY' The meaning of the individual parameters is the same as above. IRES: Indicates whether the boundary ray terminating on the considered profile was found (IRES=1) or not (IRES=0). MOUT=2: This option has meaning only in the two-point ray-tracing mode. For each code the results concerning all of the rays ariving to the considered profile are stored. All the results are stored in the routine RECEIV. The results for the rays specified by declinations from the basic system of declinations (see line 14 in input data file LIN) are stored in one line as follows: IND,IND1,KMAH,R,T,DEC. The results for the rays in the iteration to a specified receiver are stored in one line as follows: 'ITERATIONS' IND,IND1,ITER,KMAH,DD,R,T,DEC,AZIM. The results for the rays iterating to the boundary or critical rays are stored in one line as follows: 'BOUNDARY (CRITICAL) RAY' IND,IND1,ITER,R,T,DEC. The results for the rays arriving to prescribed receivers, for boundary and critical rays are stored in the same form as for MOUT=1. In addition to that information the results of the dynamic ray tracing are stored in the form of two 3x3 matrices XDYN YDYN The first matrix contains quantities related to derivatives of the Cartesian coordinates with respect to the ray coordinates, the second one quantities related to derivatives of components of the slowness vector with respect to the ray coordinates. In fact, XDYN and YDYN are obtained from the dynamic ray tracing equations, see Appendix A of [5], with the initial conditions (A.5), at the same reference. The meaning of other symbols in the above results is the same as explained for the parameter MOUT=1. MOUT=3: Similar information as above is stored for all rays in the iterative procedure. The results stored in the routine PROFIL are marked by asterisk * in the beginning of line to distinguish them from results stored in the routine RECEIV. The use of the option MOUT=3 is not recommended. Description of COMMON block/RAY/ -------------------------------- COMMON/RAY/ AY(28,2000),DS(20,20),KINT(20),HHH(3,3),TMAX,PS(3,7,20), IS(8,20),DINC,DTOLD,N,IREF,IND,IND1 AY(I,N)... values of quantities along the ray. N: successive number of a point on the ray. I=1: travel time I=2-4: x,y and z coordinates of the point of the ray. I=5-7: x,y and z components of the slowness vector. IANI.NE.0 I=8-28: density-normalized elastic parameters of an anisotropic medium IANI=0 I=8-11: P wave velocity and its first derivatives with respect to x, y and z coordinates. I=12-15: S wave velocity and its first derivatives with respect to x, y and z coordinates; I=16: the density I=17-19: components of the e_1 vector I=20-22: components of the e_2 vector I=23-28 unspecified. DS(I,IREF)... values of quantities at the point of incidence at an interface. IREF: successive number of the point of incidence. I=1-3: x,y and z components of the unit normal to the interface. I=4-6: x,y and z components of the projection of the slowness vector of the incident wave into the plane tangent to the interface. I=7: positive value of the projection of the group velocity vector of the incident wave to the normal to the interface. I=8: value of density on the side of incident wave. I=9: no meaning. I=10: positive value of the projection of the group velocity vector of the generated wave to the normal to the interface. I=11: value of density on the side of generated wave. I=12: no meaning. KINT(N)... successive numbers of points along the ray at points of incidence. N: successive number of the element of the ray. For compound element of the ray KINT=0. HHH(I,J)... I-th component of the J-th vector of the ray-centered vector basis e1, e2, t. TMAX... maximum value of travel time; if it is exceeded, the calculation of the ray terminates. PS(I,J,K)... slowness vectors of the incident and generated waves at points of incidence at interfaces; I: index specifying the x-, y- and z-component of the slowness vector (I=1,2,3). J: index specifying the type of the considered wave; J=1-3: reflected wave, J=4-6: transmitted wave; J=7: incident wave. K: successive number of the point of incidence. IS(I,J)... values of quantities related to the J-th point of incidence at an interface. J: successive number of the point of incidence. IS(1,J): number of the layer, in which generated wave propagates. IS(2,J): specifies behaviour of the wave at an interface, specified in TRANSL, see variable ITRANS IS(2,J)=0: reflection, IS(2,J)=1: transmission, IS(2,J)=999: surface conversion. IS(3,J): number of the layer in which incident wave propagates. IS(4,J): dummy parameter IS(5,J): number of inhomogeneous waves generated in the layer where reflected wave propagates. IS(6,J): number of inhomogeneous waves generated in the layer where transmitted wave propagates. IS(7,J): number of the layer in which reflected wave would propagate. IS(8,J): number of the layer in which transmitted wave would propagate. N... curent number of a point along the ray. IREF... curent number of the point of incidence. IND... parameter describing the history of the ray, see the description of the output to the file LOU. IND1... ABS(IND1) - number of the last interface hit by the ray on its way to the receiver, see the description of the output to the file LOU. References ---------- [1] Gajewski,D., Psencik,I., 1986. Numerical modelling of seismic wave fields in 3-D laterally varying layered anisotropic structures - Program ANRAY86, Internal Report Inst.of Earth and Planet.Phys., University of Alberta, Edmonton. [2] Gajewski,D., Psencik,I.,1989. Ray synthetic seismograms in 3-D laterally inhomogeneous anisotropic structures - Program ANRAY89, Internal Report Centre for Computational Seismology, LBL, Berkeley. [3] Cerveny,V., Psencik,I.,1984. Documentation of earthquake algorithms. SEIS83 - Numerical modeling of seismic wave fields in 2-D laterally varying layered structures by the ray method. E.R.Engdahl, edit., Report SE-35, Boulder, 36-40. [4] Cline,A.K.,1981. FITPACK, Dept. of Computer Sciences, Univ. of Texas at Austin. [5] Psencik,I., Teles,T.N., 1996. Point source radiation in inhomogeneous anisotropic structures, PAGEOPH, 148, 591-623.