C                                                            <HTML><PRE>
C Program HDERTEST to test the calculation of the first-order and
C second-order spatial derivatives of the Hamiltonian function
C calculated by subroutine HDER of file <A HREF="hder.for" TYPE="text/html">hder.for</A>
C
C Program HDERTEST reads the complex-valued density-reduced stiffness
C tensor composed of complex-valued elastic moduli,
C and its first-order and second-order partial derivatives with respect
C to 3 spatial coordinates.
C Program HDERTEST also reads a given slowness vector.
C For the given density-reduced stiffness tensor and slowness vector,
C program HDERTEST calculates the first-order and second-order
C phase-space derivatives of the Hamiltonian function using subroutines
C HDERR and HDER.  The kind of Hamiltonian function is specified by
C input SEP <A HREF="hder.for#SEP" TYPE="text/html">parameters</A> of subroutine HDER.
C The given density-reduced stiffness tensor is then extrapolated
C using the Taylor expansion of the second order.  The first-order
C phase-space derivatives of the Hamiltonian function are estimated
C from the values of the Hamiltonian function using the symmetric finite
C differences.  The second-order phase-space derivatives of the
C Hamiltonian function are estimated from the values of the first-order
C phase-space derivatives of the Hamiltonian function using the
C symmetric finite differences.
C The first-order and second-order phase-space derivatives of the
C Hamiltonian function obtained by subroutine HDER are then compared
C with the same derivatives obtained numerically by finite differences.
C
C Version: 8.10
C Date: 2022, November 2
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/klimes.htm
C
C.......................................................................
C                                                    <A NAME="data"></A>
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     <A NAME="SEP"></A>
C Input data file 'SEP':
C     File 'SEP' has the form of the <A HREF="../forms/sep.htm">SEP</A>
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Data specifying input file:
C     ANIDER='string'... Name of the input file containing the
C             density-reduced stiffness tensor (elastic moduli)
C             and its first-order and second-order spatial derivatives.
C             <A HREF="#ANIDER">Description of file ANIDER</A>.
C             No default.  ANIDER must be specified and cannot be blank.
C Data specifying output file:
C     HDERTEST='string'... Name of the output file containing
C             the first-order and second-order phase-space derivatives
C             of the Hamiltonian function obtained by subroutine HDER,
C             and the same derivatives obtained numerically by finite
C             differences.
C             <A HREF="#HDERTEST">Description of file HDERTEST</A>.
C             Default: HDERTEST='hdertest.out'
C Given slowness vector:
C     P1=real, P2=real, P3=real... Components of the direction of
C             a given slowness vector.  The length of this vector may be
C             arbitrary.
C             Default: P1=0.0, P2=0.0, P3=1.0
C Data specifying the Hamiltonian function:
C     ICB=integer... Positive value specifies the P wave,
C             negative value specifies an S wave.
C             Default: ICB=-1
C     DEGREE=real... Degree of the homogeneous complex-valued
C             Hamiltonian function.
C             Simultaneously the degree of the homogeneous Hamiltonian
C             function to be arithmetically averaged for common S-wave
C             ray tracing.
C             DEGREE must equal 2, 1, -1 or -2.
C             Default: DEGREE=-1
C     <A HREF="hder.for#SEP" TYPE="text/html">Input SEP parameters of subroutine HDER.</A>
C Numerical parameters for finite differencing:
C     DX=real... Step in spatial coordinates for symmetric finite
C             differences.
C             Default: DX=0.002
C     DP=real... Step in slowness-vector components for symmetric finite
C             differences.
C             Default: DP=0.001
C     NFD=integer... Order of finite differences.
C             NFD=2, NFD=4 or NFD=6.
C             Default: NFD=4
C     NPOWER=integer... Number of subsequent finite differences with
C             steps DX and DP increased by factors 2, 4, ..., 2**NPOWER.
C             Default: NPOWER=0
C
C                                                  <A NAME="ANIDER"></A>
C Input file ANIDER with the density-reduced elastic moduli
C and their first-order and second-order spatial derivatives:
C (1) 21 real parts of the elastic moduli read at once in order
C     A1111,A1122,A1133,A1123,A1113,A1112,
C           A2222,A2233,A2223,A2213,A2212,
C                 A3333,A3323,A3313,A3312,
C                       A2323,A2313,A2312,
C                             A1313,A1312,
C                                   A1212
C (2) First-order derivatives of 21 real parts of the elastic moduli
C     with respect to the X1 coordinate axis read at once in the same
C     order as above.
C (3) First-order derivatives of 21 real parts of the elastic moduli
C     with respect to the X2 coordinate axis read at once in the same
C     order.
C (4) First-order derivatives of 21 real parts of the elastic moduli
C     with respect to the X3 coordinate axis read at once in the same
C     order.
C (5) Homogeneous second-order derivatives of 21 real parts of the
C     elastic moduli with respect to the X1 coordinate axis read at once
C     in the same order.
C (6) Mixed second-order derivatives of 21 real parts of the elastic
C     moduli with respect to the X1 and X2 coordinate axes read at once
C     in the same order.
C (7) Homogeneous second-order derivatives of 21 real parts of the
C     elastic moduli with respect to the X2 coordinate axis read at once
C     in the same order.
C (8) Mixed second-order derivatives of 21 real parts of the elastic
C     moduli with respect to the X1 and X3 coordinate axes read at once
C     in the same order.
C (9) Mixed second-order derivatives of 21 real parts of the elastic
C     moduli with respect to the X2 and X3 coordinate axes read at once
C     in the same order.
C (10) Homogeneous second-order derivatives of 21 real parts of the
C     elastic moduli with respect to the X3 coordinate axis read at once
C     in the same order.
C (11)-(20) Analogous to (1)-(10), but for the imaginary parts of the
C     elastic moduli.
C
C                                                <A NAME="HDERTEST"></A>
C Output file HDERTEST contains
C (1) Copy of input data  P1, P2, P3, ICB, DEGREE.
C (2) Determined real-valued referenc polarization vectors
C     E1j, E2j, E3j, j=1,2,3,
C     real parts HP and HS of P-wave and S-wave Hamiltonian functions,
C     and imaginary part HI of the Hamiltonian functions corresponding
C     to ICB, calculated by
C (2.1) subroutine HDERR,
C (2.2) subroutine HDER.
C (3) 6 header lines specifying the phase-space derivatives of the
C Hamiltonian function, each followed by corresponding data lines.
C Each header line 'H1  H1j  H2j  H3j  H4j  H5j  H6j', j=1,2,3,4,5,6,
C corresponds to symmetric finite differencing H, H1, H2, H3, H4, H5, H6
C with respect to Xj, where X1, X2, X3, X4=P1, X5=P2 and X6=P3 are the
C phase-space coordinates.  The first line below each header line
C contains the values calculated by subroutine HDERR.  The second line
C below each header line contains the values calculated by subroutine
C HDER.  The subsequent NPOWER+1 lines contain the values estimated by
C the symmetric finite differencing the output of subroutine HDER.
C
C=======================================================================
C
      EXTERNAL ERROR,RSEP1,RSEP3T,EIGEN
