C 
C Subroutine file 'apvar.for': Applications and processing of the
C results of complete ray tracing --- Part2: Travel-time variations
C
C By Ludek Klimes
C
C This file consists of the following external procedures:
C     AP28... Subroutine designed to perform the numerical quadrature of
C             the set of given functions along a ray.  It has to be
C             called once at each point along the ray in which the
C             computed quantities are stored, i.e. after each invocation
C             of the subroutine AP00 which reads the quantities into the
C             common block /POINTC/.
C             AP28
C     AP29... Subroutine designed to evaluate the variations of the
C             travel time with respect to the model coefficients.
C             It has to be called once at each point along the ray at
C             which the computed quantities are stored, i.e. after each
C             invocation of the subroutine AP00 which reads the
C             quantities into the common block /POINTC/.
C             AP29
C     AP29A...Auxiliary subroutine to AP29.
C             AP29A
C
C Date: 1994, January 15
C Coded by Ludek Klimes
C
C=======================================================================
C
C     
C
      SUBROUTINE AP28(NSUM,SUM,IX,NDER,STEP,
     *                                NFUN1,IFUN1,FUN1,NFUN2,IFUN2,FUN2)
      INTEGER NSUM,IX,NDER,NFUN1,IFUN1(*),NFUN2,IFUN2(NFUN2)
      REAL SUM(NSUM),STEP,FUN1(*),FUN2(NFUN2*NDER)
C
C This subroutine performs the numerical quadrature of the set of given
C functions along a whole ray.  It has to be called once at each point
C along the ray in which the computed quantities are stored, i.e. after
C each invocation of the subroutine AP00 which reads the quantities into
C the common block /POINTC/.
C
C Input:
C     NSUM... Total number of the functions to be numerically integrated
C             along the ray.
C     SUM...  Array of dimension at least NSUM, in which the integrals
C             of the given functions are accumulated.  Its elements are
C             set to zeros at the initial point of the ray by this
C             subroutine.
C     IX...   Specifies the independent variable along the ray:
C             IX=0 independent variable is X, i.e. the same as for the
C               ray tracing.
C             IX=1 independent variable is the travel time.
C     NDER... NDER=1 if just the functional values of the integrated
C               functions are submitted.  Then the relative error of the
C               numerical quadrature is proportional to the third power
C               of the step along the ray (see the parameter store in
C               the input data (3) for the file 'ray.for').
C               When integrating a B-spline in a regular grid, the error
C               is about 0.01 for the step of half the size of the grid
C               interval.
C             NDER=2 if both the functional values and first derivatives
C               of the integrated functions are submitted.  Then the
C               relative error of the numerical quadrature is
C               proportional to the fourth power of the step along the
C               ray (see the parameter STORE in the input data (3) for
C               the file 'ray.for').
C               When integrating a B-spline in a regular grid, the error
C               is about 0.01 for the step of the size of the grid
C               interval.
C     STEP... Step in the independent variable along the ray (see the
C             parameter STORE in the input data (3) for the file
C             'ray.for').  Required just if NDER=1.  If NDER=1 and STEP
C             has not the correct value, the relative error of the
C             numerical quadrature is proportional to the second power
C             of the actual step along the ray.  When integrating a
C             B-spline with NDER=1 and STEP=0, in a regular grid, the
C             error is about 0.01 for the step of the size of 0.4 grid
C             interval.
C     NFUN1...Number of functions having nonzero values (or nonzero
C             first derivatives if NDER=2) at the previous point along
C             the ray.
C     IFUN1...IFUN(1:NFUN1)... Indices in the array SUM corresponding to
C             the functions having nonzero values (or nonzero first
C             derivatives if NDER=2) at the previous point along the
C             ray.
C     FUN1... FUN(1:NFUN1)... Values of the functions having nonzero
C               values (or nonzero first derivatives if NDER=2) at the
C               previous point along the ray.
C             FUN(NFUN1+1:2*NFUN1)... For NDER=2, first derivatives with
C               respect to the independent variable along the ray at the
C               previous point along the ray, of the functions having
C               nonzero values or nonzero first derivatives.
C             If this subroutine is invoked at the first point after the
C             initial point of the ray, the input values of NFUN1, IFUN
C             and FUN1 correspond to the initial (zero) point of the
C             ray.  At the subsequent points along the ray, the input
C             values of NFUN1, IFUN and FUN1 are the output from the
C             previous invocation of this subroutine.
C     NFUN2...Number of functions having nonzero values (or nonzero
C             first derivatives if NDER=2) at the current point along
C             the ray.
C     IFUN2...Indices in the array SUM corresponding to the functions
C             having nonzero values (or nonzero first derivatives if
C             NDER=2) at the current point along the ray.
C     FUN2... FUN(1:NFUN2)... Values of the functions having nonzero
C               values (or nonzero first derivatives if NDER=2) at the
C               current point along the ray.
C             FUN(NFUN2+1:2*NFUN2)... For NDER=2, first derivatives with
C               respect to the independent variable along the ray at the
C               previous point along the ray, of the functions having
C               nonzero values or nonzero first derivatives.
C
C Output:
C     SUM...  Integrals of the given functions with respect to the
C             independent variable along the ray, from the initial point
C             of the ray to the current point along the ray (stored in
C             the common block /POINTC/).
C     NFUN1,IFUN1,FUN1... Copies of the input values of NFUN2, IFUN2 and
C             FUN2.
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C No subroutines and external functions required.
C
C Date: 1994, January 23
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
C
      REAL X1,X2,W1,W2
      INTEGER ISRF1,ISRF2,ISUM,I
      SAVE X1,ISRF1
