C
C Program TIAXIS to determine the reference symmetry axis of a generally C anisotropic medium which is approximately transversely isotropic C C Version: 7.20 C Date: 2015, May 4 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 ANIMOD='string'... Name of the input file containing the C density-reduced elastic moduli. C Description of file ANIMOD. C No default. ANIMOD must be specified and cannot be blank. C Data specifying output file: C TIAXIS='string'... Name of the output file containing C the components of the reference symmetry vector, C and the ratio of the norm of the derivative of the C stiffness tensor with respect to the angle of rotation C and the norm of the stiffness tensor. C Default: TIAXIS='tiaxis.out' C C C Input file ANIMOD with the density-reduced elastic moduli: 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 C======================================================================= C EXTERNAL ERROR,RSEP1,RSEP3T,EIGEN C C Filenames and parameters: CHARACTER*80 FILE INTEGER LU PARAMETER (LU=1) C C Stiffness tensor REAL A11,A12,A13,A14,A15,A16 REAL A21,A22,A23,A24,A25,A26 REAL A31,A32,A33,A34,A35,A36 REAL A41,A42,A43,A44,A45,A46 REAL A51,A52,A53,A54,A55,A56 REAL A61,A62,A63,A64,A65,A66 REAL A(3,3,3,3) C C Other variables: INTEGER I,J,K,L,M,N,MN REAL EPS(3,3,3),D(3,3,3,3,3),B(6),E(9),ANORM,RHO1,RHO2,RHO3 C C Debugging: * INTEGER IJ,KL C C----------------------------------------------------------------------- C C Reading main input data: WRITE(*,'(A)') '+TIAXIS: Enter input filename: ' FILE=' ' READ(*,*) FILE IF(FILE.EQ.' ') THEN C TIAXIS-01 CALL ERROR('TIAXIS-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)') '+TIAXIS: Working... ' C C Reading stiffness tensor or density reduced elastic moduli: CALL RSEP3T('ANIMOD',FILE,' ') IF(FILE.EQ.' ') THEN C TIAXIS-02 CALL ERROR('TIAXIS-02: No input stiffness tensor specified') C Name ANIMOD of the input file with the stiffness tensor C must be specified and cannot be blank. END IF OPEN(LU,FILE=FILE,STATUS='OLD') READ(LU,*) A11,A12,A13,A14,A15,A16 * ,A22,A23,A24,A25,A26 * ,A33,A34,A35,A36 * ,A44,A45,A46 * ,A55,A56 * ,A66 CLOSE(LU) C C Completing the stiffness tensor A(I,J,K,L) of elastic moduli A21=A12 A31=A13 A41=A14 A51=A15 A61=A16 A32=A23 A42=A24 A52=A25 A62=A26 A43=A34 A53=A35 A63=A36 A54=A45 A64=A46 A65=A56 A(1,1,1,1)=A11 A(2,2,1,1)=A21 A(3,3,1,1)=A31 A(2,3,1,1)=A41 A(3,2,1,1)=A41 A(1,3,1,1)=A51 A(3,1,1,1)=A51 A(1,2,1,1)=A61 A(2,1,1,1)=A61 A(1,1,2,2)=A12 A(2,2,2,2)=A22 A(3,3,2,2)=A32 A(2,3,2,2)=A42 A(3,2,2,2)=A42 A(1,3,2,2)=A52 A(3,1,2,2)=A52 A(1,2,2,2)=A62 A(2,1,2,2)=A62 A(1,1,3,3)=A13 A(2,2,3,3)=A23 A(3,3,3,3)=A33 A(2,3,3,3)=A43 A(3,2,3,3)=A43 A(1,3,3,3)=A53 A(3,1,3,3)=A53 A(1,2,3,3)=A63 A(2,1,3,3)=A63 A(1,1,2,3)=A14 A(2,2,2,3)=A24 A(3,3,2,3)=A34 A(2,3,2,3)=A44 A(3,2,2,3)=A44 A(1,3,2,3)=A54 A(3,1,2,3)=A54 A(1,2,2,3)=A64 A(2,1,2,3)=A64 A(1,1,3,2)=A14 A(2,2,3,2)=A24 A(3,3,3,2)=A34 A(2,3,3,2)=A44 A(3,2,3,2)=A44 A(1,3,3,2)=A54 A(3,1,3,2)=A54 A(1,2,3,2)=A64 A(2,1,3,2)=A64 A(1,1,1,3)=A15 A(2,2,1,3)=A25 A(3,3,1,3)=A35 A(2,3,1,3)=A45 A(3,2,1,3)=A45 A(1,3,1,3)=A55 A(3,1,1,3)=A55 A(1,2,1,3)=A65 A(2,1,1,3)=A65 A(1,1,3,1)=A15 A(2,2,3,1)=A25 A(3,3,3,1)=A35 A(2,3,3,1)=A45 A(3,2,3,1)=A45 A(1,3,3,1)=A55 A(3,1,3,1)=A55 A(1,2,3,1)=A65 A(2,1,3,1)=A65 A(1,1,1,2)=A16 A(2,2,1,2)=A26 A(3,3,1,2)=A36 A(2,3,1,2)=A46 A(3,2,1,2)=A46 A(1,3,1,2)=A56 A(3,1,1,2)=A56 A(1,2,1,2)=A66 A(2,1,1,2)=A66 A(1,1,2,1)=A16 A(2,2,2,1)=A26 A(3,3,2,1)=A36 A(2,3,2,1)=A46 A(3,2,2,1)=A46 A(1,3,2,1)=A56 A(3,1,2,1)=A56 A(1,2,2,1)=A66 A(2,1,2,1)=A66 C C Levi-Civita symbols EPS(I,J,K) DO 23 K=1,3 DO 22 J=1,3 DO 21 I=1,3 EPS(I,J,K)=0. 21 CONTINUE 22 CONTINUE 23 CONTINUE EPS(1,2,3)=1. EPS(3,1,2)=1. EPS(2,3,1)=1. EPS(3,2,1)=-1. EPS(1,3,2)=-1. EPS(2,1,3)=-1. C C Derivative D(I,J,K,L,M) of elastic moduli with respect to rotation ANORM=0. DO 36 I=1,3 DO 35 J=1,3 DO 34 K=1,3 DO 33 L=1,3 ANORM=ANORM+A(I,J,K,L)*A(I,J,K,L) DO 32 M=1,3 D(I,J,K,L,M)=0. DO 31 N=1,3 D(I,J,K,L,M)=D(I,J,K,L,M)+EPS(M,I,N)*A(N,J,K,L) * +EPS(M,J,N)*A(I,N,K,L) * +EPS(M,K,N)*A(I,J,N,L) * +EPS(M,L,N)*A(I,J,K,N) 31 CONTINUE 32 CONTINUE 33 CONTINUE 34 CONTINUE 35 CONTINUE 36 CONTINUE C C Quadratic B(MN) form of derivatives D(I,J,K,L,M) MN=0 DO 46 N=1,3 DO 45 M=1,N MN=MN+1 B(MN)=0. DO 44 L=1,3 DO 43 K=1,3 DO 42 J=1,3 DO 41 I=1,3 B(MN)=B(MN)+D(I,J,K,L,M)*D(I,J,K,L,N) 41 CONTINUE 42 CONTINUE 43 CONTINUE 44 CONTINUE 45 CONTINUE 46 CONTINUE C C Eigenvalues and eigenvectors of quadratic form B(MN) CALL EIGEN(B,E,3,0) RHO1=SQRT(B(1)/ANORM) RHO2=SQRT(B(3)/ANORM) RHO3=SQRT(B(6)/ANORM) C C Output file with the reference symmetry vector: CALL RSEP3T('TIAXIS',FILE,'tiaxis.out') OPEN(LU,FILE=FILE) WRITE(LU,'(6(A,I1,A,F8.6))') 'TI',1,'=',E(7), * ' TI',2,'=',E(8), * ' TI',3,'=',E(9), * ' RHO',3,'=',RHO3, * ' RHO',2,'=',RHO2, * ' RHO',1,'=',RHO1 C C Derivative of the stiffness tensor for debugging * DO 52 KL=1,6 * IF(KL.LE.3) THEN * K=KL * L=KL * ELSE * K=(9-KL-1)/2 * L=(9-KL+2)/2 * END IF * DO 51 IJ=1,KL * IF(IJ.LE.3) THEN * I=IJ * J=IJ * ELSE * I=(9-IJ-1)/2 * J=(9-IJ+2)/2 * END IF * WRITE(LU,'(A,4I1,A,2E14.6)') ' a''',I,J,K,L,'=', * * D(I,J,K,L,1)*E(7)+D(I,J,K,L,2)*E(8)+D(I,J,K,L,3)*E(9) * 51 CONTINUE * 52 CONTINUE C CLOSE(LU) WRITE(*,'(A)') '+TIAXIS: 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 C C======================================================================= C