C
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 hder.for
C
C Program HDERTEST reads the density-reduced stiffness tensor composed
C of elastic moduli, and its first-order and second-order partial
C derivatives with respect 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 subroutine
C HDER.  The kind of Hamiltonian function is specified by input SEP
C parameters 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: 7.20
C Date: 2015, June 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                                                    
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                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
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             Description of file ANIDER.
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             Description of file HDERTEST.
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     Input SEP parameters of subroutine HDER.
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     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                                                  
C Input file ANIDER with the density-reduced elastic moduli
C and their first-order and second-order spatial derivatives:
C (1) 21 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 elastic moduli with respect to the
C     X1 coordinate axis read at once in the same order as above.
C (3) First-order derivatives of 21 elastic moduli with respect to the
C     X2 coordinate axis read at once in the same order.
C (4) First-order derivatives of 21 elastic moduli with respect to the
C     X3 coordinate axis read at once in the same order.
C (5) Homogeneous second-order derivatives of 21 elastic moduli with
C     respect to the X1 coordinate axis read at once in the same order.
C (6) Mixed second-order derivatives of 21 elastic moduli with respect
C     to the X1 and X2 coordinate axes read at once in the same order.
C (7) Homogeneous second-order derivatives of 21 elastic moduli with
C     respect to the X2 coordinate axis read at once in the same order.
C (8) Mixed second-order derivatives of 21 elastic moduli with respect
C     to the X1 and X3 coordinate axes read at once in the same order.
C (9) Mixed second-order derivatives of 21 elastic moduli with respect
C     to the X2 and X3 coordinate axes read at once in the same order.
C (10) Homogeneous second-order derivatives of 21 elastic moduli with
C     respect to the X3 coordinate axis read at once in the same order.
C
C                                                
C Output file HDERTEST contains 6 header lines specifying the
C phase-space derivatives of the Hamiltonian function.
C The header lines correspond to symmetric finite differencing with
C respect to the X1, X2, X3, X4=P1, X5=P2 and X6=P3 phase-space
C coordinates.  The line below each header line contains the values
C calculated by subroutine HDER.  The subsequent NPOWER+1 lines contain
C the values estimated by the symmetric finite differences.
C
C=======================================================================
C
      EXTERNAL ERROR,RSEP1,RSEP3T,EIGEN
C
C     Filenames and parameters:
      CHARACTER*80 FILE
      INTEGER LU
      PARAMETER (LU=1)
C
C     Stiffness tensor with derivatives
      REAL A(10,21),B(10,21)
C
C     Arguments of HDER:
      INTEGER ICB,IERR
      REAL P1,P2,P3,HP,HS,H,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     Other variables:
      INTEGER I,J,K,NPOWER
      REAL DX,DP
      REAL X,D,D1,D2,D3,D4,D5,D6
C
C-----------------------------------------------------------------------
C
C     Reading main input data:
      WRITE(*,'(A)') '+HDERTEST: Enter input filename: '
      FILE=' '
      READ(*,*) FILE
      IF(FILE.EQ.' ') THEN
C       HDERTEST-01
        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       HDERTEST-02
        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 10 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)
   10 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     Reading numerical parameters of the test:
      CALL RSEP3I('ICB',ICB,-1)
      CALL RSEP3R('DX',DX,0.002)
      CALL RSEP3R('DP',DP,0.001)
      CALL RSEP3I('NPOWER',NPOWER,0)
C
C     Calculating phase-space derivatives of the Hamiltonian function:
      CALL HDER(ICB,A,P1,P2,P3,HP,HS,H,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       HDERTEST-03
        CALL ERROR('HDERTEST-03: Error in HDER')
      END IF
C     Normalizing the given slowness vector:
      P1=P1/SQRT(H)
      P2=P2/SQRT(H)
      P3=P3/SQRT(H)
C
C     Opening the output file:
      CALL RSEP3T('HDERTEST',FILE,'hdertest.out')
      OPEN(LU,FILE=FILE)
C
C     Spatial derivatives:
      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)') 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)') 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)') H3,H13,H23,H33,H34,H35,H36
        END IF
C
C       Loop over finite differences:
        DO 24 K=0,NPOWER
          X=DX*2.**(1.0*FLOAT(K))
          DO 22 J=1,21
            DO 21 I=1,10
              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
            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
            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
            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
            GO TO 28
          END IF
          CALL HDER(ICB,B,P1,P2,P3,HP,HS,G,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           HDERTEST-04
            CALL ERROR('HDERTEST-04: 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
            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
            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
            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
            GO TO 28
          END IF
          CALL HDER(ICB,B,P1,P2,P3,HP,HS,F,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           HDERTEST-05
            CALL ERROR('HDERTEST-05: Error in HDER')
          END IF
          D =0.25*(       G -        F )/X
          D1=0.5*(     G *G1-     F *F1)/X
          D2=0.5*(     G *G2-     F *F2)/X
          D3=0.5*(     G *G3-     F *F3)/X
          D4=0.5*(SQRT(G)*G4-SQRT(F)*F4)/X
          D5=0.5*(SQRT(G)*G5-SQRT(F)*F5)/X
          D6=0.5*(SQRT(G)*G6-SQRT(F)*F6)/X
          WRITE(LU,'(7F11.6)') D,D1,D2,D3,D4,D5,D6
   24   CONTINUE
   28   CONTINUE
   29 CONTINUE
C
C     Derivatives withrespect 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)') 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)') 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)') H6,H16,H26,H36,H46,H56,H66
        END IF
C
C       Loop over finite differences:
        DO 54 K=0,NPOWER
          X=DP*2.**(1.0*FLOAT(K))
          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,Q1,Q2,Q3,HP,HS,G,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(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,Q1,Q2,Q3,HP,HS,F,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)
          D =0.25*(       G -        F )/X
          D1=0.5*(     G *G1-     F *F1)/X
          D2=0.5*(     G *G2-     F *F2)/X
          D3=0.5*(     G *G3-     F *F3)/X
          D4=0.5*(SQRT(G)*G4-SQRT(F)*F4)/X
          D5=0.5*(SQRT(G)*G5-SQRT(F)*F5)/X
          D6=0.5*(SQRT(G)*G6-SQRT(F)*F6)/X
          WRITE(LU,'(7F11.6)') D,D1,D2,D3,D4,D5,D6
   54   CONTINUE
   59 CONTINUE
C
      CLOSE(LU)
      WRITE(*,'(A)') '+HDERTEST: Done.                 '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'eigen.for'
C     eigen.for
      INCLUDE 'hder.for'
C     hder.for
C
C=======================================================================
C