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 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 FORM 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=0: C INIPAR.LE.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 INIPAR.GE.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 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 (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 INTEGER INIDIM,INIPAR REAL X1INI,X2INI,X3INI,TTINI COMMON/INISC/INIDIM,INIPAR,X1INI,X2INI,X3INI,TTINI SAVE /INISC/ INTEGER MSRFCA PARAMETER (MSRFCA=128) INTEGER ISB1I,ICB1I REAL YLI(6),YI(25),FSRFCA(MSRFCA) COMMON/INITC/ISB1I,ICB1I,YLI,YI,FSRFCA SAVE /INITC/ END C ------------------------------------------------------------------ C INIDIM,INIPAR... INPUT DATA (2). C X1INI,X2INI,X3INI,TTINI... INPUT DATA (3A) IF THEY ARE DEFINED. C ISB1I,ICB1I... INDICES OF A SIMPLE AND A COMPLEX BLOCKS IN WHICH C THE INITIAL POINT OF THE RAY IS SITUATED (SEE C.R.T.6.1). C YLI... ARRAY CONTAINING THE VALUES OF THE QUANTITIES YL(1)-YL(6) C (SEE C.R.T.5.5.4) DESCRIBING THE LOCAL PROPERTIES OF THE C MODEL AT THE INITIAL POINT OF THE RAY, SEE C.R.T.6.1. C THEY MUST NOT BE CHANGED OUTSIDE THE SUBROUTINE INIT2. C YI... ARRAY CONTAINING THE FOLLOWING QUANTITIES DESCRIBING THE C PROPERTIES OF THE RAYS AND OF THE TRAVEL-TIME FIELD, SEE C C.R.T.6.1: C YI(1)...INITIAL TRAVEL TIME. C YI(2)...INITIAL IMAGINARY PART OF THE COMPLEX TRAVEL TIME. C YI(3)-YI(5)... COORDINATES OF THE INITIAL POINT OF THE RAY. C YI(6)-YI(8)... COVARIANT COMPONENTS OF THE INITIAL SLOWNESS C VECTOR. C YI(9)-YI(11)... COVARIANT COMPONENTS OF THE FIRST BASIS VECTOR OF C THE RAY-CENTRED COORDINATE SYSTEM AT THE INITIAL POINT OF C THE RAY (PERPENDICULAR TO THE SLOWNESS VECTOR C YI(6)-YI(8)). C YI(12),YI(16) QR11,QR12 C YI(13),YI(17) QR21,QR22 C YI(14),YI(18) PR11,PR12 C YI(15),YI(19)... PR21,PR22 C ELEMENTS OF THE RAY GEOMETRICAL SPREADING MATRIX QR, AND C OF THE MATRIX PR (SEE C.R.T.,EQ.(5.13)) AT THE INITIAL C POINT OF THE RAY. C YI(20),YI(21)... TAKE-OFF PARAMETERS OF THE RAY. C IN ADDITION TO THE ABOVE QUANTITIES DESCRIBING THE PROPERTIES C DEFINED FOR A SINGLE RAY, THERE ARE ALSO QUANTITIES DESCRIBING C THE PROPERTIES OF THE DISCRETE SYSTEM OF COMPUTED RAYS IN THE C VICINITY OF THE COMPUTED RAY. THESE QUANTITIES Y(22)-Y(25) ARE C NOT DEFINED IN THE SUBROUTINE INIT2: C YI(22)..AREA OF THE ELEMENT OF THE RAY-PARAMETER SURFACE, C CORRESPONDING TO THE RAY, SEE C.R.T.,EQ.(6.1). C YI(23),YI(24),YI(25)... COMPONENTS 11, 12, 22 OF THE SYMMETRIC C MATRIX INVERSE TO THE SPECIFIC MOMENT OF THE ELEMENT OF C THE RAY-PARAMETER SURFACE CORRESPONDING TO THE RAY, SEE C C.R.T.,EQ.(6.2). C COMMON BLOCK /INISC/ IS INCLUDED IN SUBROUTINES INIS1 AND INIS2 C IN ORDER TO COMMUNICATE THE INPUT DATA TO THE SUBROUTINE INIS2. C COMMON BLOCK /INITC/ IS INCLUDED IN EXTERNAL PROCEDURES INIT2, C OUTP (SEE RAYCB), AND MAY BE INCLUDED IN ANY OTHER SUBROUTINE. C THE INDEX OF THE LAST ALLOCATED NUMERIC UNIT OF ARRAY FSRFCA IS C NAMED MSRFCA. IF MSRFCA IS CHANGED, IT MUST BE ADJUSTED IN ALL C SUBROUTINES WHICH INCLUDE COMMON BLOCK /INITC/. 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 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: 1994, JANUARY 15 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/: INTEGER MSRFCA PARAMETER (MSRFCA=128) INTEGER ISB1I,ICB1I REAL YLI(6),YI(25),FSRFCA(MSRFCA) COMMON/INITC/ISB1I,ICB1I,YLI,YI,FSRFCA 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 74... RAY OF THE GENERATED WAVE IS NOT REAL-VALUED. C C COMMON BLOCK /DCRT/ (SEE SUBROUTINE FILE 'RAY.FOR'): INTEGER MEND,MSTOR PARAMETER (MEND=128) PARAMETER (MSTOR=128) INTEGER KSTORE,NEXPS,NHLF,MODCRT REAL STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT,BOUNDR(7) INTEGER NSRFCA,NEND,KEND(MEND),NSTOR,KSTOR(MSTOR) COMMON/DCRT/ KSTORE,NEXPS,NHLF,MODCRT,STORE,STEP,UEB,UEBPP,UEBPH, * UEBHH,UEBDRT,BOUNDR,NSRFCA,NEND,KEND,NSTOR,KSTOR C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C COMMON BLOCK /INITC/: INTEGER MSRFCA PARAMETER (MSRFCA=128) INTEGER ISB1I,ICB1I REAL YLI(6),YI(25),FSRFCA(MSRFCA) COMMON/INITC/ISB1I,ICB1I,YLI,YI,FSRFCA 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 C DATE: 1992, SEPTEMBER 1 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS FOR LOCAL MODEL PARAMETERS: REAL G(12),GAMMA(18),GSQRD,FAUX(10) REAL UP(10),US(10),RO,QP,QS,VP,VS,VD(10),QL C THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A C COMMON BLOCK. THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY C COMMON BLOCK /AUXMOD/ BUT TO SHARE THE MEMORY. COMMON/AUXMOD/ G,GAMMA,GSQRD,FAUX,UP,US,RO,QP,QS,VP,VS,VD,QL C C G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'. C FAUX... AUXILIARY ARRAY TO STORE A FUNCTIONAL VALUE AND ITS C DERIVATIVES. C UP,US,RO,QP,QS... SEE SUBROUTINE PARM2 OF THE FILE 'PARM.FOR'. C UP,US,QP,QS,VP,VS,VD,QL... SEE SUBROUTINE VELOC OF THE FILE C 'MODEL.FOR'. C C....................................................................... C C AUXILIARY STORAGE LOCATIONS: INTEGER KODIND,ICBNEW,ISB2,ICB2,I,NY,INIDIM,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) IF(E11.EQ.0..AND.E12.EQ.0..AND.E22.EQ.0.) THEN C INITIAL POINT: INIDIM=0 ELSE IF(E11*E22-E12*E12.LE.0.001) THEN C INITIAL LINE: INIDIM=1 ELSE C INITIAL SURFACE: INIDIM=2 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 IF(F(1,1).NE.YI(3).OR.F(1,2).NE.YI(4).OR.F(1,3).NE.YI(5)) THEN ICB1I=0 END IF YI(3)=F(1,1) YI(4)=F(1,2) YI(5)=F(1,3) C C MATERIAL PARAMETERS (DEFINING ISB1I, ICB1I, YL(1) TO YL(6)): IF(ICB1I.EQ.0) THEN CALL BLOCK(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 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 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 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 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 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/: INTEGER INIDIM,INIPAR REAL X1INI,X2INI,X3INI,TTINI COMMON/INISC/INIDIM,INIPAR,X1INI,X2INI,X3INI,TTINI 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, JANUARY 15 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS: INTEGER LFUNC,MFUNC CHARACTER*80 TEXTI CHARACTER*3 TFUNCT DATA TFUNCT/' '/ C C TEXTI...STRING OF 80 CHARACTERS FOR VARIOUS PURPOSES. 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 C (1) READING THE NAME OF THE INPUT DATA READ(LUN,*) TEXTI C C (2) QUANTITIES DEFINING THE KIND OF INITIAL CONDITIONS READ(LUN,*) INIDIM,INIPAR C C (3) DATA DESCRIBING THE INITIAL POINT, CURVE OR SURFACE. IF(INIDIM.LT.0) THEN TEXTI='SRC.DAT' READ(LUN,*) TEXTI CLOSE(LUN) OPEN(LUN,FILE=TEXTI,STATUS='OLD') READ(LUN,*) (TEXTI,I=1,20) TTINI=0. READ(LUN,*) TEXTI,X1INI,X2INI,X3INI,TTINI 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) INTEGER NFUNC 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 C COMMON BLOCK /INISC/: INTEGER INIDIM,INIPAR REAL X1INI,X2INI,X3INI,TTINI COMMON/INISC/INIDIM,INIPAR,X1INI,X2INI,X3INI,TTINI 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: 1990, NOVEMBER 18 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 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) 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.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) 21 CONTINUE CALL VAL2(3,NFUNC-LFUNC,1,R,FF,DUMMY) DO 22 I=1,6 F(I,NFUNC)=FF(I,1) 22 CONTINUE 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) 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 A11, A12, A22, A13, A23, A33 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) INTEGER INIPAR 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.LE.1: POLAR-LIKE SPHERICAL COORDINATES (COLATITUDE, C LONGITUDE). C INIPAR.GE.2: GEOGRAPHIC-LIKE SPHERICAL COORDINATES C (LONGITUDE, LATITUDE). 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 C NO SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED. C C DATE: 1989, DECEMBER 19 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS: REAL C1,C2,S1,S2 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 C1=COS(PAR1) S1=SIN(PAR1) C2=COS(PAR2) S2=SIN(PAR2) IF(INIPAR.LE.1) THEN C POLAR-LIKE SPHERICAL COORDINATES FF(1,1)=S1*C2 FF(1,2)=S1*S2 FF(1,3)=C1 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 FF(1,1)=C2*C1 FF(1,2)=C2*S1 FF(1,3)=S2 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 RETURN END C C======================================================================= C