Program ANRAY - Version 4.00 **************************** 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 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; - 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 elastic parameters specified in grid points of a 3-D grid; - 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; - 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. 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 elastic parameter 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) "Observer" Cartesian coordinate system This system is used for the display of polarization vectors at the model surface. The system is also right-handed. Its y and z axes are oriented in the opposite sense when compared with the "model" coordinate system, i.e. the positive y axis in this coordinate system corresponds to the negative y axis in the "model" system and the vertical (z) axis is positive upwards. C) "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. D) "Elastic tensor" Cartesian coordinate system A convenient and frequently used coordinate system for the description of the anisotropic properties of a material is the system whose axes coincide with the principal axes of the elasticity tensor. Such a system is very practical since in it many of the 21 elastic parameters are zero. In this program the data related to the elastic parameters may be specified in such a coordinate system. Additional information is necessary specifying how the "elasticity tensor" coordinate system is to be rotated to coincide with the "model" coordinate system, see description of input data sub 7a. 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 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 elastic parameters are taken to have the dimensions of velocity squared, which corresponds to the elastic parameters normally associated with an anisotropic medium divided by the density of the medium. The P and S wave velocities in an isotropic layer are, for simplicity, also referred to as elastic parameters. Two methods of the determination of the elastic parameters in individual layers may be used. (a) Linear interpolation of elastic parameters along vertical lines between interfaces representing isosurfaces of elastic parameters. In this approach, the distribution of elastic parameters is laterally varying within the layers if their boundaries are curved or inclined. The vertical gradient of the 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 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 or m and m/s must be used consistently. The density should be always specified in g/cm**3. The point source may be situated at any point of the model. 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 positive component into the z axis in the "observer" coordinates (where the z axis is positive upwards). 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 can be speeded up by using paraxial ray approximation. Description of input and output files ------------------------------------- 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. Output data for plotting ray diagrams, travel time curves and amplitude-distance curves in the program ANRAYPL are stored in the file LU1. Output data for computing ray synthetic seismograms in the program SYNTAN are stored in the file LU2. Specification of the files LIN, LOU, LU1 andd LU2 can be made through the routine SERV, which is part of this program package. If the routine SERV is not used (put 'C' in front of 'CALL SERV' in the beginning of the main program), the files LIN and LOU are automatically specified LIN=5, LOU=6, and the numbers LU1 and LU2 are read in from the file LIN, see below. Input data in the file LIN -------------------------- 1) MTEXT FORMAT(20A4) MTEXT... description of the problem under consideration, (alphanumeric text). 2) LU1,LU2,INULL,IPLOT FORMAT(17I5) LU1... the number of the output file with the data for plotting in the program ANRAYPL. It has no meaning if the routine SERV is used. LU2... the number of the output file with the data for construction of synthetic seismograms in the program SYNTAN. It has no meaning if the routine SERV is used. 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 value (INULL=4). IPLOT... the switch deciding about the mode in which is the program ANRAY going to work. IPLOT=0... the program ANRAY is used for computing rays, travel times and ray amplitudes. IPLOT=1... the program ANRAY is used for generating file for plotting plane sections of slowness, phase velocity or group velocity surfaces. 3) MPRINT,NINT FORMAT(14I5) 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. 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 elastic parameters. NINT... number of interfaces in the model (including first interface representing the model surface and the last interface representing the model bottom). 4) The lines 4a-4d repeat NINT times, each for one interface, from the top to the bottom of the model. 4a)MX,MY FORMAT(14I5) MX,MY... number of grid lines x=const and y=const for approximation of the interfaces. 4b)SX(1),...,SX(MX) FORMAT(8F10.5) SX(I)... x-coordinates of grid lines for approximation of the interfaces. 4c)SY(1),...SY(MY) FORMAT(8F10.5) SY(I)... y-coordinates of grid lines for approximation of the interfaces. 4d)Z(1),...,Z(MX*MY) FORMAT(8F10.5) 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 FORMAT(8F10.5) 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. 6) Cards controlling the approximation of elastic parameters and density distribution 6a) ISQRT,IRHO FORMAT(14I5) ISQRT... controls the approximation of elastic parameters. ISQRT=0: Approximation in elastic parameters. ISQRT=1: Approximation in square roots of 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 elastic parameter A11 (which in an isotropic layer corresponds to the P-wave velocity). 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) FORMAT(8F10.5) RHO(I)... density in the I-th layer. It has meaning only if IRHO=1. If RHO(I) is not specified, RHO(I)=1.0 by default. 7) Description of distribution elastic parameters in individual layers. The input data differ for isosurface interpolation and B-spline approximation. ********************** B-spline approximation ********************** The 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 FORMAT(I10,F10.5) IANI... specifies properties of the layer. IANI=0: The layer is isotropic. IANI=1: The layer is anisotropic. 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. 7b)NX,NY,NZ FORMAT(26I3) 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. 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. 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. 7c)X1(1),X1(2),...,X1(NX) FORMAT(8F10.5) X1(I)... coordinate of the I-th grid line in the x-direction. 7c)Y1(1),Y1(2),...,Y1(NY) FORMAT(8F10.5) Y1(I)... coordinate of the I-th grid line in the y-direction. 7c)Z1(1),Z1(2),...,Z1(NZ) FORMAT(8F10.5) Z1(I)... coordinate of the I-th grid line in the z-direction. 7d) Values of elastic parameters at the grid points of the 3-D grid. In an isotropic layer, when IANI=0, values of the P and S wave velocities are read in successively. In case of an anisotropic layer, individual elastic parameters are read in successively. The elastic parameters are ordered in the following way: A11,A12,A13,A14,A15,A16, A22,A23,A24,A25,A26,A33,A34,A35,A36,A44,A$%,A46,A55,A56,A66. Each elastic parameter is read in in the following way (((W(I,J,K),K=1,NZ),J=1,NY),I=1,NX) FORMAT(8F10.5) W(I,J,K) contains the value of II-th elastic parameter at the grid point (X(I),Y(J),Z(K)). Note: The maximum numbers of grid lines to specify the distribution of 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. 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 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) FORMAT(I10,6F10.5) IANI... specifies properties of the layer. IANI=0: The layer is isotropic. IANI=1: The layer is anisotropic. ANGU(1),...,ANGU(3)... the rotation angles, by which the "elasticity tensor" coordinate system is to be rotated into the "model" coordinate system, specified at the upper (U) interface bounding the layer. The angles are specified in degrees. ANGU(1) is an angle between positive direction of the "model" x-axis and horizontal projection of positive part of "crystal" z-axis. ANGU(1) is positive if measured from positive direction of "model" x-axis towards positive "model" y-axis. ANGU(2) is an angle between the positive directions of the "model" z-axis and "crystal" z-axis. ANGU(3) is a rotation angle in the plane perpendicular to the "crystal" z-axis. For ANGU(3)=0., "crystal" y-axis is horizontal, x-axis is perpendicular to it and points down. For ANGU(3) non-zero, the "crystal" x- and y-axes are rotated positvely from their positions for ANGU(3)=0. ANGL(1),...,ANGL(3)... the same as for ANGU, but at the lower (L) interface bounding the layer. Note 1: The values of rotation angles inside of a layer are determined in the same manner as the 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) FORMAT(6F10.5) A66U... elastic parameters (i.e. actual parameters divided by density) 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) FORMAT(6F10.5) 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) FORMAT(6F10.5) A66L... elastic parameters (i.e. actual parameters divided by density) 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) FORMAT(6F10.5) 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 IPLOT=1 (see line 2). NPAR FORMAT(16I5) NPAR... switch indicating a section which wave surface is to be plotted. NPAR=1: slowness surface. NPAR=2: phase velocity surface. NPAR=3: group velocity surface. 9) This line is read only if IPLOT=1 (see line 2). X0,Y0,Z0,DDELTA FORMAT(4F10.5) 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. DDELTA... step in the angle of the phase normal 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. 10) ICONT,MEP,MOUT,MDIM,METHOD,MREG,ITMAX,IPOL,IPREC,IRAYPL, IPRINT,IAMP,MTRNS,ICOEF,IRT,ILOC,MCOD,MORI FORMAT(24I3) 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. 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. 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.0: 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. 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. 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. 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. Note 2: Paraxial ray approximation can be used only if results of dynamic ray tracing are available. Therefore, for MEP.NE.0, MDIM is set automatically equal to 1, MDIM=1. MREG... controls the evaluation of amplitudes at the termination points of the rays. MREG=0: The conversion coefficients are applied if the termination point is situated on the surface. MREG=1: The conversion coefficients are not applied. Note: MREG is set to be equal to 1 automatically for VSP computations and for receivers situated along internal interfaces (when ILOC.gt.0). In these cases conversion coefficients are not considered. 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. 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. 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. 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. 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. 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. 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 FORMAT(8F10.5) 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. 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 FORMAT(8F10.5) 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). For MEP.GT.1, NDST=MEP. PROF,RMIN,RSTEP,XPRF,YPRF FORMAT(8F10.5) 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. 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 FORMAT(8F10.5) XVSP,YVSP... x and y coordinates of the vertical profile in case of VSP computations. 13) XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,REPS,PREPS FORMAT(8F10.5) XSOUR,YSOUR,ZSOUR... (x,y,z) coordinates of the source. TSOUR... initial time. 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. 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. 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 FORMAT(8F10.5) 15) This line is read only when MCOD.EQ.0 BMIN,BSTEP,BMAX FORMAT(8F10.5) AMIN,ASTEP,AMAX,BMIN,BSTEP,BMAX... define the basic system of initial angles (in radians) for ray tracing. The 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) FORMAT(26I3) KC... KC.EQ.0: Ray of direct unconverted wave with no reflections and no conversions at interfaces. 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. 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 and 2 give the same 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 FORMAT(8F10.5) BMIN,BSTEP,BMAX FORMAT(8F10.5) The meaning of the above parameters is the same as in lines 12 and 13. Termination of computations --------------------------- The sequence of lines 16-17 can be repeated any number times, separately for each considered elementary wave. 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 (if LU1.NE.0.AND.MEP.NE.0) for plotting in program ANRAYPL, 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 No.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) VP,VS,RO... P-wave velocity, S-wave velocity and 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... dummy parameter. Note: Maximum value of N is 400. 8) (AY(2,I),AY(3,I),AY(4,I),I=1,N) FORMAT(3E15.5) The (x,y,z) coordinates of the points along the rays, starting from the source (I=1) and ending with the termination point (I=N). (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. 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. 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 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 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. (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 400 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 "observer" 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 derivatives of the Cartesian coordinates with respect to the ray coordinates, the second one the derivatives of the components of the slowness vector with respect to the ray coordinates. 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,400),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: 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-28 unspecified. DS(I,N)... values of quantities at the point of incidence at an interface. N: 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 n-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 incident 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. 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.