C
C     Filenames and parameters:
      CHARACTER*80 FILE
      INTEGER LU,MFD
      PARAMETER (LU=1,MFD=3)
C
C     Stiffness tensor with derivatives
      REAL A(10,21),B(10,21)
C     Imaginary part of the stiffness tensor with derivatives
      REAL Q(10,21),R(10,21)
C
C     Arguments of HDER:
      INTEGER ICB,IERR
      REAL P1,P2,P3,HP,HS,H,HI,H1,H2,H3,H4,H5,H6
      REAL H11,H12,H22,H13,H23,H33,H14,H24,H34,H44
      REAL H15,H25,H35,H45,H55,H16,H26,H36,H46,H56,H66,E(9)
      REAL Q1,Q2,Q3,      G,G1,G2,G3,G4,G5,G6
      REAL G11,G12,G22,G13,G23,G33,G14,G24,G34,G44
      REAL G15,G25,G35,G45,G55,G16,G26,G36,G46,G56,G66
      REAL                F,F1,F2,F3,F4,F5,F6
      REAL F11,F12,F22,F13,F23,F33,F14,F24,F34,F44
      REAL F15,F25,F35,F45,F55,F16,F26,F36,F46,F56,F66
C
C     Arguments of HDERR:
      REAL         EP,ES,EE,E1,E2,E3,E4,E5,E6
      REAL E11,E12,E22,E13,E23,E33,E14,E24,E34,E44
      REAL E15,E25,E35,E45,E55,E16,E26,E36,E46,E56,E66
