C 
C Program INV1TT to evaluate the derivatives of the travel time with
C respect to the model coefficients.
C
C Version: 5.20
C Date: 1998, August 24
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C Just a preliminary demo version, illustrating the usage of the
C routines designed to calculate the variations of travel times with
C respect to model parameters (FORTRAN77 source code files 'var.for',
C 'spsp.for', 'ap.for').
C
C Program INV1TT assumes all model parameters (coefficients) stored in
C the common block /VALC/ as in the submitted versions of user-defined
C model specification FORTRAN77 source code files 'srfc.for', 'parm.for'
C  and 'val.for'.  Thus, unlike the other parts of the complete ray
C tracing, the INV1TT program cannot work with user's modifications of
C the subroutines SRFC1, SRFC2, PARM1, and PARM2.
C
C.......................................................................
C
C                                                    
C Description of data files:
C
C Main input data file read from the interactive device (*):
C (1) 'INV1TT'
C     'INV1TT' name of the input file described below.
C     Default: 'INV1TT'='inv1tt.dat'.
C
C                                                  
C Input data INV1TT:
C     This main data file contains the names of other input and output
C     files.  The data are read in by the list directed input (free
C     format).  In the list of input data below, each numbered paragraph
C     indicates the beginning of a new input operation (new READ
C     statement).  All input variables are of the type CHARACTER.  Only
C     the first 80 characters of the strings are significant.
C (1) 'MODEL','POINTS','FTT','DATA','INV1LOG',/
C     'MODEL'... Input data file containing the model parameters.
C             See the file 'model.for' of the package 'MODEL'.
C             Description of MODEL
C     'POINTS'... Input data file containing the coordinates of shot and
C             receiver points.  Its structure is described below.
C             Ignored (not opened) if the DATA filename below is
C             blank.
C             Description of file POINTS
C     'FTT'... Input data file containing the travel times from the
C             field measurements.  Its structure is described below.
C             Ignored (not opened) if the DATA filename below is
C             blank.
C             Description of file FTT
C     'DATA'..Output file containing the computed travel times and their
C             derivatives with respect to the model coefficients.
C             Its structure is described below.
C             If the filename is blank, no file with objective prior
C             information is generated.
C             Description of file DATA
C     'INV1LOG'... Output log file.  Just for your information.
C             Not opened and generated if the DATA filename above is
C             blank.
C             Description of file IN1LOG
C     /...    An obligatory slash at the end of line.
C     Default: 'MODEL'='model.dat', 'DATA'='data.out',
C             'INV1LOG'='inv1log.out'.
C (2) DIST,VPOWER,/
C     DIST... Maximum distance between the source or receiver point and
C             the initial or end point of a synthetic ray.
C     VPOWER... Velocity power for the computation of the travel-time
C             check sum.  If the VPOWER-th velocity power is expressed,
C             in all blocks of the model, in terms of explicit functions
C             of model coordinates, linearly homogeneous in their
C             coefficients (e.g. B-splines), the travel time minus its
C             check sum (see the log output file) should be zero within
C             rounding errors.  Otherwise, the check sum may have no
C             sense.
C             VPOWER=0: no check sum is evaluated and printed.
C     /...    An obligatory slash at the end of line.
C     Default: VPOWER=0.0.
C (3) Any times the following data:
C     'CRT-R','CRT-S','CRT-I'
C     'CRT-R'... File with the quantities stored along rays (see
C             C.R.T.5.5.1).
C             Description of file CRT-R
C     'CRT-S'... File with the quantities at the points of
C             intersection of rays with the specified surface at which
C             the receivers are situated for the case of two-point ray
C             tracing (see C.R.T.5.5.2).
C             If this filename is not blank, just the two-point rays
C             with minimum travel time at each receiver are considered.
C             If this filename is blank, all rays are taken into
C             account.
C             Attention:  All rays taken into account must start in some
C             of the specified sources and terminate in some of the
C             specified receivers, see the input file FTT.
C             Description of file CRT-S
C     'CRT-I'... File with the quantities at the initial points
C             of rays, corresponding to the above file rays or points of
C             intersection (see C.R.T.6.1).
C             Description of file CRT-I
C (4) / (a slash).
C Example of data INV1TT
C
C                                                  
C Input data file POINTS:
C     This data file contains the coordinates of shot and receiver
C     points.  The data are read in by the list directed input
C     (free format).  In the list of input data below, each numbered
C     paragraph indicates the beginning of a new input operation (new
C     READ statement).  The CHARACTER strings are explicitly mentioned
C     in this description.  Otherwise, if the first letter of the
C     symbolic name of the input variable is I-N, the corresponding
C     value in input data must be of the type INTEGER.  Otherwise, the
C     input parameter is of the type REAL.
C (1) Several strings terminated by / (a slash).
C (2) List of the sources and receivers: Any times the following data:
C (2.1) POINT,X1,X2,X3
C     POINT...CHARACTER*11 string containing the name of the source or
C             receiver point.
C     X1,X2,X3... Coordinates of the source or receiver point.
C (3) / (a slash) or the end of file.
C Example of data POINTS
C
C                                                     
C Input data file FTT:
C     This data file contains the travel times from the field
C     measurements.  The data are read in by the list directed input
C     (free format).  In the list of input data below, each numbered
C     paragraph indicates the beginning of a new input operation (new
C     READ statement).  The CHARACTER strings are explicitly mentioned
C     in this description.  Otherwise, if the first letter of the
C     symbolic name of the input variable is I-N, the corresponding
C     value in input data must be of the type INTEGER.  Otherwise, the
C     input parameter is of the type REAL.
C (1) Several strings terminated by / (a slash).
C (2) List of the travel times: Any times the following data (2.1):
C (2.1) SOURCE,RECEIVER,TFIELD,TERR
C     SOURCE..CHARACTER*11 string containing the name of the source
C             point.
C     RECEIVER... CHARACTER*11 string containing the name of the
C             receiver point.  The source and receiver points may be
C             mutually interchanged.
C     TFIELD..Travel time from a field measurement.
C     TERR... Error of the travel time from a field measurement.
C (3) / (a slash) or the end of file.
C Example of data FTT
C
C                                                    
C Output file DATA:
C (1) ND,NM,(INDM(I),I=1,NM)
C     ND...   The number of data (i.e. the number of equations).
C     NM...   Number of unknown model coefficients.
C     INDM... Indices of the model coefficients influencing the travel
C             times.  The indices correspond to the relative location in
C             the memory.
C             B-spline coefficients are listed in the same order as the
C             grid velocities in the file 'MODEL'.
C (2) ND-times the following data (2.1):
C (2.1) KD,RD,ED,WD,(GD(I),I=1,NM)
C     KD...   Index of the field travel time within the file 'ftt'.
C             The field travel times are indexed consecutively 1,2,3,...
C     RD...   Field travel time minus the computed synthetic travel
C             time.  In the case of multiple two-point rays, the first
C             arrival of them is considered.
C     ED...   Prior error of the above travel time difference.
C             Identical to TERR, file 'FTT', part (3).  I.e. the
C             synthetic travel time is assumed to be sufficiently
C             accurate.
C     WD...   Field travel time (stored for the purpose to assess a
C             posterior data misfit covariance weighting matrix).
C     GD...   Derivatives of the synthetic travel time with respect to
C             the model coefficients INDM.
C
C                                                 
C Output log file INV1LOG:
C (1) For each considered ray:
C (1.1) SOURCE,RECEIVER,TFIELD,TDIF,SDIST,RDIST,CHECK
C     SOURCE..Name of the source point.
C     RECEIVER... Name of the receiver point.
C     TFIELD..Travel time from a field measurement.
C     TDIF... Field travel time minus the minimum synthetic travel time.
C     SDIST...Distance between the source and the initial point of the
C             synthetic ray.
C     RDIST...Distance between the receiver and the end point of the
C             synthetic ray.
C     CHECK...Synthetic travel time minus the travel time resulting from
C             the derivatives of the theoretical travel time with
C             respect to the model coefficients.  This quantity should
C             not exceed in order the numerical error of the synthetic
C             travel time.
C             In this version defined just for the models described in
C             the terms of velocity.
C
C-----------------------------------------------------------------------
C
C Common block /VALC/:
      INCLUDE 'val.inc'
