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