C
C     X1...   Independent variable along the ray at the previous point.
C     X2...   Independent variable along the ray at the current point.
C     ISRF1...Index of the surface at which the previous point is
C             situated, zero inside a complex block.
C     ISRF2...Index of the surface at which the current point is
C             situated, zero inside a complex block.
C     W1,W2...Weighting coefficients of the numerical quadrature.
C     ISUM... Index in the array SUM corresponding to the function under
C             consideration.
C     I...    Loop variable.
C
C.......................................................................
C
C     Integrals are set to zeros at the initial point of the ray
      IF(IPT.LE.1) THEN
        DO 10 ISUM=1,NSUM
          SUM(ISUM)=0.
   10   CONTINUE
        ISRF1=1
        IF(IX.LE.0) THEN
          X1=0.
        ELSE
          X1=YI(IX)
        END IF
      END IF
      IF(NYF.GT.0) THEN
        ISRF2=ISRFF
        IF(IX.LE.0) THEN
          X2=XF
        ELSE
          X2=YF(IX)
        END IF
      ELSE
        ISRF2=ISRF
        IF(IX.LE.0) THEN
          X2=X
        ELSE
          X2=Y(IX)
        END IF
      END IF
C
C     Numerical quadrature
      IF(X2.NE.X1) THEN
        W1=(X2-X1)/2.
        W2=W1
        IF(NDER.EQ.1) THEN
          IF(ISRF1.NE.0) THEN
            IF(ISRF2.EQ.0) THEN
C             First interval of the ray element
              W1=W1-(STEP*STEP)/(12.*(X2-X1))
              W2=W2+(STEP*STEP)/(12.*(X2-X1))
            END IF
          ELSE
            IF(ISRF2.NE.0) THEN
C             Last interval of the ray element
              W1=W1+(STEP*STEP)/(12.*(X2-X1))
              W2=W2-(STEP*STEP)/(12.*(X2-X1))
            END IF
          END IF
        END IF
        DO 21 I=1,NFUN1
          ISUM=IFUN1(I)
          SUM(ISUM)=SUM(ISUM)+W1*FUN1(I)
   21   CONTINUE
        DO 22 I=1,NFUN2
          ISUM=IFUN2(I)
          SUM(ISUM)=SUM(ISUM)+W2*FUN2(I)
   22   CONTINUE
        IF(NDER.EQ.2) THEN
          W1=((X2-X1)**2)/12.
          W2=-W1
          DO 31 I=1,NFUN1
            ISUM=IFUN1(I)
            SUM(ISUM)=SUM(ISUM)+W1*FUN1(NFUN1+I)
   31     CONTINUE
          DO 32 I=1,NFUN2
            ISUM=IFUN2(I)
            SUM(ISUM)=SUM(ISUM)+W2*FUN2(NFUN2+I)
   32     CONTINUE
        END IF
      END IF
