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.30 C Date: 2016, March 25 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 files: 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 ANITI='string'... Optional name of the output file containing the C density-reduced elastic moduli of the reference C transversely isotropic medium. C This output file has the same form as input file ANIMOD. C Default: ANITI=' ' (no output) 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 Variables related to optional reference TI stiffness tensor: REAL T(3),Z(3,3),C(3,3),S(3,3),ATI(3,3,3,3) 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 C----------------------------------------------------------------------- C C Optional reference TI stiffness tensor: C C Reference symmetry vector T(1)=E(7) T(2)=E(8) T(3)=E(9) C C Generating transformation matrices Z, C and S DO 62 J=1,3 DO 61 I=1,3 Z(I,J)=T(I)*T(J) IF(I.EQ.J) THEN C(I,J)=1.-Z(I,J) ELSE C(I,J)= -Z(I,J) K=6-I-J IF(I.LT.J) THEN S(I,J)=T(K) ELSE S(I,J)=-T(K) END IF IF(K.EQ.2) THEN S(I,J)=-S(I,J) END IF END IF 61 CONTINUE 62 CONTINUE C C Initializing transversely isotropic stiffness tensor ATI DO 74 L=1,3 DO 73 K=1,3 DO 72 J=1,3 DO 71 I=1,3 ATI(I,J,K,L)=0. 71 CONTINUE 72 CONTINUE 73 CONTINUE 74 CONTINUE C C Creating transversely isotropic stiffness tensor ATI CALL TRANSA(1.000,Z,Z,Z,Z,A,ATI) C CALL TRANSA(0.500,Z,Z,C,C,A,ATI) CALL TRANSA(0.500,Z,C,Z,C,A,ATI) CALL TRANSA(0.500,Z,C,C,Z,A,ATI) CALL TRANSA(0.500,C,Z,Z,C,A,ATI) CALL TRANSA(0.500,C,Z,C,Z,A,ATI) CALL TRANSA(0.500,C,C,Z,Z,A,ATI) C CALL TRANSA(0.500,Z,Z,S,S,A,ATI) CALL TRANSA(0.500,Z,S,Z,S,A,ATI) CALL TRANSA(0.500,Z,S,S,Z,A,ATI) CALL TRANSA(0.500,S,Z,Z,S,A,ATI) CALL TRANSA(0.500,S,Z,S,Z,A,ATI) CALL TRANSA(0.500,S,S,Z,Z,A,ATI) C CALL TRANSA(0.125,C,C,S,S,A,ATI) CALL TRANSA(0.125,C,S,C,S,A,ATI) CALL TRANSA(0.125,C,S,S,C,A,ATI) CALL TRANSA(0.125,S,C,C,S,A,ATI) CALL TRANSA(0.125,S,C,S,C,A,ATI) CALL TRANSA(0.125,S,S,C,C,A,ATI) C CALL TRANSA(0.375,C,C,C,C,A,ATI) CALL TRANSA(0.375,S,S,S,S,A,ATI) C C Relative Frobenius norm of the difference RNORM=0. DO 84 L=1,3 DO 83 K=1,3 DO 82 J=1,3 DO 81 I=1,3 RNORM=RNORM+(ATI(I,J,K,L)-A(I,J,K,L))**2 81 CONTINUE 82 CONTINUE 83 CONTINUE 84 CONTINUE RNORM=SQRT(RNORM/ANORM) WRITE(LU,'(A,F8.6)') 'Relative difference=',RNORM CLOSE(LU) C C Writing transversely isotropic stiffness tensor ATI CALL RSEP3T('ANITI',FILE,' ') IF(FILE.EQ.' ') THEN GO TO 90 END IF OPEN(LU,FILE=FILE) WRITE(LU,'( 6(F9.6,1X))')ATI(1,1,1,1),ATI(1,1,2,2),ATI(1,1,3,3) * ,ATI(1,1,2,3),ATI(1,1,1,3),ATI(1,1,1,2) WRITE(LU,'(10X,5(F9.6,1X))') ATI(2,2,2,2),ATI(2,2,3,3) * ,ATI(2,2,2,3),ATI(2,2,1,3),ATI(2,2,1,2) WRITE(LU,'(20X,4(F9.6,1X))') ATI(3,3,3,3) * ,ATI(3,3,2,3),ATI(3,3,1,3),ATI(3,3,1,2) WRITE(LU,'(30X,3(F9.6,1X))')ATI(2,3,2,3),ATI(2,3,1,3),ATI(2,3,1,2) WRITE(LU,'(40X,2(F9.6,1X))') ATI(1,3,1,3),ATI(1,3,1,2) WRITE(LU,'(50X,1(F9.6,1X))') ATI(1,2,1,2) CLOSE(LU) C 90 CONTINUE WRITE(*,'(A)') '+TIAXIS: Done. ' STOP END C C======================================================================= C SUBROUTINE TRANSA(FACT,B1,B2,B3,B4,A1,A2) REAL FACT,B1(3,3),B2(3,3),B3(3,3),B4(3,3),A1(3,3,3,3),A2(3,3,3,3) C C----------------------------------------------------------------------- C INTEGER I1,I2,I3,I4,J1,J2,J3,J4 REAL B34,B234,AUX C DO 80 J4=1,3 DO 70 J3=1,3 DO 60 J2=1,3 DO 50 J1=1,3 AUX=0. DO 40 I4=1,3 DO 30 I3=1,3 B34=B3(J3,I3)*B4(J4,I4) DO 20 I2=1,3 B234=B2(J2,I2)*B34 DO 10 I1=1,3 AUX=AUX+B1(J1,I1)*B234*A1(I1,I2,I3,I4) 10 CONTINUE 20 CONTINUE 30 CONTINUE 40 CONTINUE A2(J1,J2,J3,J4)=A2(J1,J2,J3,J4)+FACT*AUX 50 CONTINUE 60 CONTINUE 70 CONTINUE 80 CONTINUE C RETURN 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