C     val.inc
C None of the storage locations of the common block are altered.
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C     Allocation of working arrays:
C
      INTEGER MP,MT
      PARAMETER (MP=1000)
      PARAMETER (MT=4000)
      CHARACTER*11 POINT(MP)
      INTEGER KS(MT),KR(MT)
      REAL COOR(3,MP),TFIELD(MT),TERR(MT),TTMIN(MT)
      EQUIVALENCE (KS    ,RAM        )
      EQUIVALENCE (KR    ,RAM(  MT+1))
      EQUIVALENCE (TFIELD,RAM(2*MT+1))
      EQUIVALENCE (TERR  ,RAM(3*MT+1))
      EQUIVALENCE (TTMIN ,RAM(4*MT+1))
      EQUIVALENCE (COOR  ,RAM(5*MT+1))
C
      INTEGER INDM(NPAR)
      REAL SUM(NPAR)
C     INDM... Indices of the unknown model parameters.
      EQUIVALENCE (SUM  , RAM(5*MT+3*MP+1))
C
C     Output data:  subjective prior information, working array in SOFT
      REAL CS(1)
C
C-----------------------------------------------------------------------
C
C     Filenames:
      CHARACTER*80 FILE2,FILE0,FILE1,FILE3,FILE4,FILE6
C
C     Logical unit numbers:
      INTEGER LU0,LU1,LU2,LU3,LU4,LU5,LU7,IU1,IU2
      PARAMETER (LU0=10)
      PARAMETER (LU1=11)
      PARAMETER (LU2=12)
      PARAMETER (LU3=13)
      PARAMETER (LU4=14)
      PARAMETER (LU5=15)
      PARAMETER (LU7=17)
