C
C Program ANIROT to rotate the tensor of elastic moduli C of an anisotropic elastic medium C C Version: 7.30 C Date: 2015, December 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 Names of the input and output files: C ANIIN='string'... Name of the input file containing the elastic C moduli before rotation. C Description of file ANIIN. C Default: ANIIN='anirot.dat' C ANIOUT='string'... Name of the output file containing the elastic C moduli after rotation. C Description of file ANIOUT. C Default: ANIOUT='anirot.out' C Data describing the rotation: C ANIDEG=real... Rotation angle in degrees (not in radians). C ANIE1=real, ANIE2=real, ANIE3=real... Three components of the C axial vector around which the anisotropic elastic medium C is rotated. The length of the axial vector does not C matter. The medium is rotated according to the right-hand C rule in a right-handed coordinate system. The medium is C thus rotated counter-clockwise when looking against the C axial vector. C Transformation matrix E_ij: C The i-th component of the j-th rotated local Cartesian C basis vector, expressed with respect to the local C Cartesian basis vectors before rotation, is C E_ij = cos(ANIDEG) delta_ij C + [1-cos(ANIDEG)] E_i/E E_j/E C - sin(ANIDEG) epsilon_ijk E_k/E C where (E_1,E_2,E_3)=(ANIE1,ANIE2,ANIE3), C E = sqrt(E_k E_k) C is the Cartesian length of the rotation vector, C delta_ij is the Kronecker delta, and C epsilon_ijk is completely skew Levi-Civita's symbol, C with epsilon_123=epsilon_231=epsilon_312=1. C Transformation of the elastic moduli: C B_ijkl = E_ia E_ib E_ic E_id A_abcd C where the Einstein summation is applied, C A_ijkl are the elastic moduli before rotation, and C B_ijkl are the elastic moduli after rotation. C Default: C ANIE1=0., ANIE2=0., ANIE3=0. (no rotation) C Optional parameters specifying the form of the quantities C written in the output formatted file: C MINDIG,MAXDIG=positive integers ... See the description in file C forms.for. C Recommended value of MAXDIG is 7 or 8 for this program. C C C Input file ANIIN with the elastic moduli before rotation: 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 Output file ANIOUT with the elastic moduli after rotation: C (1) Elastic moduli C B1111,B1122,B1133,B1123,B1113,B1112 C (2) Elastic moduli C B2222,B2233,B2223,B2213,B2212 C (3) Elastic moduli C B3333,B3323,B3313,B3312 C (4) Elastic moduli C B2323,B2313,B2312 C (5) Elastic moduli C B1313,B1312 C (6) Elastic modulus C B1212 C C======================================================================= C C Filenames and parameters: CHARACTER*80 FSEP,FIN,FOUT INTEGER LU1 PARAMETER (LU1=1) C C Other variables: REAL DEG,E1,E2,E3 REAL RAD,C,S,BMIN,BMAX,AUX1,AUX2,AUX3,AUX4 REAL AV(6,6),A(3,3,3,3),B(3,3,3,3),E(3,3) INTEGER I12,I34,I1,I2,I3,I4,J1,J2,J3,J4 CHARACTER*17 FORMAT C C....................................................................... C C Reading main input data: WRITE(*,'(A)') '+ANIROT: Enter input filename: ' FSEP=' ' READ (*,*) FSEP IF(FSEP.EQ.' ') THEN C ANIROT-01 CALL ERROR('ANIROT-01: No input 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(LU1,FSEP) WRITE(*,'(A)') '+ANIROT: Working... ' C C Reading input SEP parameters: CALL RSEP3T('ANIIN' ,FIN ,'anirot.dat') CALL RSEP3T('ANIOUT',FOUT,'anirot.out') CALL RSEP3R('ANIDEG',DEG,0.) CALL RSEP3R('ANIE1' ,E1 ,0.) CALL RSEP3R('ANIE2' ,E2 ,0.) CALL RSEP3R('ANIE3' ,E3 ,0.) C C Transformation matrix E_ij: RAD=3.141592654*DEG/180. AUX1=SQRT(E1**2+E2**2+E3**2) IF(RAD.EQ.0..OR.AUX1.EQ.0.) THEN E(1,1)=1. E(2,1)=0. E(3,1)=0. E(1,2)=0. E(2,2)=1. E(3,2)=0. E(1,3)=0. E(2,3)=0. E(3,3)=1. ELSE E1=E1/AUX1 E2=E2/AUX1 E3=E3/AUX1 C=COS(RAD) S=SIN(RAD) AUX1=1.-C E(1,1)=AUX1*E1*E1+C E(2,1)=AUX1*E2*E1+S*E3 E(3,1)=AUX1*E3*E1-S*E2 E(1,2)=AUX1*E1*E2-S*E3 E(2,2)=AUX1*E2*E2+C E(3,2)=AUX1*E3*E2+S*E1 E(1,3)=AUX1*E1*E3+S*E2 E(2,3)=AUX1*E2*E3-S*E1 E(3,3)=AUX1*E3*E3+C END IF C C Reading the elastic moduli in Voight notation: OPEN(LU1,FILE=FIN,STATUS='OLD') READ(LU1,*) AV(1,1),AV(1,2),AV(1,3),AV(1,4),AV(1,5),AV(1,6) * ,AV(2,2),AV(2,3),AV(2,4),AV(2,5),AV(2,6) * ,AV(3,3),AV(3,4),AV(3,5),AV(3,6) * ,AV(4,4),AV(4,5),AV(4,6) * ,AV(5,5),AV(5,6) * ,AV(6,6) CLOSE(LU1) C C Transforming the Voight notation to the tensor: DO 14 I4=1,3 DO 13 I3=1,3 IF(I3.EQ.I4) THEN I34=I3 ELSE I34=9-I3-I4 END IF DO 12 I2=1,3 DO 11 I1=1,3 IF(I1.EQ.I2) THEN I12=I1 ELSE I12=9-I1-I2 END IF IF(I12.LE.I34) THEN A(I1,I2,I3,I4)=AV(I12,I34) ELSE A(I1,I2,I3,I4)=AV(I34,I12) END IF 11 CONTINUE 12 CONTINUE 13 CONTINUE 14 CONTINUE C C Transformation of the elastic moduli: BMIN= 1.0E20 BMAX=-1.0E20 DO 34 J4=1,3 DO 33 J3=1,3 DO 32 J2=1,3 DO 31 J1=1,3 AUX1=0. DO 24 I4=1,3 AUX4=E(J4,I4) DO 23 I3=1,3 AUX3=AUX4*E(J3,I3) DO 22 I2=1,3 AUX2=AUX3*E(J2,I2) DO 21 I1=1,3 AUX1=AUX1+AUX2*E(J1,I1)*A(I1,I2,I3,I4) 21 CONTINUE 22 CONTINUE 23 CONTINUE 24 CONTINUE B(J1,J2,J3,J4)=AUX1 BMIN=AMIN1(BMIN,AUX1) BMAX=AMAX1(BMAX,AUX1) 31 CONTINUE 32 CONTINUE 33 CONTINUE 34 CONTINUE C C Setting the output format: FORMAT='( 6(F00.0,1X))' CALL FORM1(BMIN,BMAX,FORMAT(8:15)) FORMAT(14:17)= '1X))' I1=ICHAR(FORMAT(9:9))-ICHAR('0') I1=10*I1+ICHAR(FORMAT(10:10))-ICHAR('0')+1 C C Writing the elastic moduli in Voight notation: OPEN(LU1,FILE=FOUT) WRITE(LU1,FORMAT) * B(1,1,1,1),B(1,1,2,2),B(1,1,3,3),B(1,1,2,3),B(1,1,1,3),B(1,1,1,2) I2=I1 I3=I2/10 I4=I2-I3*10 FORMAT(2:2)=CHAR(ICHAR('0')+I3) FORMAT(3:3)=CHAR(ICHAR('0')+I4) FORMAT(4:5)='X,' WRITE(LU1,FORMAT) * ,B(2,2,2,2),B(2,2,3,3),B(2,2,2,3),B(2,2,1,3),B(2,2,1,2) I2=I2+I1 I3=I2/10 I4=I2-I3*10 FORMAT(2:2)=CHAR(ICHAR('0')+I3) FORMAT(3:3)=CHAR(ICHAR('0')+I4) WRITE(LU1,FORMAT) * ,B(3,3,3,3),B(3,3,2,3),B(3,3,1,3),B(3,3,1,2) I2=I2+I1 I3=I2/10 I4=I2-I3*10 FORMAT(2:2)=CHAR(ICHAR('0')+I3) FORMAT(3:3)=CHAR(ICHAR('0')+I4) WRITE(LU1,FORMAT) * ,B(2,3,2,3),B(2,3,1,3),B(2,3,1,2) I2=I2+I1 I3=I2/10 I4=I2-I3*10 FORMAT(2:2)=CHAR(ICHAR('0')+I3) FORMAT(3:3)=CHAR(ICHAR('0')+I4) WRITE(LU1,FORMAT) * ,B(1,3,1,3),B(1,3,1,2) I2=I2+I1 I3=I2/10 I4=I2-I3*10 FORMAT(2:2)=CHAR(ICHAR('0')+I3) FORMAT(3:3)=CHAR(ICHAR('0')+I4) WRITE(LU1,FORMAT) * ,B(1,2,1,2) CLOSE(LU1) WRITE(*,'(A)') '+ANIROT: 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 'forms.for' C forms.for C C======================================================================= C