C
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 Date: 2000, May 19 C Coded by Ludek Klimes C C....................................................................... C C This file consists of the following subroutines: 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 INIT1 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 INIT2 C INIT5...Subroutine determining whether ray tracing should be C repeated for another source point. 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 INIS1 C INIS2...Sample subroutine returning the functional values and C their first and second derivatives, of the functions C describing the initial surface. C INIS2 C SQRT3...Subroutine evaluating the square root of the given C real-valued positive-semidefinite symmetric 3*3 matrix. C SQRT3 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 SPHERE 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....................................................................... C C Description of data files: C C Input parameters taken from the input SEP parameter file for the C Complete Ray Tracing program C INIDIM C INIPAR C ADVANC C are described in file crtin.for. C C C Input file SRC containing the coordinates of the initial point: C The data are read only if INIDIM.LE.0. C If there are several source points, ray tracing is repeated for C all source points, repeating all elementary waves specified in C data set CODE for each source point. C Data in data sets RPAR and C WRIT are not repeated, they should be C specified for all elementary waves of all source points. C (1) Several strings terminated by / (a slash). C The simplest way is to submit just the /. C (2) For each source point data (2.1): C (2.1) '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: X2INI=0, X3INI=0, TTINI=0. C (3) / (a slash) or end of file. C Example of data set SRC C C C Input data INIT for the initial points of rays: C The data are read only if INIDIM.GT.0. 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) Data describing the initial point, curve or surface. C For INIDIM=1: C (2B) 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 (2B) 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 (2A), (2B): C (2A) X1INI,X2INI,X3INI... Coordinates of a given point, C see the description of the input data (1). C (2B) 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 The structure of the input data (2B) is given by the subroutine VAL1 C and is described below. C Examples of data set INIT: C C Plane wave incident at the model bottom C C Zero-offset rays (exploding reflector) C C Above mentioned input data (2B) 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 Input data for one function 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 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....................................................................... C C Storage in the memory: C Input data INIT (2A) 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 include file 'initd.inc'. C initd.inc C C======================================================================= 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 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 INIS1 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 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 IEND 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 dcrt.inc C None of the storage locations of the common block are altered. C C Common block /INISC/: INCLUDE 'initd.inc' C initd.inc C None of the storage locations of the common block are altered. C C Common block /INITC/: INCLUDE 'initc.inc' C 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 Date: 2000, May 19 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 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 IY(11)=0 IF(MODCRT.GE.3) THEN C 613 CALL ERROR('613 in INIT2: 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. 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 C 617 CALL ERROR('617 in INIT2: Strange common initial point') C Initial conditions determined by subroutine INIS2 do not C resemble common initial point. Contact the authors. END IF ELSE IF(E11*E22-E12*E12.LE.0.001) THEN C Initial line: IF(INIDIM.NE.1) THEN C 618 CALL ERROR('618 in INIT2: Strange initial line') C Initial conditions determined by subroutine INIS2 do not C resemble initial line. Contact the authors. END IF ELSE C Initial surface: IF(INIDIM.NE.2) THEN C 619 CALL ERROR('619 in INIT2: Strange initial surface') C Initial conditions determined by subroutine INIS2 do not C resemble initial surface. Contact the authors. 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(YI(3),0,0,I,ISB1I,ICB2) ELSE ICB2=IABS(ICB1I) ENDIF 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 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 C 611 WRITE(*,'(A,I2)') ' Subroutine CODE reports IEND=',IEND CALL ERROR('611 in INIT2: 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. END IF IF(ICB2.NE.IABS(ICBNEW)) THEN C 612 CALL ERROR('612 in INIT2: 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. 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 C IY(11) is initiallized in the beginning of this subroutine. IY(12)=0 RETURN END C C======================================================================= C C C SUBROUTINE INIT5(IRAY) INTEGER IRAY C C Subroutine INIT5 determines whether ray tracing should be repeated for C another source point (common initial point of rays). C C No input. C C Output: C IRAY... IRAY=0: There is another source point. The source point C has replaced the previous source point in locations C X1INI, X2INI, X3INI, TTINI of common block /INISC/ C IRAY=-1: There is not another source point. C C Common block /INISC/: INCLUDE 'initd.inc' C initd.inc C C No subroutines and external functions required. C C Date: 1999, September 23 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I C IF(NSRC.GT.0) THEN C There is another source point X1INI=X1SRC(1) X2INI=X2SRC(1) X3INI=X3SRC(1) TTINI=TTSRC(1) NSRC=NSRC-1 DO 10 I=1,NSRC X1SRC(I)=X1SRC(I+1) X2SRC(I)=X2SRC(I+1) X3SRC(I)=X3SRC(I+1) TTSRC(I)=TTSRC(I+1) 10 CONTINUE IRAY=0 ELSE IRAY=-1 END IF RETURN END C C======================================================================= 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 May be zero for a point source. 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 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 RSEP3I,RSEP3R,RSEP3T,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 Date: 1999, September 23 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER LUSRC PARAMETER (LUSRC=9) REAL UNDEF PARAMETER (UNDEF=-999999.) INTEGER LFUNC,MFUNC,I 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 Quantities defining the kind of initial conditions CALL RSEP3I('INIDIM',INIDIM,0) CALL RSEP3I('INIPAR',INIPAR,2) CALL RSEP3R('ADVANC',ADVANC,0.) C NSRC=0 IF(INIDIM.LE.0) THEN C Point source CALL RSEP3T('SRC',TEXTI,'src.dat') OPEN(LUSRC,FILE=TEXTI,STATUS='OLD') READ(LUSRC,*) (TEXTI,I=1,20) X2INI=0. X3INI=0. TTINI=0. READ(LUSRC,*) TEXTI,X1INI,X2INI,X3INI,TTINI DO 18 I=1,MSRC X1SRC(I)=UNDEF X2SRC(I)=0. X3SRC(I)=0. TTSRC(I)=0. READ(LUSRC,*,END=19) TEXTI,X1SRC(I),X2SRC(I),X3SRC(I),TTSRC(I) IF(X1SRC(I).EQ.UNDEF) THEN GO TO 19 END IF 18 CONTINUE C 604 CALL ERROR('604 in INIS1: Too many source points') C Dimension MSRC of arrays X1SRC, X2SRC, X3SRC, TTSRC in file C initd.inc is smaller than the number C of source points. 19 CONTINUE NSRC=I-1 CLOSE(LUSRC) ELSE IF(LUN.LE.0) THEN C 603 CALL ERROR('603 in INIS1: No input device') C For INIDIM.GT.0, there must be specified input file INIT, C see the main input data for CRT. END IF C C (1) Reading the name of the input data READ(LUN,*) TEXTI C C (3) Data describing the initial point, curve or surface. 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 C 602 CALL ERROR('602 in INIS1: Wrong value of INIPAR') C For INIDIM=1, there must be INIPAR=1 or INIPAR=2. END IF END IF READ(LUN,*) MFUNC IF(NFUNC-LFUNC.LE.MFUNC) THEN CALL VAL1(LUN,3,MFUNC,1,TFUNCT) ELSE C 601 CALL ERROR('601 in INIS1: 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. END IF C END IF RETURN END C C======================================================================= 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 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 contravariant 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 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 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 614 CALL ERROR * ('614 in SQRT3: Matrix is not positive-semidefinite') C Input matrix B is not positive-semidefinite. 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 C 615 CALL ERROR('615 in SQRT3: 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 matrix can C be evaluated. END IF RETURN END C C======================================================================= 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 If SIN(colatitude).LT.0, the ray is reported out of the C ray-parameter domain: IEND=71. C INIPAR.GE.2: Geographic-like spherical coordinates C (longitude, latitude). C If COS(latitude).LT.0, the ray is reported out of the C ray-parameter domain: IEND=71. 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.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: 1999, April 29 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 C Azimuthal equidistant projection R=SQRT(PAR1*PAR1+PAR2*PAR2) IF(R.GT.3.2) THEN IEND=71 RETURN ELSE IF(R.GT.0.) THEN S1=SIN(R) IF(S1.LT.0.) THEN IEND=71 RETURN END IF C1=COS(R) R1=PAR1/R R2=PAR2/R FF(1,1)= S1*R1 FF(1,2)= S1*R2 FF(1,3)= C1 IF(S1.LT.0.004.AND.R.GT.3.1) THEN C Correction for nearly parallel F(2,*) and F(3,*) near C the singularity at R=pi S1=0.004 END IF 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 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) ELSE FF(1,1)= 0. FF(1,2)= 0. FF(1,3)= 1. FF(2,1)= 1. FF(2,2)= 0. FF(2,3)= 0. FF(3,1)= 0. FF(3,2)= 1. FF(3,3)= 0. FF(4,1)= 0. FF(4,2)= 0. FF(4,3)=-1. FF(5,1)= 0. FF(5,2)= 0. FF(5,3)= 0. FF(6,1)= 0. FF(6,2)= 0. FF(6,3)=-1. END IF 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