C
C     Input data:
      CHARACTER*1  TEXT(10)
      CHARACTER*11 SRC,REC
      INTEGER NP,NT
      REAL DIST,DIST2,VPOWER
C     POINT...Names of the source and receiver points.
C     NP...   Number of source and receiver points.
C     NT...   Number of field travel times.
C     KS(I)...Index of the source point corresponding to the I-th field
C             travel time.
C     KR(I)...Index of the receiver point corresponding to the I-th
C             field travel time.
C     DIST... Maximum distance between the source or receiver point and
C             the initial or end point of a synthetic ray.
C     DIST2...DIST**2
C     VPOWER... Velocity power for the computation of the travel-time
C             check SUM.
C     COOR... Coordinates of the source or receiver points.
C     TFIELD..Field travel times.
C     TERR... Field travel time errors.
C
C     Output data:  variations of the synthetic travel time:
      INTEGER NSUM,NM
C     NM...   Number of the unknown model parameters.
C
C     Auxiliary storage locations:
      INTEGER IS,IR,IT,ND,IRAYTT,I,J,K
      REAL TT,TTAUX,TDIF,XI1,XI2,XI3,XE1,XE2,XE3,PI(6),PE(6)
      REAL SI,SE,RI,RE,SDIST,RDIST
C     IS...   Index of the source point.
C     IR...   Index of the receiver point.
C     IT...   Index of the field travel time.
C     ND...   Number of synthetic travel times corresponding to the
C             field travel times.
C     IRAYTT..Index of the last processed ray.
C     I,J,K...Temporary storage locations.
C     TTMIN...Minimum synthetic travel times corresponding to the
C             individual field travel times.
C     TT...   Synthetic travel time.
C     TTAUX...Temporary storage location.
C     TDIF... Field travel time minus the minimum synthetic travel time.
C     XI1,XI2,XI3,XE1,XE2,XE3... Coordinates of the initial and end
C             points of the last processed ray.
C     PI,PE...Slowness vectors at the initial and end points of the
C             last processed ray.
C     SI,SE,RI,RE... Squares of the distances between the source or
C             receiver points and the initial or end points of the ray.
C     SDIST...Distance between the source and the initial point of the
C             synthetic ray.
C     RDIST...Distance between the receiver and the end point of the
C             synthetic ray.
C
C.......................................................................
C
C     Opening data files and reading the input data:
C
C     Main input data file read from the interactive device (*):
      WRITE(*,'(A)') ' Enter the name of the main input data file: '
      FILE0='inv1.dat'
      READ(*,*) FILE0
      WRITE(*,'(A)') '+                                            '