C
C     Other variables:
      INTEGER I,J,K,IFD,NFD,NPOWER,IDER
      REAL DEGREE,DX,DP
      REAL X,D(MFD),D1(MFD),D2(MFD),D3(MFD),D4(MFD),D5(MFD),D6(MFD)
C
C-----------------------------------------------------------------------
C
C     Reading main input data:
      WRITE(*,'(A)') '+HDERTEST: Enter input filename: '
      FILE=' '
      READ(*,*) FILE
      IF(FILE.EQ.' ') THEN
C       <A NAME="error01"></A><FONT COLOR="RED">HDERTEST-01</FONT>
        CALL ERROR('HDERTEST-01: No input SEP file specified')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      END IF
      CALL RSEP1(LU,FILE)
      WRITE(*,'(A)') '+HDERTEST: Working...            '
C
C     Reading density reduced stiffness tensor (elastic moduli):
      CALL RSEP3T('ANIDER',FILE,' ')
      IF(FILE.EQ.' ') THEN
C       <A NAME="error02"></A><FONT COLOR="RED">HDERTEST-02</FONT>
        CALL ERROR('HDERTEST-02: No input stiffness tensor specified')
C       Name ANIDER of the input file with the stiffness tensor
C       and its first-order and second-order spatial derivatives
C       must be specified and cannot be blank.
      END IF
      OPEN(LU,FILE=FILE,STATUS='OLD')
      DO 11 I=1,10
        READ(LU,*) A(I,1),A(I,2),A(I,4),A(I, 7),A(I,11),A(I,16)
     *                   ,A(I,3),A(I,5),A(I, 8),A(I,12),A(I,17)
     *                          ,A(I,6),A(I, 9),A(I,13),A(I,18)
     *                                 ,A(I,10),A(I,14),A(I,19)
     *                                         ,A(I,15),A(I,20)
     *                                                 ,A(I,21)
   11 CONTINUE
      DO 13 I=1,10
C       Default: no attenuation
        DO 12 J=1,21
          Q(I,J)=0.
   12   CONTINUE
        READ(LU,*) Q(I,1),Q(I,2),Q(I,4),Q(I, 7),Q(I,11),Q(I,16)
     *                   ,Q(I,3),Q(I,5),Q(I, 8),Q(I,12),Q(I,17)
     *                          ,Q(I,6),Q(I, 9),Q(I,13),Q(I,18)
     *                                 ,Q(I,10),Q(I,14),Q(I,19)
     *                                         ,Q(I,15),Q(I,20)
     *                                                 ,Q(I,21)
   13 CONTINUE
      CLOSE(LU)
C
C     Reading the direction of the slowness vector:
      CALL RSEP3R('P1',P1,0.0)
      CALL RSEP3R('P2',P2,0.0)
      CALL RSEP3R('P3',P3,1.0)
C
C     Data specifying the Hamiltonian function:
      CALL RSEP3I('ICB',ICB,-1)
      CALL RSEP3R('DEGREE',DEGREE,-1)
