C Subroutine file 'init.for' to read the input data for the initial C surface, and to define initial values for complete ray tracing. C C By Vlastislav Cerveny, Ludek Klimes, Ivan Psencik C C This file consists of the following subroutines: C INITB...Block data subroutine defining common block /INITC/ to C store the important quantities at the initial point of the C ray, and common block /INISC/ to store some of the input C data. C INIT1...Interface subroutine to INIS1 reading the input data for C the four functions specifying the initial conditions for C the computed rays. C INIT2...Subroutine evaluating, for given ray take-off parameters, C the values of the computed quantities at the initial point C of the ray, and storing the important quantities at the C initial point of the ray in the common block /INITC/. C INIS1...Sample subroutine designed to read the input data for the C initial points of rays. A two-parametric system of rays C of each elementary wave is assumed. A ray of the C elementary wave is specified by its two take-off C parameters. The computed rays may start from a single C initial point common to all rays, from a curve along which C an initial travel time is defined, from an initial surface C along which an initial travel time is defined, etc. C INIS2...Sample subroutine returning the functional values and C their first and second derivatives, of the functions C describing the initial surface. C SQRT3...Subroutine evaluating the square root of the given C real-valued positive-semidefinite symmetric 3*3 matrix. C SPHERE..Subroutine transforming spherical coordinates PAR1, PAR2 C into the Cartesian coordinates of the corresponding point C on the unit sphere. It also evaluates the first and C second derivatives of the Cartesian coordinates with C respect to PAR1 and PAR2. C Subroutines INIS1 and INIS2, defining the common initial point, C initial curve or initial surface, call subroutines VAL1 and VAL2 which C must be appended. In addition, subroutines CURVN1 (or its alternative C CURVB1), CURV2D (or its alternative CURVBD), SURFB1, SURFBD, VAL3B1, C VAL3BD, VGEN, TERMS, SNHCSH, TRIDEC, TRISOL, DSPLNZ, INTRVL from the C subroutine package 'FITPACK' by Alan Kaylor Cline, Department of C Computer Sciences, University of Texas at Austin, are used. In the C complete ray tracing, subroutines INIS1 and INIS2 may be replaced by C any user-defined package containing subroutines INIS1 and INIS2 with C the same number, type and meaning of their parameters as in this C file. C C Input data for the initial points of rays: C The data are read in by the list directed input (free format). In C the list of input data below, each numbered paragraph indicates C the beginning of a new input operation (new READ statement). If C the first letter of the symbolic name of the input variable is C I-N, the corresponding value in input data must be of the type C INTEGER. Otherwise (except TEXTI), the input parameter is of the C type real. C (1) TEXTI C String describing the data. Only the first 80 characters of the C string are significant. C (2) INIDIM,INIPAR,ADVANC,/ C Quantities defining the kind of initial conditions (the kind of C the source). C INIDIM..Determines the dimensionality of the source: C -1...Single initial point (point source), its coordinates C are read from a separate file. C 0...Single initial point (point source), C 1...Initial curve (line source), C 2...Initial surface. C INIPAR..Determines the parametrization of rays: C For INIDIM.LE.0: C INIPAR.LT.0: The same as for IABS(INIPAR), but with C the unit vector (T1,T2,T3) tangent to the ray C changed to (T1,-T2,-T3). C INIPAR.EQ.1: Ray parameters are polar-like spherical C coordinates (colatitude,longitude) connected with the C local Cartesian coordinate system which basis vectors C are given by the square root of the metric tensor at C the initial point. C If SIN(colatitude).LT.0, the ray is reported out of C the ray-parameter domain: IEND=71 in subroutine INIT2. C INIPAR.EQ.2: Ray parameters are geographic-like C spherical coordinates (longitude,latitude) connected C with the local Cartesian coordinate system which basis C vectors are given by the square root of the metric C tensor at the initial point. C If COS(latitude).LT.0, the ray is reported out of C the ray-parameter domain: IEND=71 in subroutine INIT2. C INIPAR.EQ.3: Azimuthal equidistant projection of a unit C sphere is parametrized by 2 Cartesian coordinates C centred at the projection point. This option is C suitable especially for reflecton seismic studies. C The unit vector tangent to the ray, expressed in the C local Cartesian coordinate system which basis vectors C are given by the square root of the metric tensor at C the initial point, is given by C T1=PAR1*SIN(R)/R C T2=PAR2*SIN(R)/R C T3= COS(R) C with R=SQRT(PAR1*PAR1+PAR2*PAR2). C If R.GT.2*PI, the ray is reported out of the C ray-parameter domain: IEND=71 in subroutine INIT2. C For INIDIM=1: C INIPAR must be 1 or 2. The INIPAR-th ray parameter is C identical with the parameter parametrizing the initial C curve. The other ray parameter is the angle between the C ray take-off plane and the normal to the interpolated C surface. The ray take-off plane is given by the tangent C to the initial line and by the slowness vector. C For INIPAR=1, the initial line is the line PAR2=0 at the C interpolated surface and is parametrized by PAR1. C For INIPAR=2, the initial line is the line PAR1=0 at the C interpolated surface and is parametrized by PAR2. C For INIDIM=2: C Ray parameters are identical with two parameters C parametrizing the initial surface. C INIPAR.LE.0: Initial surface is described in terms of C functions specifying the dependence of general C coordinates (X1,X2,X3) on two parameters of the C initial surface. C INIPAR.EQ.1: Initial surface is specified in the C polar-like spherical coordinates (colatitude, C longitude, radius) connected with the local Cartesian C coordinate system which basis vectors are given by the C square root of the metric tensor at the given point. C Colatitude and longitude are the parameters, and the C initial surface is determined by a function specifying C the dependence of the radius on these parameters C (colatitude and longitude). C INIPAR.GE.2: Initial surface is specified in the C geographic-like spherical coordinates (longitude, C latitude, radius) connected with the local Cartesian C coordinate system which basis vectors are given by the C square root of the metric tensor at the given point. C Longitude and latitude are the parameters, and the C initial surface is determined by a function specifying C the dependence of the radius on these parameters C (longitude and latitude). C ADVANC..Initial point of the ray is shifted by distance ADVANC C perpendicularly to the initial surface or line, C or tangentially to the ray for the single initial point. C All initial and other quantities (except for the metric C tensor) are then evaluated at the shifted initial point. C Finally, the initial travel time is linearly updated C using the initial slowness vector. This option may be C usefull if the source is situated close to a structural C interface. C Default: C INIDIM=-1, INIPAR=2, ADVANC=0. C (3) Data describing the initial point, curve or surface. C For INIDIM=-1: C (3A) 'FSRC' C 'FSRC'... Name of the input file containing the C coordinates of a single initial point and the initial C value of the travel time. C For INIDIM=0: C (3A) X1INI,X2INI,X3INI,TTINI C X1INI,X2INI,X3INI... Coordinates of a single initial C point. C TTINI... Initial value of the travel time. C For INIDIM=1: C (3B) Input data for NFUNC=4 functions X1(.,.),X2(.,.), C X3(.,.), TT(.,.) of two variables. The INIPAR-th one of C the 2 variables being simultaneously the ray parameter). C For INIDIM=2, INIPAR.LE.0: C (3B) Input data for NFUNC=4 functions X1(.,.),X2(.,.), C X3(.,.), TT(.,.) of two variables (two initial surface C parameters, being simultaneously the ray parameters). C For INIDIM=2, INIPAR.GE.1: the following two data sets (3A), (3B): C (3A) X1INI,X2INI,X3INI... Coordinates of a given point, C see the description of the input data (1). C (3B) Input data for NFUNC=2 functions R(.,.),TT(.,.) of C two variables (two initial surface parameters, being C simultaneously the ray parameters). C R(.,.) describes the radius, see input data (1), C TT(.,.) is the initial travel time. C Default: 'FSRC'='src.dat', TTINI=0. C The structure of the input data (3B) is given by the subroutine VAL1 C and is described later on. C C Input file 'FSRC' containing the coordinates of the initial point: C This file is used only if INIDIM=-1, see the input data (2) above. C (1) Several strings terminated by / (a slash). C The simlest way is to submit just the /. C (2) 'NAME',X1INI,X2INI,X3INI,TTINI,/ C 'NAME'... String reserved for the name of the initial point. C No meaning here, anything in apostrophes may be submitted. C X1INI,X2INI,X3INI... Coordinates of a single initial point. C TTINI... Initial value of the travel time. C Default: TTINI=0. C C Storage in the memory: C The input data (2) and (3A) are stored in the common block C /INISC/. The important quantities at the initial point of the ray C (see C.R.T.6.1) are stored in the common block /INITC/. These C common blocks are defined in the following subroutine: C ------------------------------------------------------------------ BLOCK DATA INITB INCLUDE 'initd.inc' INCLUDE 'initc.inc' END C ------------------------------------------------------------------ C C Above mentioned input data (3B) for the initial curve or for the C initial surface are read in by the subroutine VAL1 and have the C following structure: C These input data define at least NFUNC individual functions C describing the initial conditions. They are read in by subroutine C VAL1 called by INIS1. The number MFUNC of all functions specified C in the input data may be greater or equal to NFUNC. The data are C read in by the list directed input (free format). C (1) MFUNC C The number of all input functions. It must be greater or equal to C the number NFUNC of the functions required to describe the C coordinates and travel time along the initial curve or surface. C The functions indexed 1 to NFUNC must be the functions describing C the coordinates and travel time along the initial curve or C surface. C (2) NFUNC-times (i.e. once for each function) input data (2A)+(2B): C (2A) TEXTF,IFUNC C Identification of the function. C TEXTF...Any string. Its first 3 characters must differ from 'END'. C IFUNC...Index of the function: C 1 to 3... coordinates and 4... travel time, or C 1... radius and 2... travel time. Amplitude and/or other C quantities may follow. C (2B) 'Input data for one function', see below. C (3) TEXTE,AUX C End of data. C TEXTE...String, the first 3 characters of which must be upper-case C 'END'. C AUX... Any number or a slash. C Remark: C If the initial surface coincides with a structural interface C (e.g., exploding-reflector initial conditions), it is reasonable C to slightly shift the initial surface from the structural C interface into the complex block into which the wave propagates. C Otherwise, because of rounding errors, there is a danger that some C parts of the initial surface are situated at the opposite side C of the interface. C C Input data for one function: C The data are read in by the list directed input (free format). In C the list of input data below, each numbered paragraph indicates C the beginning of a new input operation (new READ statement). If C the first letter of the symbolic name of the input variable is C I-N, the corresponding value in input data must be of the type C INTEGER. Otherwise, the input parameter is of the type REAL. C (1) IVAR1,IVAR2,0,SIGMA C The form of the function. C IVAR1,IVAR2... Denote the form of the function. The function must C be of the form C F(G1,G2) = W(A1,A2)-B1-B2 . C G1, G2 are the ray parameters. Each of A1, A2, B1, B2 must C be either: (a) one of ray parameters G1, G2, or (b) must C be left out. At most 2 of parameters A1-B2 may be of kind C (a). Note that IVAR1 controls the type of A1 and B1, C IVAR2 controls the type of A2 and B2. C For IVAR1.EQ.0: A1, B1 are empty (left out), C for IVAR1.EQ.1: A1=G1, B1 is empty, C for IVAR1.EQ.2: A1=G2, B1 is empty, C for IVAR1.EQ.-1: B1=G1, A1 is empty, C for IVAR1.EQ.-2: B1=G2, A1 is empty, C the meaning of the parameter IVAR2 is similar. C Examples: C IVAR1: IVAR2: the form of the function: C 1 2 F(G1,G2)=W(G1,G2) C 2 1 F(G1,G2)=W(G2,G1) C 1 0 F(G1,G2)=W(G1) C 1 -2 F(G1,G2)=W(G1)-G2 C Function W is interpolated by means of splines under C tension. C SIGMA...Is the tension factor (its sign is ignored). This value C indicates the curviness desired. If ABS(SIGMA) is nearly C zero (e.g. 0.001), the resulting surface is approximately C the tensor product of cubic splines. If ABS(SIGMA) is C large (e.g. 50.), the resulting surface is approximately C tri-linear. If SIGMA equals zero, tensor products of C cubic splines result. A recommended value for SIGMA is C approximately 1. In absolute value. C (2) NX(1),...,NX(NVAR) C The numbers of grid coordinates for the interpolation. C This input is performed if at least one of IVAR1, IVAR2 is C positive. C Each of NX(1),...,NX(NVAR) corresponds to one positive value of C IVAR1, IVAR2 and specifies the number of grid coordinates C corresponding to that independent variable of function W, see (1). C The sign of NX(1),...,NX(NVAR) is ignored. NVAR (.LE.2) is the C number of positive values of the above quantities IVAR1, IVAR2, C i.e. the number of independent variables of function W, see (1). C (3) X1(1),...,X1(NX(1)) C The grid coordinates corresponding to the first independent C variable of function W, see (1). C This input is performed if NX(1) is specified, see (2), and is not C zero. The grid coordinates may be specified in any order. C (4) X2(1),...,X2(NX(2)) C The grid coordinates corresponding to the second independent C variable of function W, see (1). C This input is performed if NX(2) is specified, see (2), and is not C zero. The grid coordinates may be specified in any order. C (5) (((W(I,J),I=1,MAX(NX(1),1)),J=1,MAX(NX(2),1)) C The values of function W at grid points. Function value W(I,J) C corresponds to point (X1(I),X2(J)). C C Date: 1996, May 9 C Coded by Ludek Klimes C C======================================================================= C SUBROUTINE INIT1(LUN) INTEGER LUN C C Subroutine INIT1 calls the subroutine INIS1 to read the input data for C the initial points of rays. C C Input: C LUN... Logical unit number of the external input device C containing the input data. C The input parameter is not altered. C C No output. C C Common block /INITC/: INCLUDE 'initc.inc' C None of the storage locations of the common block, except ICB1I, are C altered. ICB1I is set to zero. C C Subroutines and external functions required: EXTERNAL INITB,INIS1 C INITB.. Block data subroutine of this file. C INIS1... This file. C C Date: 1993, December 18 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER NFUNC PARAMETER (NFUNC=4) C CALL INIS1(LUN,NFUNC) ICB1I=0 RETURN END C C======================================================================= C SUBROUTINE INIT2(PAR1,PAR2,YL,Y,YY,IY,IEND,IWAVE0,IKODE) REAL PAR1,PAR2,YL(6),Y(35),YY(5) INTEGER IY(12),IEND,IWAVE0,IKODE C C Subroutine INIT2 evaluates, for given ray take-off parameters, the C values of the computed quantities at the initial point of the ray, and C stores the important quantities at the initial point of the ray in the C common block /INITC/. C C Input: C PAR1,PAR2... Ray take-off parameters. C YL,Y,YY,IY... Arrays of dimensions at least 6, 35, 5, 12, C respectively. C IWAVE0..Index of the already computed elementary wave having the C most numerous common elements with the current elementary C wave. Need not be defined if IKODE=0. C IKODE...The length of the common part of the codes of the IWAVE-th C and IWAVE0-th elementary waves. C The input parameters PAR1, PAR2 are not altered. C C Output. C YL... Array containing local quantities at the initial point of C the ray (see C.R.T.5.5.4). The quantities are listed in C the subroutine file 'ray.for'. C Y... Array containing basic quantities computed along the ray C at the initial point of the ray (see C.R.T.5.2.1). The C quantities are listed in the subroutine file 'ray.for'. C YY... Array containing real auxiliary quantities computed along C the ray at the initial point of the ray (see C.R.T.5.2.2). C The quantities are listed in the subroutine file C 'ray.for'. C IY... Array containing integer auxiliary quantities computed C along the ray at the initial point of the ray (see C C.R.T.5.2.2). The quantities are listed in the subroutine C file 'ray.for'. C IEND... Information on the initial point of the ray: C 0... Computation of the ray may follow. C 71... There is no ray corresponding to the given ray C parameters. E.g., the given parameters do not belong to C the domain of the initial surface. C 72... Initial point of the ray is not situated in the C required complex block. C 73... Initial point of the ray is not situated in the C computational volume. C 74... Ray of the generated wave is not real-valued. C C Common block /DCRT/ (see subroutine file 'ray.for'): INCLUDE 'dcrt.inc' C None of the storage locations of the common block are altered. C C Common block /INISC/: INCLUDE 'initd.inc' C None of the storage locations of the common block are altered. C C Common block /INITC/: INCLUDE 'initc.inc' C All the storage locations of the common block are defined in this C subroutine. ICB1I must be zero before the first invocation of this C subroutine, other storage locations may be undefined. C C Subroutines and external functions required: EXTERNAL NSRFC,BLOCK,VELOC,KOOR,METRIC,PARM2,SMVPRD,CODE,INIS2 INTEGER NSRFC,KOOR C NSRFC,BLOCK,VELOC... File 'model.for' of the package 'MODEL'. C KOOR,METRIC... File 'metric.for' of the package 'MODEL'. C PARM2... File 'parm.for' of the package 'MODEL'. C SMVPRD... File 'means.for' of the package 'MODEL'. C CODE... File 'code.for'. C INIS2... This file. C C Error messages: C 611... Wrong function of subroutine CODE: C Other reason of the termination of the ray computation C than 22 should not be reported by the subroutine CODE when C referenced by the subroutine INIT2. Contact the authors. C 612... Wrong function of subroutine CODE: C Subroutine CODE requires the first ray element to be C situated in another complex block than the initial point. C This error should not appear. Contact the authors. C 613... This program branch is not coded: C Interpolation mode of the complete ray tracing program has C not been enabled yet. See 'ray.for', description of input C data, (2), MODCRT. C 617... Strange common initial point: C Initial conditions determined by subroutine INIS2 do not C resemble common initial point. Contact the authors. C 618... Strange initial line: C Initial conditions determined by subroutine INIS2 do not C resemble initial line. Contact the authors. C 619... Strange initial surface: C Initial conditions determined by subroutine INIS2 do not C resemble initial surface. Contact the authors. C C Date: 1996, September 30 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations for local model parameters: FAUX(10), C G(12),GAMMA(18),GSQRD, UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL: INCLUDE 'auxmod.inc' C C....................................................................... C C Auxiliary storage locations: INTEGER KODIND,ICBNEW,ISB2,ICB2,I,NY,NFUNC PARAMETER (NFUNC=4) REAL E11,F21,F31,C11,B11,FF11,FFE11,CB11,ECB11,U1,T1,BT1 REAL E12,F22,F32,C12,B12,FF12,FFE12,CB12,ECB12,U2,T2,BT2 REAL E22,F23,F33,C22,B22,FF21,FFE21,CB21,ECB21 REAL E(3),F(6,NFUNC), FF22,FFE22,CB22,ECB22 EQUIVALENCE (E(1),E11),(E(2),E12),(E(3),E22) REAL V0,DETC,TRECE,Z1,Z2,Z3,AUX1,AUX2,AUX3 SAVE V0 C C KODIND..Position in the code of elementary wave. C ICBNEW..The index of the complex block in which the initial point C of the ray should be situated, supplemented by the sign C '+' for P wave or '-' for S wave. Output from subroutine C CODE. C ISB2... Index of the simple block in which the initial point is C situated. C ICB2... Index of the complex block in which the initial point is C situated. C I... Auxiliary and loop variable. C NY... Number of defined locations of array Y. C INIDIM..Distinguishes between initial point, line, or surface. C E=(E11,E12,E22)... Projection matrix, see the subroutine INIS2. C F...Functions describing the initial surface and their C derivatives, see the subroutine INIS2. C F21,F22,F23, F31,F32,F33... Covariant components of the vectors C F(2,1),F(2,2),F(2,3), F(3,1),F(3,2),F(3,3) tangent to the C initial surface. C C11,C12,C22... Components of matrix C of space or space-time C scalar products, eventually of its square root. C B11,B12,B22... Components of inverse square root of matrix C. C FF11,FF12,FF21,FF22... Components of the auxiliary matrix FF. C FFE11,FFE21,FFE12,FFE22... Components of the matrix FF*E. C CB11,CB21,CB12,CB22... Components of the matrix 1-C*B. C ECB11,ECB21,ECB12,ECB22... Components of the matrix E*CB/Tr(E*C). C U1,U2... Slowness derivatives with respect to PAR1,PAR2. C T1,T2... Velocity*derivatives of the travel time along the initial C surface with respect to PAR1,PAR2. C BT1,BT2... Components of the vector (B+ECB)*T. C V0... Propagation velocity at the initial point. C DETC... Determinant of matrix C, eventually of its square root. C TRECE... Trace of the matrix E*C*E. C Z1,Z2,Z3... Unit normal to the initial surface - covariant comp. C AUX1,AUX2,AUX3... Auxiliary storage locations. C C....................................................................... C IF(MODCRT.GE.3) THEN PAUSE 'Error 613 in INIT2: This program branch is not coded' STOP END IF C CALL INIS2(NFUNC,PAR1,PAR2,E,F,IEND) IF(IEND.NE.0) THEN RETURN END IF IF(F(1,1).LT.BOUNDR(1).OR.BOUNDR(2).LT.F(1,1).OR. * F(1,2).LT.BOUNDR(3).OR.BOUNDR(4).LT.F(1,2).OR. * F(1,3).LT.BOUNDR(5).OR.BOUNDR(6).LT.F(1,3).OR. * BOUNDR(7).LT.F(1,4)) THEN IEND=73 RETURN END IF IF(E11.EQ.0..AND.E12.EQ.0..AND.E22.EQ.0.) THEN C Initial point: IF(INIDIM.GT.0) THEN PAUSE 'Error 617 in INIT2: Strange common initial point' STOP END IF ELSE IF(E11*E22-E12*E12.LE.0.001) THEN C Initial line: IF(INIDIM.NE.1) THEN PAUSE 'Error 618 in INIT2: Strange initial line' STOP END IF ELSE C Initial surface: IF(INIDIM.NE.2) THEN PAUSE 'Error 619 in INIT2: Strange initial surface' STOP END IF END IF C C Initial travel time YI(1)=F(1,4) C Initial imaginary travel time YI(2)=0. C C Coordinates of the initial point of the ray YI3=YI(3) YI4=YI(4) YI5=YI(5) YI(3)=F(1,1) YI(4)=F(1,2) YI(5)=F(1,3) C C Local coordinate system: C Covariant components of vectors tangent to the initial surface IF(KOOR().NE.0) THEN CALL METRIC(YI(3),GSQRD,G,GAMMA) CALL SMVPRD(G,F(2,1),F(2,2),F(2,3),F21,F22,F23) CALL SMVPRD(G,F(3,1),F(3,2),F(3,3),F31,F32,F33) ELSE GSQRD=1. F21=F(2,1) F22=F(2,2) F23=F(2,3) F31=F(3,1) F32=F(3,2) F33=F(3,3) END IF C Scalar products of vectors tangent to the initial surface C11=F(2,1)*F21+F(2,2)*F22+F(2,3)*F23 C12=F(2,1)*F31+F(2,2)*F32+F(2,3)*F33 C22=F(3,1)*F31+F(3,2)*F32+F(3,3)*F33 DETC=C11*C22-C12*C12 C Unit normal to the initial surface - covariant components AUX1=GSQRD/SQRT(DETC) Z1=(F(2,2)*F(3,3)-F(2,3)*F(3,2))*AUX1 Z2=(F(2,3)*F(3,1)-F(2,1)*F(3,3))*AUX1 Z3=(F(2,1)*F(3,2)-F(2,2)*F(3,1))*AUX1 C C Modification of the coordinates of the initial point of the ray: IF(ADVANC.NE.0.) THEN C Shifting the initial point in the direction of the C unit normal to the initial surface - contravariant components IF(KOOR().NE.0) THEN YI(3)=F(1,1)+ADVANC*(G(1)*Z1+G(2)*Z2+G(4)*Z3) YI(4)=F(1,2)+ADVANC*(G(2)*Z1+G(3)*Z2+G(5)*Z3) YI(5)=F(1,3)+ADVANC*(G(4)*Z1+G(5)*Z2+G(6)*Z3) ELSE YI(3)=F(1,1)+ADVANC*Z1 YI(4)=F(1,2)+ADVANC*Z2 YI(5)=F(1,3)+ADVANC*Z3 END IF END IF C C Material parameters (defining ISB1I, ICB1I, YL(1) to YL(6)): IF(ICB1I.NE.0) THEN IF(YI(3).NE.YI3.OR.YI(4).NE.YI4.OR.YI(5).NE.YI5) THEN ICB1I=0 END IF END IF IF(ICB1I.EQ.0) THEN CALL BLOCK(.TRUE.,YI(3),0,ISB1I,I,ISB2,ICB2,FAUX) ISB1I=ISB2 C Defining locations of the array FSRFCA of the common /INITC/: DO 11 I=1,NSRFCA CALL SRFC2(NSRFC()+I,YI(3),FAUX) FSRFCA(I)=FAUX(1) 11 CONTINUE ELSE ICB2=IABS(ICB1I) ENDIF IF(ICB2.EQ.0) THEN IEND=72 RETURN ENDIF IY(2)=0 IY(6)=0 IY(8)=ICB2 CALL CODE(IY,KODIND,ICBNEW,IEND) IF(IEND.EQ.22) THEN IEND=72 RETURN ELSE IF(IEND.NE.0) THEN WRITE(*,'(A,I2)') ' Subroutine CODE reports IEND=',IEND PAUSE 'Error 611 in INIT2: Wrong function of subroutine CODE' STOP END IF IF(ICB2.NE.IABS(ICBNEW)) THEN PAUSE 'Error 612 in INIT2: Wrong function of subroutine CODE' STOP END IF IF(ICB1I.EQ.0.OR.ICB1I.NE.ICBNEW) THEN ICB1I=ICBNEW CALL PARM2(IABS(ICBNEW),YI(3),UP,US,YLI(3),QP,QS) CALL VELOC(ICBNEW,UP,US,QP,QS,YLI(1),YLI(2),VD,QL) V0=VD(1) YLI(4)=VD(2) YLI(5)=VD(3) YLI(6)=VD(4) ENDIF C C Important quantities at the initial point of the ray (C.R.T.6.1): C Slowness derivatives with respect to ray parameters AUX1=-(YLI(4)*F(2,1)+YLI(5)*F(2,2)+YLI(6)*F(2,3))/(V0*V0) AUX2=-(YLI(4)*F(3,1)+YLI(5)*F(3,2)+YLI(6)*F(3,3))/(V0*V0) U1=E11*AUX1+E12*AUX2 U2=E12*AUX1+E22*AUX2 C Slowness vector AUX1=( C22*F(2,4)-C12*F(3,4))/DETC AUX2=(-C12*F(2,4)+C11*F(3,4))/DETC AUX3=V0**(-2)-AUX1*F(2,4)-AUX2*F(3,4) IF(AUX3.LE.0.) THEN C Evanescent wave IEND=74 RETURN END IF AUX3=SQRT(AUX3) YI(6)=F21*AUX1+F31*AUX2+Z1*AUX3 YI(7)=F22*AUX1+F33*AUX2+Z2*AUX3 YI(8)=F23*AUX1+F33*AUX2+Z3*AUX3 C Space-time scalar products of vectors tangent to the surface T1=F(2,4)*V0 T2=F(3,4)*V0 C11=C11-T1*T1 C12=C12-T1*T2 C22=C22-T2*T2 DETC=SQRT(C11*C22-C12*C12) IF(INIDIM.NE.1) THEN C Initial surface or initial point: C Square root of the matrix C AUX1=SQRT(C11+C22+DETC+DETC) C11=(C11+DETC)/AUX1 C12=C12/AUX1 C22=(C22+DETC)/AUX1 C Inverse square root of the matrix C B11= C22/DETC B12=-C12/DETC B22= C11/DETC C First basis vector of ray-centred coordinate system AUX3=V0*(B11*T1+B12*T2) YI(9) =F21*B11+F31*B12-YI(6)*AUX3 YI(10)=F22*B11+F32*B12-YI(7)*AUX3 YI(11)=F23*B11+F33*B12-YI(8)*AUX3 C Geometrical spreading matrix Q YI(12)=C11*E11+C12*E12 YI(13)=C12*E11+C22*E12 YI(16)=C11*E12+C12*E22 YI(17)=C12*E12+C22*E22 C Matrix P FF11=F(4,4)-YI(6)*F(4,1)-YI(7)*F(4,2)-YI(8)*F(4,3)-T1*U1 FF12=F(5,4)-YI(6)*F(5,1)-YI(7)*F(5,2)-YI(8)*F(5,3) FF21=FF12-T2*U1 FF12=FF12-T1*U2 FF22=F(6,4)-YI(6)*F(6,1)-YI(7)*F(6,2)-YI(8)*F(6,3)-T2*U2 YI(14)=B11*FF11+B12*FF21 YI(15)=B12*FF11+B22*FF21 YI(18)=B11*FF12+B12*FF22 YI(19)=B12*FF12+B22*FF22 ELSE C Initial line: C Infinite part of the inverse square root of the matrix C B11=(1.-E11)/DETC B12= -E12 /DETC B22=(1.-E22)/DETC C Matrix CB=1-C*B CB11=1.-C11*B11+C12*B12 CB21= -C12*B11+C22*B12 CB12= -C11*B12+C12*B22 CB22=1.-C12*B12+C22*B22 C E-projection of the finite part of the inverse square root of C TRECE=SQRT(E11*C11+2.*E12*C12+E22*C22) ECB11=(E11*CB11+E12*CB21)/TRECE ECB21=(E12*CB11+E22*CB21)/TRECE ECB12=(E11*CB12+E12*CB22)/TRECE ECB22=(E12*CB12+E22*CB22)/TRECE C First basis vector of ray-centred coordinate system AUX1=B11+ECB11 AUX2=B12+ECB12 BT1=AUX1*T1+AUX2*T2 BT2=(B12+ECB21)*T1+(B22+ECB22)*T2 AUX3=V0*BT1 YI(9) =F21*AUX1+F31*AUX2-YI(6)*AUX3 YI(10)=F22*AUX1+F32*AUX2-YI(7)*AUX3 YI(11)=F23*AUX1+F33*AUX2-YI(8)*AUX3 C Geometrical spreading matrix Q YI(12)=E11*TRECE YI(13)=E12*TRECE YI(16)=E12*TRECE YI(17)=E22*TRECE C Matrix P FF11=F(4,4)-YI(6)*F(4,1)-YI(7)*F(4,2)-YI(8)*F(4,3) FF12=F(5,4)-YI(6)*F(5,1)-YI(7)*F(5,2)-YI(8)*F(5,3) FF22=F(6,4)-YI(6)*F(6,1)-YI(7)*F(6,2)-YI(8)*F(6,3) FFE11=FF11*E11+FF12*E12 FFE21=FF12*E11+FF22*E12 FFE12=FF11*E12+FF12*E22 FFE22=FF12*E12+FF22*E22 YI(14)=B11*FF11+B12*FF12+ECB11*FFE11+ECB12*FFE21-BT1*U1 YI(15)=B12*FF11+B22*FF12+ECB12*FFE11+ECB22*FFE21-BT2*U1 YI(18)=B11*FF12+B12*FF22+ECB11*FFE12+ECB12*FFE22-BT1*U2 YI(19)=B12*FF12+B22*FF22+ECB12*FFE12+ECB22*FFE22-BT2*U2 END IF C Take-off parameters YI(20)=PAR1 YI(21)=PAR2 C C Modification of the initial travel time: IF(ADVANC.NE.0.) THEN YI(1)=YI(1)+(YI(3)-F(1,1))*YI(6) * +(YI(4)-F(1,2))*YI(7) * +(YI(5)-F(1,3))*YI(8) END IF C C C Initial values for the complete ray tracing (C.R.T.6.2): DO 21 I=1,6 YL(I)=YLI(I) 21 CONTINUE DO 22 I=1,11 Y(I)=YI(I) 22 CONTINUE IF(ICB1I.GE.0) THEN NY=27+2 ELSE NY=27+8 ENDIF DO 23 I=12,NY Y(I)=0.0 23 CONTINUE Y(12)=1.0 Y(17)=1.0 Y(22)=1.0 Y(27)=1.0 Y(28)=1.0 IF(NY.GE.34) Y(34)=1.0 YY(1)=0.0 YY(2)=UEB YY(3)=0.0 YY(4)=0.0 YY(5)=0.0 IY(1)=NY IY(2)=KODIND IY(3)=0 IY(4)=ISB1I IY(5)=ICB1I IY(6)=0 C Note: IY(7),IY(8) may be undefined IY(7)=0 IY(8)=0 IY(9)=0 IY(10)=0 IY(11)=0 IY(12)=0 RETURN END C C======================================================================= C SUBROUTINE INIS1(LUN,NFUNC) INTEGER LUN,NFUNC C C Subroutine INIS1 reads the input data for the initial points of rays C and stores them in common block /INISC/, and if required, calls the C subroutine VAL1 to read the input data for the interpolated functions C of two variables (ray parameters), to determine the coefficients C necessary to compute an interpolatory function on a two dimensional C rectangular grid, and to store them in the memory. The functions C determined can be represented as a tensor product of splines under C tension. For actual mapping of points it is necessary to call the C subroutine INIS2, which also returns the first and second partial C derivatives. C C Input: C LUN... Logical unit number of the external input device C containing the input data. C NFUNC...Number of the functions required to be defined during the C current invocation of INIS1. Since the functions C specified in the input data do not coincide with the C required functions but are transformed to them, NFUNC need C not equal the number of functions specified in the input C data. C None of the input parameters are altered. C C No output. C C Common block /INISC/: INCLUDE 'initd.inc' C All the storage locations of the common block are defined in this C subroutine. C C Subroutines and external functions required: EXTERNAL VAL1 C VAL1, SORTV, READV... File 'val.for' of the package 'MODEL'. C CURVN1 or CURVB1 (alternatives), SURFB1, VAL3B1, SNHCSH, VGEN, C TERMS, TRIDEC, TRISOL... Subroutine package 'FITPACK' C (file 'fit.for' of the package 'MODEL'). C C Error messages: C 601... Small number of input functions: C The number of input functions is less than the number of C functions necessary to describe coordinates and travel C time along the initial surface. C 602... Wrong value of INIPAR: C For INIDIM=1, there must be INIPAR=1 or INIPAR=2. C C Date: 1994, February 26 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER LUSRC PARAMETER (LUSRC=9) INTEGER LFUNC,MFUNC CHARACTER*80 TEXTI CHARACTER*3 TFUNCT DATA TFUNCT/' '/ C C LUSRC...Logical unit number of the source-point file. The file is C opened and closed here. C LFUNC...Difference between the number NFUNC of the required C functions and the number of input functions specifying C them. C MFUNC...Number of functions specified in the input data. C TEXTI...String of 80 characters for various purposes. C C (1) Reading the name of the input data READ(LUN,*) TEXTI C C (2) Quantities defining the kind of initial conditions INIDIM=-1 INIPAR= 2 ADVANC= 0. READ(LUN,*) INIDIM,INIPAR,ADVANC C C (3) Data describing the initial point, curve or surface. IF(INIDIM.LT.0) THEN TEXTI='src.dat' READ(LUN,*) TEXTI OPEN(LUSRC,FILE=TEXTI,STATUS='OLD') READ(LUSRC,*) (TEXTI,I=1,20) TTINI=0. READ(LUSRC,*) TEXTI,X1INI,X2INI,X3INI,TTINI CLOSE(LUSRC) ELSE IF(INIDIM.EQ.0) THEN TTINI=0. READ(LUN,*) X1INI,X2INI,X3INI,TTINI ELSE IF(INIDIM.EQ.2.AND.INIPAR.GT.0) THEN READ(LUN,*) X1INI,X2INI,X3INI LFUNC=2 ELSE LFUNC=0 IF(INIDIM.EQ.1.AND.(INIPAR.LT.1.OR.INIPAR.GT.2)) THEN PAUSE 'Error 602 in INIS1: Wrong value of INIPAR' STOP END IF END IF READ(LUN,*) MFUNC IF(NFUNC-LFUNC.LE.MFUNC) THEN CALL VAL1(LUN,3,MFUNC,1,TFUNCT) ELSE PAUSE 'Error 601 in INIS1: Small number of input functions' STOP END IF END IF C RETURN END C C======================================================================= C SUBROUTINE INIS2(NFUNC,PAR1,PAR2,E,F,IEND) INTEGER NFUNC,IEND REAL PAR1,PAR2,E(3),F(6,NFUNC) C C Subroutine INIS2 evaluates the functional values and the derivatives C of the functions describing the initial surface. The first three C functions of given ray parameters are coordinates of the point C corresponding to the given ray parameters, the fourth function is the C initial value of the travel time. The single initial point common to C all rays or the initial line are treated as singular limiting cases of C the initial surface. The input data specifying the functions must C have been read by the subroutine INIS1. C C Input: C NFUNC...Number of functions required. It is assumed to be 4 in C this version (three coordinates and the travel time). C PAR1,PAR2... Ray parameters. C E... Array of the dimension at least 3. C F... Array of the dimension at least 6*NFUNC. C None of the input parameters except E and F are altered. C C Output: C E... Array containing the components E11, E12, E22 of the 2*2 C symmetric projection matrix onto the tangent space to the C ray parameter's manifold. For a non-degenerate initial C surface, E is the identity matrix. For a single initial C point, E is the zero matrix. For the initial line, E is C the projection matrix of the rank 1. Note that a C projection matrix E satisfies the relation E*E=E. C F(1:6,I)... For a non-degenerate initial surface, the value and C the first and second partial derivatives F, F1, F2, F11, C F12, F22 of the I-th function F(PAR1,PAR2). Note that C F1 = E11,E12 * F1 C F2 E12,E22 F2 , C and C F11,F12 = E11,E12 * F11,F12 * E11,E12 C F12,F22 E12,E22 F12,F22 E12,E22 . C Thus, in a degenerate case (i.e. if E is not the identity C matrix), the first derivatives are modified in the C following way, C F1 = F1 + F31 - E11,E12 * F31 C F2 F2 F32 E12,E22 F32 , C and second derivatives are modified as follows: C F11,F12 = F11,F12 + F311,F312 - E11,E12*F311,F312*E11,E12 C F12,F22 F12,F22 F312,F322 E12,E22 F312,F322 E12,E22. C Here F31, F32, F311, F312 and F322 are the derivatives of C F1, F2, F11, F12 and F22 with respect to the small C parameter (e.g. a radius) which shrinks to zero upon an C initial line or at a single initial point. C IEND... Information on the initial point of the ray: C 0... Computation of the ray may follow. C 71... There is no ray corresponding to the given ray C parameters. E.g., the given parameters do not belong to C the domain of the initial surface. C C Common block /INISC/: INCLUDE 'initd.inc' C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL METRIC,VAL2,SPHERE,SQRT3 C METRIC..File 'metric.for' of the package 'MODEL'. C VAL2... File 'val.for' of the package 'MODEL'. C CURV2D OR CURVBD (alternatives), SURFBD, VAL3BD, SNHCSH, DSPLNZ, C INTRVL... Subroutine package 'FITPACK' (file 'fit.for' of C the package 'MODEL'). C SPHERE,SQRT3... This file. C C Date: 1996, May 9 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I,I1,I2,I11,I22,LFUNC REAL AUX1,AUX2,DUMMY,R(3),G(12),FF(6,3) C C I... Auxiliary loop variable. C I1,I11..Array subscripts of the first and second derivatives with C respect to INIPAR-th ray parameter. C I2,I22..Array subscripts of the first and second derivatives with C respect to the other than INIPAR-th ray parameter. C LFUNC...Difference between the numbers of required (NFUNC) and C interpolated functions C AUX1,AUX2... Auxiliary storage locations. C DUMMY...Dummy storage location. C R... Array used for general coordinates or ray parameters. C G... G(1)-G(6)... Covariant components of the symmetric 3*3 C metric tensor, or contravariant components of the C symmetric 3*3 matrix of the basis vectors of the local C Cartesian coordinate system (i.e. the square root of the C conravariant metric tensor). C G(7)-G(12)... Contravariant components of the symmetric C 3*3 metric tensor. C FF... General-purpose auxiliary array. Used to store local C Cartesian coordinates and their derivatives with respect C to the ray parameters. Temporarily used also as dummy C storage location for Christoffel symbols, for the last C interpolated function, e.t.c. C C....................................................................... C IEND=0 IF(INIDIM.LE.1.OR.INIPAR.GE.1) THEN C Matrix E of local Cartesian coordinate system basis vectors R(1)=X1INI R(2)=X2INI R(3)=X3INI CALL METRIC(R,DUMMY,G,FF) END IF IF(INIDIM.LE.0) THEN C Single initial point: C Projection matrix E(1)=0. E(2)=0. E(3)=0. C Contravariant components of the symmetric 3*3 matrix of the C basis vectors of the local Cartesian coordinate system CALL SQRT3(G(7),G) C Mapping of the ray parameters onto a unit sphere CALL SPHERE(INIPAR,PAR1,PAR2,FF,IEND) IF(IEND.NE.0) THEN RETURN END IF C Required functions (3 general coordinates and a travel time) F(1,1)=X1INI F(1,2)=X2INI F(1,3)=X3INI F(1,4)=TTINI DO 14 I=2,6 F(I,1)=G(1)*FF(I,1)+G(2)*FF(I,2)+G(4)*FF(I,3) F(I,2)=G(2)*FF(I,1)+G(3)*FF(I,2)+G(5)*FF(I,3) F(I,3)=G(4)*FF(I,1)+G(5)*FF(I,2)+G(6)*FF(I,3) F(I,4)=0. 14 CONTINUE ELSE C Initial line or initial surface: C Interpolated functions R(1)=PAR1 R(2)=PAR2 R(3)=0. IF(INIDIM.EQ.1) THEN R(3-INIPAR)=0. END IF IF(INIDIM.EQ.2.AND.INIPAR.GT.0) THEN LFUNC=2 ELSE LFUNC=0 END IF DO 21 I=LFUNC+1,NFUNC-1 CALL VAL2(3,I-LFUNC,1,R,F(1,I),DUMMY) F(4,I)=F(5,I) F(5,I)=F(6,I) F(6,I)=F(1,I+1) 21 CONTINUE CALL VAL2(3,NFUNC-LFUNC,1,R,FF,DUMMY) F(1,NFUNC)=FF(1,1) F(2,NFUNC)=FF(2,1) F(3,NFUNC)=FF(3,1) F(4,NFUNC)=FF(5,1) F(5,NFUNC)=FF(6,1) F(6,NFUNC)=FF(1,2) IF(INIDIM.EQ.1) THEN C Initial line: C Covariant components of the vector tangent to the initial line I1=1+INIPAR FF(6,1)=G(1)*F(I1,1)+G(2)*F(I1,2)+G(4)*F(I1,3) FF(6,2)=G(2)*F(I1,1)+G(3)*F(I1,2)+G(5)*F(I1,3) FF(6,3)=G(4)*F(I1,1)+G(5)*F(I1,2)+G(6)*F(I1,3) C Contravariant unit vector tangent to the initial line AUX2=F(I1,1)*FF(6,1)+F(I1,2)*FF(6,2)+F(I1,3)*FF(6,3) AUX1=SQRT(AUX2) FF(1,1)=F(I1,1)/AUX1 FF(1,2)=F(I1,2)/AUX1 FF(1,3)=F(I1,3)/AUX1 C Derivative of the unit vector tangent to the initial line I11=2*I1 AUX2=(F(I11,1)*FF(6,1)+F(I11,2)*FF(6,2)+F(I11,3)*FF(6,3))/AUX2 FF(2,1)=(F(I11,1)-FF(1,1)*AUX2)/AUX1 FF(2,2)=(F(I11,2)-FF(1,2)*AUX2)/AUX1 FF(2,3)=(F(I11,3)-FF(1,3)*AUX2)/AUX1 C Covariant vector normal to the interpolated surface I2=5-I1 FF(3,1)=FF(1,2)*F(I2,3)-FF(1,3)*F(I2,2) FF(3,2)=FF(1,3)*F(I2,1)-FF(1,1)*F(I2,3) FF(3,3)=FF(1,1)*F(I2,2)-FF(1,2)*F(I2,1) C Derivative of the vector normal to the interpolated surface FF(4,1)=FF(2,2)*F(I2,3)-FF(2,3)*F(I2,2)+ * FF(1,2)*F(5,3) -FF(1,3)*F(5,2) FF(4,2)=FF(2,3)*F(I2,1)-FF(2,1)*F(I2,3)+ * FF(1,3)*F(5,1) -FF(1,1)*F(5,3) FF(4,3)=FF(2,1)*F(I2,2)-FF(2,2)*F(I2,1)+ * FF(1,1)*F(5,2) -FF(1,2)*F(5,1) C Contravariant components FF(5,1)=G(7) *FF(3,1)+G(8) *FF(3,2)+G(10)*FF(3,3) FF(5,2)=G(8) *FF(3,1)+G(9) *FF(3,2)+G(11)*FF(3,3) FF(5,3)=G(10)*FF(3,1)+G(11)*FF(3,2)+G(12)*FF(3,3) FF(6,1)=G(7) *FF(4,1)+G(8) *FF(4,2)+G(10)*FF(4,3) FF(6,2)=G(8) *FF(4,1)+G(9) *FF(4,2)+G(11)*FF(4,3) FF(6,3)=G(10)*FF(4,1)+G(11)*FF(4,2)+G(12)*FF(4,3) C Required functions (3 general coordinates and a travel time) E(2)=0. IF(INIPAR.LE.1) THEN E(1)=1. E(3)=0. AUX1=COS(PAR2) AUX2=SIN(PAR2) ELSE E(1)=0. E(3)=1. AUX1=COS(PAR1) AUX2=SIN(PAR1) END IF I22=10-I11 DO 24 I=1,3 F(I22,I)=-AUX2*F(I2,I)+AUX1*FF(5,I) F(I2,I) = AUX1*F(I2,I)+AUX2*FF(5,I) F(5,I) = AUX1*F(5,I) +AUX2*FF(6,I) 24 CONTINUE F(I22,4)=0. F(I2,4)=0. F(5,4)=0. ELSE C Initial surface: C Projection matrix E(1)=1. E(2)=0. E(3)=1. C Required functions (3 general coordinates and a travel time) IF(INIPAR.GE.1) THEN C Contravariant components of the symmetric 3*3 matrix of the C basis vectors of the local Cartesian coordinate system CALL SQRT3(G(7),G) C Mapping of the ray parameters onto a unit sphere CALL SPHERE(INIPAR,PAR1,PAR2,FF,IEND) IF(IEND.NE.0) THEN RETURN END IF C Local Cartesian coordinates DO 33 I=1,3 FF(6,I)=F(6,3)*FF(1,I)+2.*F(3,3)*FF(3,I)+F(1,3)*FF(6,I) FF(5,I)=F(5,3)*FF(1,I)+F(3,3)*FF(2,I)+F(3,3)*FF(1,I) * +F(1,3)*FF(5,I) FF(4,I)=F(4,3)*FF(1,I)+2.*F(2,3)*FF(2,I)+F(1,3)*FF(4,I) FF(3,I)=F(3,3)*FF(1,I)+F(1,1)*FF(3,1) FF(2,I)=F(2,3)*FF(1,I)+F(1,1)*FF(2,1) FF(1,I)=F(1,3)*FF(1,I) 33 CONTINUE C General coordinates DO 34 I=1,6 F(I,1)=G(1)*FF(I,1)+G(2)*FF(I,2)+G(4)*FF(I,3) F(I,2)=G(2)*FF(I,1)+G(3)*FF(I,2)+G(5)*FF(I,3) F(I,3)=G(4)*FF(I,1)+G(5)*FF(I,2)+G(6)*FF(I,3) 34 CONTINUE F(1,1)=F(1,1)+X1INI F(1,2)=F(1,2)+X2INI F(1,3)=F(1,3)+X3INI END IF END IF END IF RETURN END C C======================================================================= C SUBROUTINE SQRT3(B,A) REAL B(6),A(6) C C Subroutine SQRT3 evaluates the square root A of the given real-valued C positive-semidefinite symmetric 3*3 matrix B. The square root is the C real-valued positive-semidefinite symmetric 3*3 matrix A satisfying C the equation A*A=B. C C Input: C B... Array of dimension at least 6, containing the components C B11, B12, B22, B13, B23, B33 of the given symmetric 3*3 C matrix B. C A... Array of dimension at least 6. C The input parameter B is not altered. C C Output. C A... Array containing the components A11, A12, A22, A13, A23, C A33 of the symmetric 3*3 matrix a which is the square root C of the given matrix B. C C No subroutines and external functions required. C C Error messages: C 614... Matrix is not positive-semidefinite: C Input matrix B is not positive-semidefinite. C 615... This program branch is not coded: C The square root of general symmetric matrix B has not been C coded. At present, only the square root of diagonal C matrix can be evaluated. C C Date: 1990, January 22 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C No auxiliary storage locations C IF(B(2).EQ.0..AND.B(4).EQ.0..AND.B(5).EQ.0.) THEN C Diagonal matrix IF(B(1).LT.0..OR.B(3).LT.0..OR.B(6).LT.0.) THEN C B is not positive-semidefinite matrix PAUSE'Error 614 in SQRT3: Matrix is not positive-semidefinite' STOP ELSE A(1)=SQRT(B(1)) A(2)=0. A(3)=SQRT(B(3)) A(4)=0. A(5)=0. A(6)=SQRT(B(6)) END IF ELSE C General symmetric matrix PAUSE 'Error 615 in SQRT3: This program branch is not coded' STOP END IF RETURN END C C======================================================================= C SUBROUTINE SPHERE(INIPAR,PAR1,PAR2,FF,IEND) INTEGER INIPAR,IEND REAL PAR1,PAR2,FF(6,3) C C Subroutine SPHERE transforms spherical coordinates PAR1, PAR2 into the C Cartesian coordinates of the corresponding point on the unit sphere. C It also evaluates the first and second derivatives of the Cartesian C coordinates with respect to PAR1 and PAR2. C C Input: C INIPAR..Determines the type of spherical coordinates: C INIPAR.LT.0: The same as for IABS(INIPAR), but with the C unit vector (T1,T2,T3) tangent to the ray changed to C (T1,-T2,-T3). C INIPAR.EQ.1: Polar-like spherical coordinates (colatitude, C longitude). C INIPAR.GE.2: Geographic-like spherical coordinates C (longitude, latitude). C INIPAR.EQ.3: The unit vector tangent to the ray, C expressed in the local Cartesian coordinate system C which basis vectors are given by the square root of C the metric tensor at the initial point, is given by C T1=PAR1*SIN(R)/R C T2=PAR2*SIN(R)/R C T3= COS(R) C with R=SQRT(PAR1*PAR1+PAR2*PAR2). C If R.GT.2*PI, the ray is reported out of the C ray-parameter domain: IEND=71. C PAR1,PAR2... Ray parameters. C FF... Array of the dimension at least 6*3. C None of the input parameters except FF are altered. C C Output: C FF(1:6,I)...I-th Cartesian coordinate of the point on the unit C sphere given by PAR1, PAR2, and its first and second C partial derivatives with respect to PAR1 and PAR2 in the C order FF, FF1, FF2, FF11, FF12, FF22. C IEND... Information on the initial point of the ray: C 0... Computation of the ray may follow. C 71... There is no ray corresponding to the given ray C parameters. I.e., the given parameters are outside the C domain allowed. C C No subroutines and external functions required. C C Date: 1995, July 31 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: REAL C1,C2,S1,S2,R,R1,R2,R11,R12,R22,R111,R112,R122,R222 C C C1,C2...Cosines of the take-off angles at a single initial point. C S1,S2...Sines of the take-off angles at a single initial point. C R,R1,R2,R11,R12,R22,R111,R112,R122,R222... C R=SQRT(PAR1*PAR1+PAR2*PAR2) and its partial derivatives. C IEND=0 IF(IABS(INIPAR).LE.2) THEN C1=COS(PAR1) S1=SIN(PAR1) C2=COS(PAR2) S2=SIN(PAR2) IF(IABS(INIPAR).LE.1) THEN C Polar-like spherical coordinates IF(S1.LT.0.) THEN IEND=71 RETURN END IF FF(1,1)= S1*C2 FF(1,2)= S1*S2 FF(1,3)= C1 IF(S1.EQ.0.) THEN C Avoiding null vectors S1=0.000001 END IF FF(2,1)= C1*C2 FF(2,2)= C1*S2 FF(2,3)=-S1 FF(3,1)=-S1*S2 FF(3,2)= S1*C2 FF(3,3)= 0. FF(4,1)=-S1*C2 FF(4,2)=-S1*S2 FF(4,3)=-C1 FF(5,1)=-C1*S2 FF(5,2)= C1*C2 FF(5,3)= 0. FF(6,1)=-S1*C2 FF(6,2)=-S1*S2 FF(6,3)= 0. ELSE C Geographic-like spherical coordinates IF(C2.LT.0.) THEN IEND=71 RETURN END IF FF(1,1)= C2*C1 FF(1,2)= C2*S1 FF(1,3)= S2 IF(C2.EQ.0.) THEN C Avoiding null vectors C2=0.000001 END IF FF(2,1)=-C2*S1 FF(2,2)= C2*C1 FF(2,3)= 0. FF(3,1)=-S2*C1 FF(3,2)=-S2*S1 FF(3,3)= C2 FF(4,1)=-C2*C1 FF(4,2)=-C2*S1 FF(4,3)= 0. FF(5,1)= S2*S1 FF(5,2)=-S2*C1 FF(5,3)= 0. FF(6,1)=-C2*C1 FF(6,2)=-C2*S1 FF(6,3)=-S2 END IF ELSE R=SQRT(PAR1*PAR1+PAR2*PAR2) R1=PAR1/R R2=PAR2/R R11= R2*R2/R R12=-R1*R2/R R22= R1*R1/R R111= -3.*R11*R1 /R R112=(-R2/R-3.*R12*R1)/R R122=(-R1/R-3.*R12*R2)/R R222= -3.*R22*R2 /R C1=COS(R) S1=SIN(R) FF(1,1)= S1*R1 FF(1,2)= S1*R2 FF(1,3)= C1 FF(2,1)= C1*R1*R1+S1*R11 FF(2,2)= C1*R1*R2+S1*R12 FF(2,3)=-FF(1,1) FF(3,1)= FF(2,2) FF(3,2)= C1*R2*R2+S1*R22 FF(3,3)=-FF(1,2) FF(4,1)=-S1*R1*R1*R1 +3.*C1*R1*R11+S1*R111 FF(4,2)=-S1*R1*R1*R2+C1*R11*R2+2.*C1*R1*R12+S1*R112 FF(4,3)=-FF(2,1) FF(5,1)= FF(4,2) FF(5,2)=-S1*R1*R2*R2+C1*R1*R22+2.*C1*R2*R12+S1*R122 FF(5,3)=-FF(2,2) FF(6,1)= FF(5,2) FF(6,2)=-S1*R2*R2*R2 +3.*C1*R2*R22+S1*R222 FF(6,3)=-FF(3,2) END IF C IF(INIPAR.LT.0) THEN DO 12 J=2,3 DO 11 I=1,6 FF(I,J)=-FF(I,J) 11 CONTINUE 12 CONTINUE END IF C RETURN END C C======================================================================= C