C
C     Copying NFUN2,IFUN2,FUN2 into NFUN1,IFUN1,FUN1
      NFUN1=NFUN2
      DO 91 I=1,NFUN2
        IFUN1(I)=IFUN2(I)
        FUN1(I)=FUN2(I)
   91 CONTINUE
      DO 92 I=NFUN2+1,NDER*NFUN2
        FUN1(I)=FUN2(I)
   92 CONTINUE
C
      X1=X2
      ISRF1=ISRF2
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE AP29(NSUM,SUM)
      INTEGER NSUM
      REAL SUM(NSUM)
C
C This subroutine evaluates variations of the travel time with respect
C to the model coefficients.  It has to be called once at each point
C along the ray at which the computed quantities are stored, i.e. after
C each invocation of the subroutine AP00 which reads the quantities into
C the common block /POINTC/.
C Subroutine PARM2 is called to evaluate the material parameters at the
C current point and, at a structural interface, also subroutine SRFC2 is
C called to evaluate the function describing the interface.  After the
C invocation of PARM2 or SRFC2, respectively, subroutine VAR6 is called
C to recall the variations of the model parameters or of the interface,
C with respect to the model coefficients.  If the user replaces the
C subroutine file 'parm.for' or 'srfc.for' by his own version, it is his
C own responsibility to call subroutines VAR1 to VAR5 (see the file
C 'var.for') in such a way that the required variations are stored when
C returning from his own subroutine PARM2 or SRFC2.
C
C Input:
C     NSUM... Total number of the coefficients describing the model.
C     SUM...  Array of dimension at least NSUM, in which the variations
C             of the travel time with respect to the model coefficients
C             are accumulated.  Its elements are set to zeros at the
C             initial point of the ray by this subroutine.
C
C Output:
C     SUM...  Variations of the travel time (from the initial point of
C             the ray to the current point along the ray) with respect
C             to the model coefficients.
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      INTEGER KOOR
      EXTERNAL KOOR,METRIC,SRFC2,VAR6,AP28
C     SMVPRD
C     PARM2,VELOC
C
C Date: 1996, September 30
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations for local model parameters:  FAUX(10),
C       G(12),GAMMA(18),GSQRD,  UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL:
      INCLUDE 'auxmod.inc'
C     auxmod.inc
C
C.......................................................................
C
C     Auxiliary storage locations:
C
      INTEGER ISR,II,IBI,MFUN
      PARAMETER (MFUN=64)
      INTEGER NFUN1,IFUN1(MFUN),NFUN2,IFUN2(MFUN)
      REAL B0I,B1I,B2I,B3I,FUN1(2*MFUN),FUN2(2*MFUN)
      REAL PIN1,PIN2,PIN3,C(3),P1,P2,P3,AUX0
      SAVE NFUN1,IFUN1,FUN1,PIN1,PIN2,PIN3