C
C     Reading numerical parameters of the test:
      CALL RSEP3R('DX',DX,0.002)
      CALL RSEP3R('DP',DP,0.001)
      CALL RSEP3I('NFD',NFD,4)
      IF(NFD.NE.2.AND.NFD.NE.4.AND.NFD.NE.6) THEN
C       <A NAME="error03"></A><FONT COLOR="RED">HDERTEST-03</FONT>
        CALL ERROR('HDERTEST-03: Incorrect order NFD of differencing')
C       Order NFD of finite differencing must be NFD=2, NFD=4 or NFD=6.
      END IF
      CALL RSEP3I('NPOWER',NPOWER,0)
C
C     Opening the output file:
      CALL RSEP3T('HDERTEST',FILE,'hdertest.out')
      OPEN(LU,FILE=FILE)
C
C     Comparison with the elastic medium:
      CALL HDERR(ICB,A,P1,P2,P3,EP,ES,EE,E1,E2,E3,E4,E5,E6,
     *           E11,E12,E22,E13,E23,E33,E14,E24,E34,E44,
     *           E15,E25,E35,E45,E55,E16,E26,E36,E46,E56,E66,
     *           E,IERR)
      IF(IERR.NE.0) THEN
C       <A NAME="error04"></A><FONT COLOR="RED">HDERTEST-04</FONT>
        CALL ERROR('HDERTEST-04: Error in HDERR')
      END IF
C
C     Comparing HDERR and HDER:
      WRITE(LU,'(7(4X,A,4X))') 'P1 ','P2 ','P3 ','ICB','DEGREE'
      WRITE(LU,'(7F11.6)')      P1  , P2  , P3  , ICB , DEGREE
      WRITE(LU,'(7(4X,A,4X))') 'E1j','E2j','E3j','HP ','HS ','HI '
      WRITE(LU,'(7F11.6)')      E(1), E(2), E(3), EP  , ES  , 0.
      WRITE(LU,'(7F11.6)')      E(4), E(5), E(6)
      WRITE(LU,'(7F11.6)')      E(7), E(8), E(9)
C
C     Calculating phase-space derivatives of the Hamiltonian function:
      CALL HDER(ICB,A,Q,P1,P2,P3,HP,HS,H,HI,H1,H2,H3,H4,H5,H6,
     *          H11,H12,H22,H13,H23,H33,H14,H24,H34,H44,
     *          H15,H25,H35,H45,H55,H16,H26,H36,H46,H56,H66,
     *          E,IERR)
      IF(IERR.NE.0) THEN
C       <A NAME="error05"></A><FONT COLOR="RED">HDERTEST-05</FONT>
        CALL ERROR('HDERTEST-05: Error in HDER')
      END IF
C
C     Comparing HDERR and HDER:
      WRITE(LU,'(7F11.6)')      E(1), E(2), E(3), HP  , HS  , HI
      WRITE(LU,'(7F11.6)')      E(4), E(5), E(6)
      WRITE(LU,'(7F11.6)')      E(7), E(8), E(9)
C
C     Normalizing the given slowness vector:
      P1=P1/SQRT(H)
      P2=P2/SQRT(H)
      P3=P3/SQRT(H)
C
C     Spatial derivatives (comparing HDERR, HDER and finite diferences):
      DO 29 IDER=1,3
C       Writing the header line and the derivatives calculated by HDER:
        IF(IDER.EQ.1) THEN
C         X1 derivative
          WRITE(LU,'(7(4X,A,4X))')
     *                        'H1 ','H11','H21','H31','H41','H51','H61'
          WRITE(LU,'(7F11.6)') E1,E11,E12,E13,E14,E15,E16
          WRITE(LU,'(7F11.6)') H1,H11,H12,H13,H14,H15,H16
        ELSE IF(IDER.EQ.2) THEN
