Program WEAKAN ************** Program WEAKAN is designed for the calculations in the quasi- isotropic (QI) approximation. This version calculates QI ray amplitudes of the qS waves propagating in continuous, inhomogeneous, weakly anisotropic media. It solves two coupled ordinary differential equations providing zero-order QI amplitudes of qS waves. The equations are solved along a ray in a background isotropic medium. The program WEAKAN uses the file LU1 generated by the program package ANRAY. WEAKAN generates file LU6, which can be used in the program FRESAN and following programs for computation of synthetic seismograms. The program yields regular results in regions, in which the program ANRAY, based on the standard ray method for anisotropic media fails, i.e. in weakly anisotropic media and in singular regions like cross and kiss singularity. It fails in vicinity of caustics. 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' 'LU6'/ 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 input data file LU1, generated by the program ANRAY. 'LU6' is the name of the output data file LU6. / is a slash recomended in batch and script files to enable future extensions. Defaults: 'LIN'='anray.dat' 'LOU'='anray.out' 'LU1'='lu1.dat' 'LU6'='lu6.dat' Example of the main input data: 'weakan.qi' / Input data consist partially of the data generated by the program ANRAY, stored in the formatted form in the file LU1, and partially of the additional input data read by list-directed input (free format), prepared by the user and stored in the file LIN. Output data describing the computations are stored in the file LOU. Output data for computing ray synthetic seismograms in the program FRESAN are stored in the file LU6. The data stored at LU1 ---------------------- The data are stored in the file LU1 in the formatted form. For details see the description of the content of the file LU1 in program ANRAY. 1) ICONT,NDST,ILOC FORMAT(26I3) 2) ROS FORMAT(8F10.5) 3) NPN,NPN,NPN FORMAT(26I3) 4) APN,APN,APN,APN,APN FORMAT(5E15.5) 5) XPRF,YPRF,0.,PROF FORMAT(8F10.5) 6) (DST(I),I=1,NDST) FORMAT(8F10.5) 7) NTOT,II FORMAT(26I3) 8) (T(J),X(J),Y(J),Z(J),P1(J),P2(J),P3(J),VP(J),VS(J),RHO(J), (E(1,K,J),K=1,3),(E(2,L,J),L=1,3),J = 1,N) FORMAT(16E15.5) 9) INDEX FORMAT(26I3) 10) (INDI(I),XCOOR(I),ZCOOR(I),TIME(I),AMPX(1,I),AMPY(1,I), AMPZ(1,I),I=1,INDEX) FORMAT(I5,9E15.5) 11) IF(INDEX.LT.0) ((AMPX(2,I),AMPY(2,I),AMPZ(2,I),I=1,INDEX) FORMAT(6E15.5) 12) (P(K),K=1,3) FORMAT(3E15.5) 13) (POL(K),K=1,3) FORMAT(3E15.5) 14) IF(INDEX.LT.0) (POL1(K),K=1,3) FORMAT(3E15.5) The sequence of data 1-14 repeats for each considered elementary wave. Indication of the end of information in LU1 is ICONT=0 in line 1. The additional input data in the file LIN ----------------------------------------- 1) MTEXT MTEXT... description of the problem under consideration, (alphanumeric text). Default 'WEAKAN'. 2) INULL,INEWB 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. INEWB... switch deciding if the QI approach is to be used with the standard weak anisotropy matrix B_{mn} or with its modified version calculated in the routine NEWB INEWB=0: standard QI approach. INEWB=1: modified weak anisotropy matrix used. Default INEWB=0. 3) MPRINT,NINT 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 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 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. 6) Cards controlling the approximation of elastic parameters and density distribution 6a) ISQRT,IRHO ISQRT... controls the approximation of elastic parameters. ISQRT=0: Approximation in elastic parameters. Default value. 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). 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 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 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. 7c)Y1(1),Y1(2),...,Y1(NY) Y1(I)... coordinate of the I-th grid line in the y-direction. 7c)Z1(1),Z1(2),...,Z1(NZ) 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 squares 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, A45,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) 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. 8) F(1),F(2),F(3) F(I)..I-th component of the force vector. 9) FL,FD,NF FL... lower limit of the considered frequency range. FD... frequency sampling step. NF... number of frequency samples (NF may not exceed 1000). Example of data LIN for isotropic background for QI computations Termination of computations --------------------------- The computation terminates when NTOT=0 is found on 7th line of LU1. Output to the file LU6 ---------------------- Data necessary for calculation of the frequency response in the program FRESAN are stored in the file LU6. The data are stored in LU6 in the following formatted form: 1) FL,FD,NF FORMAT(2F10.5,I5) FL... lower limit of the considered frequency range. FD... frequency sampling step. NF... number of frequency samples. 2) (E(1,K,NTOT),K=1,3),(E(2,L,NTOT),L=1,3) FORMAT(16E15.5) E(I,J,N)... J-th component of the I-th ray-centered basis vector at the termination point of the ray (NTOT). 3) FF,W FORMAT(16E15.5) FF... frequency W(I)... (I=1,2) complex-valued I-th ray amplitude in the QI approximation for the frequency FF. The line 3 repeats NF-times for each ray read from LU1. The lines 2 and 3 repeat for each ray. Output to the file LOU ---------------------- The file LOU contains a copy of the input data. Additional data stored in LOU is information about the model. It is controlled by the switch MPRINT. 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.