C
C     ISR...  Index of the surface covering the interface.
C     II...   Loop variable (sequential number of the required
C             variation).
C     IBI...  Absolute index of the function coefficient.
C     B0I,B1I,B2I,B3I... Variation of the functional value and the three
C             first derivatives, with respect to the IBI-th coefficient
C             of the model.
C     NFUN1...Number of functions having nonzero values or nonzero first
C             derivatives at the previous point along the ray.
C     IFUN1...IFUN(1:NFUN1)... Indices in the array SUM corresponding to
C             the functions having nonzero values or nonzero first
C             derivatives at the previous point along the ray.
C     FUN1... FUN(1:NFUN1)... Values of the functions having nonzero
C               values or nonzero first derivatives at the previous
C               point along the ray.
C             FUN(NFUN1+1:2*NFUN1)... First derivatives with respect to
C               the independent variable along the ray at the previous
C               point along the ray, of the functions having nonzero
C               values or nonzero first derivatives.
C             At the first point after the initial point of the ray,
C             the values of NFUN1, IFUN and FUN1 correspond to the
C             initial (zero) point of the ray.
C     NFUN2...Number of functions having nonzero values or nonzero
C             first derivatives at the current point along the ray.
C     IFUN2...Indices in the array SUM corresponding to the functions
C             having nonzero values or nonzero first derivatives at the
C             current point along the ray.
C     FUN2... FUN(1:NFUN2)... Values of the functions having nonzero
C               values or nonzero first derivatives at the current point
C               along the ray.
C             FUN(NFUN2+1:2*NFUN2)... First derivatives with respect to
C               the independent variable along the ray at the previous
C               point along the ray, of the functions having nonzero
C               values or nonzero first derivatives.
C     PIN1,PIN2,PIN3... Contravariant components of the slowness vector
C             at the point of incidence.
C     C...    Coordinates.
C     P1,P2,P3... Contravariant components of the slowness vector.
C     AUX0... Temporary storage location.
C
C.......................................................................
C
C     Initial point of the ray:
      IF(IPT.LE.1) THEN
        PIN1=0.
        PIN2=0.
        PIN3=0.
        CALL AP29A(ICB1I,0,YI,ISR,C,P1,P2,P3,MFUN,NFUN1,IFUN1,FUN1)
      END IF
C
C     Another point of the ray:
      IF(NYF.GT.0) THEN
        CALL AP29A(ICB1F,ISRFF,YF,ISR,C,P1,P2,P3,MFUN,NFUN2,IFUN2,FUN2)
      ELSE
        CALL AP29A(ICB1, ISRF ,Y ,ISR,C,P1,P2,P3,MFUN,NFUN2,IFUN2,FUN2)
      END IF
C
C     Numerical quadrature:
      CALL AP28(NSUM,SUM,1,2,0.,NFUN1,IFUN1,FUN1,NFUN2,IFUN2,FUN2)
C
C     Structural interface:
      IF(ISR.NE.0) THEN
        IF(PIN1.EQ.0..AND.PIN2.EQ.0..AND.PIN3.EQ.0.) THEN
C         incident ray:
          PIN1=P1
          PIN2=P2
          PIN3=P3
        ELSE
C         Reflected/transmitted ray:
C         Including the variation of the travel time with respect to the
C         structural interface
          CALL SRFC2(IABS(ISR),C,VD)
          IF(KOOR().NE.0) THEN
            CALL METRIC(C,GSQRD,G,GAMMA)
            AUX0=VD(2)*(G(7)*VD(2)+2.*(G(8)*VD(3)+G(10)*VD(4))) +
     *           VD(3)*(G(9)*VD(3)+2.*G(11)*VD(4)) + VD(4)*G(12)*VD(4)
          ELSE
            AUX0=VD(2)*VD(2)+VD(3)*VD(3)+VD(4)*VD(4)
          END IF
          AUX0=( VD(2)*(P1-PIN1)+VD(3)*(P2-PIN2)+VD(4)*(P3-PIN3) )/AUX0
          II=0
   30     CONTINUE
            II=II+1
            CALL VAR6(1,II,NFUN2,IBI,B0I,B1I,B2I,B3I)
            IF(II.LE.NFUN2) THEN
              SUM(IBI)=SUM(IBI)+AUX0*B0I
            END IF
          IF(II.LT.NFUN2) GO TO 30
          PIN1=0.
          PIN2=0.
          PIN3=0.
        END IF
      END IF
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C
C     
C
      SUBROUTINE AP29A(ICB1,ISRF,Y,ISR,C,P1,P2,P3,MFUN,NFUN,IFUN,FUN)
      INTEGER ICB1,ISRF,ISR,MFUN,NFUN,IFUN(MFUN)
      REAL Y(8),C(3),P1,P2,P3,FUN(2*MFUN)