C         X2 derivative
          WRITE(LU,'(7(4X,A,4X))')
     *                        'H2 ','H12','H22','H32','H42','H52','H62'
          WRITE(LU,'(7F11.6)') E2,E12,E22,E23,E24,E25,E26
          WRITE(LU,'(7F11.6)') H2,H12,H22,H23,H24,H25,H26
        ELSE
C         X3 derivative
          WRITE(LU,'(7(4X,A,4X))')
     *                        'H3 ','H13','H23','H33','H43','H53','H63'
          WRITE(LU,'(7F11.6)') E3,E13,E23,E33,E34,E35,E36
          WRITE(LU,'(7F11.6)') H3,H13,H23,H33,H34,H35,H36
        END IF
C
C       Loop over finite differences:
        DO 27 K=0,NPOWER
          X=DX*2.**(1.0*FLOAT(K))
          DO 25 IFD=2,NFD,2
            IF(IFD.GT.2) THEN
              X=X+X
            END IF
            DO 22 J=1,21
              DO 21 I=1,10
                B(I,J)=A(I,J)
                B(I,J)=A(I,J)
   21         CONTINUE
              IF(IDER.EQ.1) THEN
C               X1 derivative
                B(1,J)=A(1,J)+A(2,J)*X+0.5*A(5,J)*X*X
                B(2,J)=A(2,J)+A(5,J)*X
                B(3,J)=A(3,J)+A(6,J)*X
                B(4,J)=A(4,J)+A(8,J)*X
                R(1,J)=Q(1,J)+Q(2,J)*X+0.5*Q(5,J)*X*X
                R(2,J)=Q(2,J)+Q(5,J)*X
                R(3,J)=Q(3,J)+Q(6,J)*X
                R(4,J)=Q(4,J)+Q(8,J)*X
              ELSE IF(IDER.EQ.2) THEN
C               X2 derivative
                B(1,J)=A(1,J)+A(3,J)*X+0.5*A(7,J)*X*X
                B(2,J)=A(2,J)+A(6,J)*X
                B(3,J)=A(3,J)+A(7,J)*X
                B(4,J)=A(4,J)+A(9,J)*X
                R(1,J)=Q(1,J)+Q(3,J)*X+0.5*Q(7,J)*X*X
                R(2,J)=Q(2,J)+Q(6,J)*X
                R(3,J)=Q(3,J)+Q(7,J)*X
                R(4,J)=Q(4,J)+Q(9,J)*X
              ELSE
C               X3 derivative
                B(1,J)=A(1,J)+A(4,J)*X+0.5*A(10,J)*X*X
                B(2,J)=A(2,J)+A(8,J)*X
                B(3,J)=A(3,J)+A(9,J)*X
                B(4,J)=A(4,J)+A(10,J)*X
                R(1,J)=Q(1,J)+Q(4,J)*X+0.5*Q(10,J)*X*X
                R(2,J)=Q(2,J)+Q(8,J)*X
                R(3,J)=Q(3,J)+Q(9,J)*X
                R(4,J)=Q(4,J)+Q(10,J)*X
              END IF
   22       CONTINUE
            IF(B(1,1).LT.0..OR.B(1,3).LT.0..OR.B(1,6).LT.0..OR.
     *         B(1,10).LT.0..OR.B(1,15).LT.0..OR.B(1,21).LT.0.) THEN
C             <A NAME="error06"></A><FONT COLOR="RED">HDERTEST-06</FONT>
              CALL ERROR('HDERTEST-06: Negative elastic modulus')
C             Negative elastic modulus A1111, A2222, A3333, A2323, A1313
C             or A1212 has been identified during finite differencing with
C             respect to one of spatial coordinates.
C             Check the relation of step DX in spatial coordinates in
C             comparison with the first-order and second-order derivatives
C             of the above mentioned elastic moduli.
            END IF
            CALL HDER(ICB,B,R,P1,P2,P3,HP,HS,G,HI,G1,G2,G3,G4,G5,G6,
     *                G11,G12,G22,G13,G23,G33,G14,G24,G34,G44,
     *                G15,G25,G35,G45,G55,G16,G26,G36,G46,G56,G66,
     *                E,IERR)
            IF(IERR.NE.0) THEN
