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