C
C Auxiliary subroutine to AP29.
C
C Input:
C     ICB1... Index of the complex block.
C     ISRF... Index of the surface covering the interface.
C     Y...    Quantities computed along a ray.
C     MFUN... Array dimension.
C
C Output:
C     ISR...  Index of the surface covering the interface.
C     C...    Coordinates.
C     P1,P2,P3... Contravariant components of the slowness vector.
C     NFUN... Number of variations.
C     IFUN... Indices of variations.
C     FUN...  FUN(1:NFUN)... Values of variations.
C             FUN(NFUN+1:2*NFUN)... First derivatives of variations
C               with respect to the independent variable along the ray.
C
C Subroutines and external functions required:
      INTEGER KOOR
      EXTERNAL KOOR,METRIC,SMVPRD,PARM2,VELOC,VAR6
C
C.......................................................................
C
C     Auxiliary storage locations for local model parameters:  FAUX(10),
C       G(12),GAMMA(18),GSQRD,  UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL:
      INCLUDE 'auxmod.inc'
C     auxmod.inc
C
C.......................................................................
C
      REAL AUX0,AUX1,AUX2,AUX3,AUX4
      INTEGER NEXPS,IVAL,II
      PARAMETER (NEXPS=0)
      REAL B0I,B1I,B2I,B3I
C
C     AUX0,AUX1,AUX2,AUX3,AUX4... Auxiliary storage locations
C             for local model parameters or temporary variables.
C     IVAL... Index of the function describing the model.
C             IVAL=1 for P-wave,
C             IVAL=2 for S-wave.
C     II...   Loop variable (sequential number of the required
C             variation).
C     B0I,B1I,B2I,B3I... Variation of the functional value and the three
C             first derivatives, with respect to the IBI-th coefficient
C             of the model.
C
C.......................................................................
C
C     Assignments:
      ISR=ISRF
      C(1)=Y(3)
      C(2)=Y(4)
      C(3)=Y(5)
      IF(KOOR().NE.0) THEN
        CALL METRIC(Y(3),GSQRD,G,GAMMA)
        CALL SMVPRD(G(7),Y(6),Y(7),Y(8),P1,P2,P3)
      ELSE
        P1=Y(6)
        P2=Y(7)
        P3=Y(8)
      END IF
C
C     Material parameters:
      CALL PARM2(IABS(ICB1),Y(3),UP,US,AUX0,AUX1,AUX2)
      CALL VELOC(ICB1,UP,US,AUX1,AUX2,AUX3,AUX4,VD,AUX0)
C     Material parameters and their variations are defined.
C
      AUX0=-VD(1)**(-NEXPS-1)
      AUX4=-VD(1)**(-NEXPS-NEXPS+1)
      AUX1=AUX4*P1
      AUX2=AUX4*P2
      AUX3=AUX4*P3
      AUX4=-FLOAT(NEXPS+1)*(VD(2)*AUX1+VD(3)*AUX2+VD(4)*AUX3)/VD(1)
      IF(ICB1.GT.0) THEN
C       P-wave:
        IVAL=1
      ELSE
C       S-wave:
        IVAL=2
      END IF
      II=0
   20 CONTINUE
        II=II+1
        CALL VAR6(IVAL,II,NFUN,IFUN(II),B0I,B1I,B2I,B3I)
        IF(II.LE.NFUN) THEN
          IF(NFUN.GT.MFUN) THEN
C           729
            PAUSE 'Error 729 in AP29: Array index out of range'
            STOP
C           Dimension MFUN of arrays IFUN1, FUN1, IFUN2, FUN2 should
C           be increased.
          END IF
          FUN(II)=AUX0*B0I
          FUN(NFUN+II)=AUX1*B1I+AUX2*B2I+AUX3*B3I+AUX4*B0I
        END IF
      IF(II.LT.NFUN) GO TO 20
      RETURN
      END
C
C=======================================================================
C