C             <A NAME="error07"></A><FONT COLOR="RED">HDERTEST-07</FONT>
              CALL ERROR('HDERTEST-07: Error in HDER')
            END IF
            DO 23 J=1,21
              IF(IDER.EQ.1) THEN
C               X1 derivative
                B(1,J)=A(1,J)-A(2,J)*X+0.5*A(5,J)*X*X
                B(2,J)=A(2,J)-A(5,J)*X
                B(3,J)=A(3,J)-A(6,J)*X
                B(4,J)=A(4,J)-A(8,J)*X
                R(1,J)=Q(1,J)-Q(2,J)*X+0.5*Q(5,J)*X*X
                R(2,J)=Q(2,J)-Q(5,J)*X
                R(3,J)=Q(3,J)-Q(6,J)*X
                R(4,J)=Q(4,J)-Q(8,J)*X
              ELSE IF(IDER.EQ.2) THEN
C               X2 derivative
                B(1,J)=A(1,J)-A(3,J)*X+0.5*A(7,J)*X*X
                B(2,J)=A(2,J)-A(6,J)*X
                B(3,J)=A(3,J)-A(7,J)*X
                B(4,J)=A(4,J)-A(9,J)*X
                R(1,J)=Q(1,J)-Q(3,J)*X+0.5*Q(7,J)*X*X
                R(2,J)=Q(2,J)-Q(6,J)*X
                R(3,J)=Q(3,J)-Q(7,J)*X
                R(4,J)=Q(4,J)-Q(9,J)*X
              ELSE
C               X3 derivative
                B(1,J)=A(1,J)-A(4,J)*X+0.5*A(10,J)*X*X
                B(2,J)=A(2,J)-A(8,J)*X
                B(3,J)=A(3,J)-A(9,J)*X
                B(4,J)=A(4,J)-A(10,J)*X
                R(1,J)=Q(1,J)-Q(4,J)*X+0.5*Q(10,J)*X*X
                R(2,J)=Q(2,J)-Q(8,J)*X
                R(3,J)=Q(3,J)-Q(9,J)*X
                R(4,J)=Q(4,J)-Q(10,J)*X
              END IF
   23       CONTINUE
            IF(B(1,1).LT.0..OR.B(1,3).LT.0..OR.B(1,6).LT.0..OR.
     *         B(1,10).LT.0..OR.B(1,15).LT.0..OR.B(1,21).LT.0.) THEN
C             <A NAME="error08"></A><FONT COLOR="RED">HDERTEST-08</FONT>
              CALL ERROR('HDERTEST-08: Negative elastic modulus')
C             Negative elastic modulus A1111, A2222, A3333, A2323, A1313
C             or A1212 has been identified during finite differencing with
C             respect to one of spatial coordinates.
C             Check the relation of step DX in spatial coordinates in
C             comparison with the first-order and second-order derivatives
C             of the above mentioned elastic moduli.
            END IF
            CALL HDER(ICB,B,R,P1,P2,P3,HP,HS,F,HI,F1,F2,F3,F4,F5,F6,
     *                F11,F12,F22,F13,F23,F33,F14,F24,F34,F44,
     *                F15,F25,F35,F45,F55,F16,F26,F36,F46,F56,F66,
     *                E,IERR)
            IF(IERR.NE.0) THEN