C
C     Input data INV1TT:
      OPEN(LU0,FILE=FILE0,STATUS='OLD')
      FILE2='model.dat'
      FILE4='data.out'
      FILE6='inv1log.out'
      READ(LU0,*) FILE2,FILE0,FILE1,FILE4,FILE6
      VPOWER=0.
      READ(LU0,*) DIST,VPOWER
      DIST2=DIST*DIST
      OPEN(LU4,FILE=FILE2,STATUS='OLD')
      CALL MODEL1(LU4)
      CLOSE(LU4)
C
C     Number of unknown model parameters:
      CALL SOFT(2,0,0,0,0,0,0,0.,NM,INDM,CS,1,CS)
      WRITE(*,'(A,I4,A)') '+',NM,' model parameters'
C
C.......................................................................
C
C     Reading source and receiver points:
C
      OPEN(LU4,FILE=FILE0,STATUS='OLD')
      NP=-1
      READ(LU4,*,END=2) TEXT
    1 CONTINUE
        NP=NP+1
        IF(NP.GE.MP) THEN
C         INV1TT-01
          CALL ERROR('INV1TT-01: Too many source and receiver points')
        END IF
        POINT(NP+1)=' '
        READ(LU4,*,END=2) POINT(NP+1),(COOR(I,NP+1),I=1,3)
      IF(POINT(NP+1).NE.' ') GO TO 1
    2 CONTINUE
      CLOSE(LU4)
C
C     Reading field travel times:
C
      OPEN(LU4,FILE=FILE1,STATUS='OLD')
      NT=0
      READ(LU4,*,END=2) TEXT
    3 CONTINUE
        NT=NT+1
        IF(NT.GT.MT) THEN
C         INV1TT-02
          CALL ERROR('INV1TT-02: Too many field travel times')
        END IF
        SRC=' '
        READ(LU4,*,END=9) SRC,REC,TFIELD(NT),TERR(NT)
        IF(SRC.EQ.' ') THEN
          GO TO 9
        END IF
        DO 4 I=1,NP
          IF(SRC.EQ.POINT(I)) THEN
            KS(NT)=I
            GO TO 5
          END IF
    4   CONTINUE
C       INV1TT-03
        CALL ERROR('INV1TT-03: Source name not recognized')
    5   CONTINUE
        DO 6 I=1,NP
          IF(REC.EQ.POINT(I)) THEN
            KR(NT)=I
            GO TO 7
          END IF
    6   CONTINUE
C       INV1TT-04
        CALL ERROR('INV1TT-04: Receiver name not recognized')
    7   CONTINUE
      GO TO 3
    9 CONTINUE
      NT=NT-1
      CLOSE(LU4)
C
C.......................................................................
C
C     Computing quantities describing objective prior information:
C
      OPEN(LU5,STATUS='SCRATCH',FORM='UNFORMATTED')
      OPEN(LU7,FILE=FILE6)
      WRITE(*,*)
C
      KS(NT+1)=NP+1
      KR(NT+1)=NP+1
      TFIELD(NT+1)=0.
      IRAY=0
      IWAVE=0
      NSUM=IPAR(IPAR(IPAR(2)))
      DO 12 I=1,NT
        TTMIN(I)=999999.
   12 CONTINUE
C
C     Loop for the files with computed rays
   20 CONTINUE
        FILE1=' '
        READ(LU0,*,END=70) FILE1,FILE2,FILE3
        IF(FILE1.EQ.' ') THEN
          GO TO 70
        END IF
        I=INDEX(FILE1,'          ')
        J=INDEX(FILE2,'          ')
        K=INDEX(FILE3,'          ')
        WRITE(*,'(''+Processing: '',3A)') FILE1(1:I),FILE2(1:J),
     *                                                        FILE3(1:K)
        WRITE(*,*)
        OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD')
        IF(FILE2.EQ.' ') THEN
          IU1=0
          IU2=LU1
        ELSE
          IU1=LU1
          IU2=LU2
          OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD')
        END IF
        OPEN(LU3,FILE=FILE3,FORM='UNFORMATTED',STATUS='OLD')
        IRAYTT=0
