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 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 AP29A...Auxiliary subroutine to AP29.
C
C Date: 1994, January 15
C Coded by Ludek Klimes
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 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...Weigting 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
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 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 Error messages:
C 729... Array index out of range:
C Dimension MFUN of arrays IFUN1, FUN1, IFUN2, FUN2 should
C be increased.
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
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
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
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
PAUSE 'Error 729 in AP29: Array index out of range'
STOP
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