C             <A NAME="error09"></A><FONT COLOR="RED">HDERTEST-09</FONT>
              CALL ERROR('HDERTEST-09: Error in HDER')
            END IF
            D (IFD/2)=0.25*(       G -        F )/X
            D1(IFD/2)=0.5*(     G *G1-     F *F1)/X
            D2(IFD/2)=0.5*(     G *G2-     F *F2)/X
            D3(IFD/2)=0.5*(     G *G3-     F *F3)/X
            D4(IFD/2)=0.5*(SQRT(G)*G4-SQRT(F)*F4)/X
            D5(IFD/2)=0.5*(SQRT(G)*G5-SQRT(F)*F5)/X
            D6(IFD/2)=0.5*(SQRT(G)*G6-SQRT(F)*F6)/X
   25     CONTINUE
          IF(NFD.EQ.4) THEN
            D (1)=(4.*D (1)-D (2))/3.
            D1(1)=(4.*D1(1)-D1(2))/3.
            D2(1)=(4.*D2(1)-D2(2))/3.
            D3(1)=(4.*D3(1)-D3(2))/3.
            D4(1)=(4.*D4(1)-D4(2))/3.
            D5(1)=(4.*D5(1)-D5(2))/3.
            D6(1)=(4.*D6(1)-D6(2))/3.
          ELSE IF(NFD.EQ.6) THEN
            D (1)=1.5*D (1)-0.6*D (2)+0.1*D (3)
            D1(1)=1.5*D1(1)-0.6*D1(2)+0.1*D1(3)
            D2(1)=1.5*D2(1)-0.6*D2(2)+0.1*D2(3)
            D3(1)=1.5*D3(1)-0.6*D3(2)+0.1*D3(3)
            D4(1)=1.5*D4(1)-0.6*D4(2)+0.1*D4(3)
            D5(1)=1.5*D5(1)-0.6*D5(2)+0.1*D5(3)
            D6(1)=1.5*D6(1)-0.6*D6(2)+0.1*D6(3)
          END IF
          WRITE(LU,'(7F11.6)') D(1),D1(1),D2(1),D3(1),D4(1),D5(1),D6(1)
   27   CONTINUE
   29 CONTINUE
C
C     Derivatives with respect to the slowness vector:
      DO 59 IDER=1,3
C       Writing the header line and the derivatives calculated by HDER
        IF(IDER.EQ.1) THEN
C         P1 derivative
          WRITE(LU,'(7(4X,A,4X))')
     *                      'H4 ','H14','H24','H34','H44','H54','H64'
          WRITE(LU,'(7F11.6)') E4,E14,E24,E34,E44,E45,E46
          WRITE(LU,'(7F11.6)') H4,H14,H24,H34,H44,H45,H46
        ELSE IF(IDER.EQ.2) THEN
C         P2 derivative
          WRITE(LU,'(7(4X,A,4X))')
     *                      'H5 ','H15','H25','H35','H45','H55','H65'
          WRITE(LU,'(7F11.6)') E5,E15,E25,E35,E45,E55,E56
          WRITE(LU,'(7F11.6)') H5,H15,H25,H35,H45,H55,H56
        ELSE
C         P3 derivative
          WRITE(LU,'(7(4X,A,4X))')
     *                      'H6 ','H16','H26','H36','H46','H56','H66'
          WRITE(LU,'(7F11.6)') E6,E16,E26,E36,E46,E56,E66
          WRITE(LU,'(7F11.6)') H6,H16,H26,H36,H46,H56,H66
        END IF
C
C       Loop over finite differences:
        DO 57 K=0,NPOWER
          X=DP*2.**(1.0*FLOAT(K))
          DO 55 IFD=2,NFD,2
            IF(IFD.GT.2) THEN
              X=X+X
            END IF
            Q1=P1
            Q2=P2
            Q3=P3
            IF(IDER.EQ.1) THEN
              Q1=P1+X
            ELSE IF(IDER.EQ.2) THEN
              Q2=P2+X
            ELSE
              Q3=P3+X
            END IF
            CALL HDER(ICB,A,Q,Q1,Q2,Q3,HP,HS,G,HI,G1,G2,G3,G4,G5,G6,
     *                G11,G12,G22,G13,G23,G33,G14,G24,G34,G44,
     *                G15,G25,G35,G45,G55,G16,G26,G36,G46,G56,G66,
     *                E,IERR)
            IF(IERR.NE.0) THEN