C
C       Loop for the points of intersection of rays with the surface
   30   CONTINUE
C         Reading the results of the complete ray tracing
          CALL AP00(IU1,IU2,LU3)
          IF(IPT.LE.1.AND.IRAYTT.NE.0)THEN
C           New ray - recording results for the last ray iraytt
C           loop for field travel times - searching for two-point ray
            DO 39 IT=1,NT
              IS=KS(IT)
              IR=KR(IT)
              SI=(COOR(1,IS)-XI1)**2+(COOR(2,IS)-XI2)**2
     *                              +(COOR(3,IS)-XI3)**2
              RI=(COOR(1,IR)-XI1)**2+(COOR(2,IR)-XI2)**2
     *                              +(COOR(3,IR)-XI3)**2
              SE=(COOR(1,IS)-XE1)**2+(COOR(2,IS)-XE2)**2
     *                              +(COOR(3,IS)-XE3)**2
              RE=(COOR(1,IR)-XE1)**2+(COOR(2,IR)-XE2)**2
     *                              +(COOR(3,IR)-XE3)**2
              IF(SE.LE.SI.AND.RI.LE.RE) THEN
C               interchanging source and receiver points
                SI=SE
                RE=RI
                IS=KR(IT)
                IR=KS(IT)
              END IF
              IF((SI.LE.DIST2.AND.RE.LE.DIST2)) THEN
C               Synthetic ray may correspond to the field travel time
C               check for ray directions near the source
C               (allowable angle deviation +-30 deg: cosine**2=0.750)
C               (allowable angle deviation +-15 deg: cosine**2=0.933)
                TTAUX=(COOR(1,IR)-COOR(1,IS))*(XE1-XI1)
     *               +(COOR(2,IR)-COOR(2,IS))*(XE2-XI2)
     *               +(COOR(3,IR)-COOR(3,IS))*(XE3-XI3)
                IF(TTAUX.GT.0..AND.TTAUX*TTAUX.GT.
     *                         0.933*((COOR(1,IR)-COOR(1,IS))**2
     *                               +(COOR(2,IR)-COOR(2,IS))**2
     *                               +(COOR(3,IR)-COOR(3,IS))**2)
     *                  *((XE1-XI1)**2+(XE2-XI2)**2+(XE3-XI3)**2) ) THEN
C                 Synthetic ray corresponds to the field travel time
                  TTAUX=TT-PI(1)*(COOR(1,IS)-XI1)
     *                    -PI(2)*(COOR(2,IS)-XI2)
     *                    -PI(3)*(COOR(3,IS)-XI3)
     *                    +PE(1)*(COOR(1,IR)-XE1)
     *                    +PE(2)*(COOR(2,IR)-XE2)
     *                    +PE(3)*(COOR(3,IR)-XE3)
                  IF(TTAUX.LT.TTMIN(IT)) THEN
                    TTMIN(IT)=TTAUX
C                   Possible minimum synthetic travel time
                    SDIST=SQRT(SI)
                    RDIST=SQRT(RE)
                    WRITE(LU5) IT,TT,TTAUX,SDIST,RDIST,(SUM(I),I=1,NSUM)
                  END IF
                END IF
              END IF
   39       CONTINUE
            IF(NT.EQ.0) THEN
              WRITE(LU5) NT+1,TT,TT,0.,0.,(SUM(I),I=1,NSUM)
            END IF
            IRAYTT=0
          END IF
          IF(IWAVE.EQ.0) THEN
            GO TO 60
          END IF
