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:
'weak.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.