C             <A NAME="error10"></A><FONT COLOR="RED">HDERTEST-10</FONT>
              CALL ERROR('HDERTEST-10: Error in HDER')
            END IF
            IF(IDER.EQ.1) THEN
              Q1=P1-X
            ELSE IF(IDER.EQ.2) THEN
              Q2=P2-X
            ELSE
              Q3=P3-X
            END IF
            CALL HDER(ICB,A,Q,Q1,Q2,Q3,HP,HS,F,HI,F1,F2,F3,F4,F5,F6,
     *                F11,F12,F22,F13,F23,F33,F14,F24,F34,F44,
     *                F15,F25,F35,F45,F55,F16,F26,F36,F46,F56,F66,
     *                E,IERR)
            IF(IERR.NE.0) THEN
C             <A NAME="error11"></A><FONT COLOR="RED">HDERTEST-11</FONT>
              CALL ERROR('HDERTEST-11: Error in HDER')
            END IF
            D (IFD/2)=0.25*(       G -        F )/X
            D1(IFD/2)=0.5*(     G *G1-     F *F1)/X
            D2(IFD/2)=0.5*(     G *G2-     F *F2)/X
            D3(IFD/2)=0.5*(     G *G3-     F *F3)/X
            D4(IFD/2)=0.5*(SQRT(G)*G4-SQRT(F)*F4)/X
            D5(IFD/2)=0.5*(SQRT(G)*G5-SQRT(F)*F5)/X
            D6(IFD/2)=0.5*(SQRT(G)*G6-SQRT(F)*F6)/X
   55     CONTINUE
          IF(NFD.EQ.4) THEN
            D (1)=(4.*D (1)-D (2))/3.
            D1(1)=(4.*D1(1)-D1(2))/3.
            D2(1)=(4.*D2(1)-D2(2))/3.
            D3(1)=(4.*D3(1)-D3(2))/3.
            D4(1)=(4.*D4(1)-D4(2))/3.
            D5(1)=(4.*D5(1)-D5(2))/3.
            D6(1)=(4.*D6(1)-D6(2))/3.
          ELSE IF(NFD.EQ.6) THEN
            D (1)=1.5*D (1)-0.6*D (2)+0.1*D (3)
            D1(1)=1.5*D1(1)-0.6*D1(2)+0.1*D1(3)
            D2(1)=1.5*D2(1)-0.6*D2(2)+0.1*D2(3)
            D3(1)=1.5*D3(1)-0.6*D3(2)+0.1*D3(3)
            D4(1)=1.5*D4(1)-0.6*D4(2)+0.1*D4(3)
            D5(1)=1.5*D5(1)-0.6*D5(2)+0.1*D5(3)
            D6(1)=1.5*D6(1)-0.6*D6(2)+0.1*D6(3)
          END IF
          WRITE(LU,'(7F11.6)') D(1),D1(1),D2(1),D3(1),D4(1),D5(1),D6(1)
   57   CONTINUE
   59 CONTINUE
C
      CLOSE(LU)
      WRITE(*,'(A)') '+HDERTEST: Done.                 '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     <A HREF="../forms/error.for" TYPE="text/html">error.for</A>
      INCLUDE 'sep.for'
C     <A HREF="../forms/sep.for" TYPE="text/html">sep.for</A>
      INCLUDE 'length.for'
C     <A HREF="../forms/length.for" TYPE="text/html">length.for</A>
      INCLUDE 'hder.for'
C     <A HREF="../forms/hder.for" TYPE="text/html">hder.for</A>
      INCLUDE 'cubic.for'
C     <A HREF="../forms/cubic.for" TYPE="text/html">cubic.for</A>
      INCLUDE 'eigen.for'
C     <A HREF="../forms/eigen.for" TYPE="text/html">eigen.for</A>
C
C=======================================================================
C                                                                 </PRE>