C         *** for future extensions (selection of two-point rays):
C         IF(IU1.NE.0) THEN
C           CALL AP30(IREC)
C           IF(IREC.EQ.0) THEN
C             IF(IPT.LE.1.) THEN
C               WRITE(*,'(''+WAVE:'',I3,''    RAY:'',I4,''    POINT:'',
C    *                                              I4)') IWAVE,IRAY,IPT
C             END IF
C             GO TO 30
C           END IF
C         END IF
C         ***
          IF(IPT.EQ.1.OR.MOD(IPT,10).EQ.0) THEN
            WRITE(*,'(''+WAVE:'',I3,''    RAY:'',I4,''    POINT:'',I4)')
     *                                                    IWAVE,IRAY,IPT
          END IF
          IRAYTT=IRAY
          XI1=YI(3)
          XI2=YI(4)
          XI3=YI(5)
          XE1=Y(3)
          XE2=Y(4)
          XE3=Y(5)
          CALL AP01(TT,TTAUX)
          CALL AP02(PI,PE)
          CALL AP29(NSUM,SUM)
        GO TO 30
C       End of the loop for points of intersection of rays with surface
C
   60   CONTINUE
        CLOSE(LU1)
        CLOSE(LU2)
      GO TO 20
C
C     All minimum travel times and their derivatives are stored in the
C     scratch file LU5.
C
C.......................................................................
C
C     Writing objective prior information:
C
   70 CONTINUE
      OPEN(LU4,FILE=FILE4)
      I=MAX0(INDEX(FILE4,'          ')-1,11)
      WRITE(*,'(''+Generating the output: '',A)') FILE4(1:I)
      ND=0
      DO 71 I=1,NT
        IF(TTMIN(I).LT.999999.) THEN
          ND=ND+1
        END IF
   71 CONTINUE
C     ND is the number of equations.
C
      REWIND(LU5)
      WRITE(LU4,'(I9,5I13)') ND,NM
      WRITE(LU4,'(I9,5I13)') (INDM(I),I=1,NM)
      WRITE(LU7,'(2A)') ' SOURCE      RECEIVER       TFIELD    ',
     *                  'TFIELD-TT     SDIST       RDIST    TT-CHECKSUM'
   73 CONTINUE
        READ(LU5,END=79) IT,TT,TTAUX,SDIST,RDIST,(SUM(I),I=1,NSUM)
        TTMIN(NT+1)=TTAUX
        IF(TTAUX.EQ.TTMIN(IT)) THEN
C         Minimum synthetic travel time
C
C         System of linear equations:
          TDIF=TFIELD(IT)-TTAUX
          WRITE(LU4,'(I9,4X,6G13.6)') IT,TDIF,TERR(IT),TFIELD(IT)
          WRITE(LU4,'(6G13.6)') (SUM(INDM(I)),I=1,NM)
C
C         Check sums and log output:
          IF(VPOWER.NE.0.) THEN
            TTAUX=0.
            DO 74 I=1,NM
              J=INDM(I)
              IF(IPAR(IPAR(IPAR(1))).LT.J) THEN
                IF(SUM(J).NE.0.) THEN
                  TTAUX=TTAUX+RPAR(J)*SUM(J)
                END IF
              END IF
   74       CONTINUE
            IS=KS(IT)
            IR=KR(IT)
            TTAUX=TT+VPOWER*TTAUX
            WRITE(LU7,'(2(1X,A),5F12.6)') POINT(IS),POINT(IR),
     *           TFIELD(IT),TDIF,SDIST,RDIST,TTAUX
          ELSE
            WRITE(LU7,'(2(1X,A),5F12.6)') POINT(IS),POINT(IR),
     *           TFIELD(IT),TDIF,SDIST,RDIST
          END IF
        END IF
      GO TO 73
C
   79 CONTINUE
      CLOSE(LU4)
      CLOSE(LU5)
      CLOSE(LU7)
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'modelv.for'
C     modelv.for
      INCLUDE 'metric.for'
C     metric.for
      INCLUDE 'srfc.for'
C     srfc.for
      INCLUDE 'parmv.for'
C     parmv.for
      INCLUDE 'valv.for'
C     valv.for
      INCLUDE 'fitv.for'
C     fitv.for
      INCLUDE 'var.for'
C     var.for
      INCLUDE 'spsp.for'
C     spsp.for
      INCLUDE 'soft.for'
C     soft.for
      INCLUDE 'means.for'
C     means.for
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'apvar.for'
C     apvar.for
C
C=======================================================================
C