C
C Program ANISRF to triangularize the phase-slowness surface C or the ray-velocity surface C C Version: 7.40 C Date: 2016, November 8 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 If ANIMOD=' ', a unit sphere is triangularized. C Description of file ANIMOD. C Default: ANIMOD=' ' C Data specifying output files: C VRTX='string'... Name of the file with vertices of the triangles. C Description of file VRTX C Default: VRTX='vrtx.out' C TRGL='string'... Name of the file describing the triangles. C Description of file TRGL C Default: TRGL='trgl.out' C ANISTRE='string'... Name of the optional output file containing C the minimum absolute value of the slowness vector Pmin in C the medium, the maximum absolute value of the slowness C vector Pmax, and the anistropy strength defined as C 2*(Pmax-Pmin)/(Pmax+Pmin). C Default: ANISTRE=' ' (no file generated) C Data specifying the triangularized surface: C KSURF=integer... Surface to be triangularized. C KSURF=0: Phase-slowness surface. C Else: Ray-velocity surface. C Default: KSURF=0 C KWAVE=integer... Selection of elastic wave. C KWAVE=1: Slower S wave. C KWAVE=2: Faster S wave. C Else: P wave. C Default: KWAVE=1 C NDIVTRGL=non-negative integer... Specification of the density of C triangle vertices along the unit sphere. C The number of vertices is 10*4**NDIVTRGL+2, C the number of triangles is 20*4**NDIVTRGL. C For NDIVTRGL=0, twelve vertices and twenty triangles C corresponding to a regular icosahedron are generated. C For NDIVTRGL=1, new vertices at the centres of the edges C of the icosahedron are created, and each triangle of the C icosahedron is divided into 4 new triangles. C For NDIVTRGL=2, new vertices at the centres of the old C edges are created, and each old triangle is divided into C 4 new triangles, and so on. C Default: NDIVTRGL=4 C Optional parameters specifying the form of the quantities C written in the output formatted file VRTX: C MINDIG,MAXDIG=positive integers ... See the description in file C forms.for. 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 Output file VRTX with the vertices: C (1) / (a slash) C (2) For each vertex data (2.1): C (2.1) 'NAME',X1,X2,X3,Z1,Z2,Z3,KURV,/ C 'NAME'... Name of the vertex. C X1,X2,X3... Coordinates of the vertex. C Z1,Z2,Z3... Normal to the surface at the vertex. C KURV... Determines the signature of the curvature matrix. C KURV=0: Convex (signature 1,1) C KURV=1: Hyperbolic (signature 1,-1) C KURV=2: Concave (signature -1,-1) C KURV=-1: Infinite or undefined curvature. C /... Slash. C (3) / (a slash) C C C Output file TRGL with the triangles: C (1) For each triangle data (1.1): C (1.1) I1,I2,I3,/ C I1,I2,I3... Indices of 3 vertices of the triangle, right-handed C with respect to the given surface normals. C The vertices in file VRTX are indexed by positive integers C according to their order. C /... List of vertices is terminated by a slash. C C======================================================================= C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C C....................................................................... C EXTERNAL ERROR,RSEP1,RSEP3I,ANISRF C C Filenames and parameters: CHARACTER*80 FILE INTEGER LU PARAMETER (LU=1) C C Other variables: INTEGER NDIVTR,NP,NT C C----------------------------------------------------------------------- C C Reading main input data: WRITE(*,'(A)') '+ANISRF: Enter input filename: ' FILE=' ' READ(*,*) FILE IF(FILE.EQ.' ') THEN C ANISRF-01 CALL ERROR('ANISRF-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(LU,FILE) WRITE(*,'(A)') '+ANISRF: Working... ' C C Allocating memory and triangularizing the surface: CALL RSEP3I('NDIVTRGL',NDIVTR,4) NP=10*4**NDIVTR+2 NT=20*4**NDIVTR IF(210+7*NP+NT.GT.MRAM) THEN C ANISRF-02 CALL ERROR('ANISRF-02: Too small array RAM') C Dimension MRAM of array RAM in include file 'ram.inc' must be at C least 90*4**NDIVTRGL+212, where NDIVTRGL is the input SEP C parameter of this program 'anisrf.for'. END IF CALL ANISRF(LU,FILE,NDIVTR,RAM(1),RAM(211),RAM(211+3*NP), * IRAM(211+6*NP),IRAM(211+7*NP)) C WRITE(*,'(A)') '+ANISRF: Done. ' STOP END C C======================================================================= C SUBROUTINE ANISRF(LU,FILE,NDIVTR,A,XP,ZP,KURV,KT) INTEGER LU,NDIVTR CHARACTER*(*) FILE REAL A(10,21),XP(3,10*4**NDIVTR+2),ZP(3,10*4**NDIVTR+2) INTEGER KURV(10*4**NDIVTR+2),KT(3,20*4**NDIVTR) C C----------------------------------------------------------------------- C EXTERNAL ERROR,RSEP3T,RSEP3I,FORM1,GAA,GABP,GAAPP,ICOSAH C C....................................................................... C C Output format: CHARACTER*36 FORMAT C C Other variables: INTEGER NP,NT,KWAVE,KSURF,I,J REAL P1,P2,P3,V1,V2,V3,Z1,Z2,Z3 REAL G,G1,G2,G3,GA1,GA2,GA3,GB1,GB2,GB3,G11,G12,G22,G13,G23,G33 REAL E11,E21,E31,E12,E22,E32,E13,E23,E33,WA,WB,XPMIN,XPMAX,DG,AUX REAL ABSP,ABSPMI,ABSPMA LOGICAL LANISO C C----------------------------------------------------------------------- C C Number of vertices, number of triangles: NP=10*4**NDIVTR+2 NT=20*4**NDIVTR C C Triangularizing a unit sphere: CALL ICOSAH(NDIVTR,KT,XP) C C Reading density reduced elastic moduli: CALL RSEP3T('ANIMOD',FILE,' ') IF(FILE.NE.' ') THEN OPEN(LU,FILE=FILE,STATUS='OLD') READ(LU,*) A(1,1),A(1,2),A(1,4),A(1, 7),A(1,11),A(1,16) * ,A(1,3),A(1,5),A(1, 8),A(1,12),A(1,17) * ,A(1,6),A(1, 9),A(1,13),A(1,18) * ,A(1,10),A(1,14),A(1,19) * ,A(1,15),A(1,20) * ,A(1,21) CLOSE(LU) DO 12 J=1,21 DO 11 I=2,10 A(I,J)=0.0 11 CONTINUE 12 CONTINUE C C Reading name of the output file with anisotropy strength CALL RSEP3T('ANISTRE',FILE,' ') IF(FILE.NE.' ') THEN LANISO=.TRUE. ELSE LANISO=.FALSE. ENDIF C C Calculating slowness vectors, ray-velocity vectors C and wave-propagation metric tensors: CALL RSEP3I('KWAVE',KWAVE,1) CALL RSEP3I('KSURF',KSURF,0) DO 13 I=1,NP P1=XP(1,I) P2=XP(2,I) P3=XP(3,I) C Eigenvalues of the Christoffel matrix CALL GAA(A,P1,P2,P3, * G1,G2,G3,E11,E21,E31,E12,E22,E32,E13,E23,E33) IF(KWAVE.EQ.1) THEN C S1 wave (the slowest wave) G=G1 IF(G.LE.0.0) THEN C ANISRF-03 CALL ERROR('ANISRF-03: Zero eigenvalue of Christ.matrix') C The eigenvalue of the Christoffel matrix corresponding C to the selected wave IWAVE is zero. This may happen, C e.g., if an S wave is selected in a liquid. END IF WA=G-G2 WB=G-G3 CALL GABP(G,E12,E22,E32,E11,E21,E31,GA1,GA2,GA3) CALL GABP(G,E13,E23,E33,E11,E21,E31,GB1,GB2,GB3) CALL GAAPP(E11,E21,E31,G11,G12,G22,G13,G23,G33) ELSE IF(KWAVE.EQ.2) THEN C S2 wave G=G2 IF(G.LE.0.0) THEN C ANISRF-04 CALL ERROR('ANISRF-04: Zero eigenvalue of Christ.matrix') C The eigenvalue of the Christoffel matrix corresponding C to the selected wave IWAVE is zero. This may happen, C e.g., if an S wave is selected in a liquid. END IF WA=G-G1 WB=G-G3 CALL GABP(G,E11,E21,E31,E12,E22,E32,GA1,GA2,GA3) CALL GABP(G,E13,E23,E33,E12,E22,E32,GB1,GB2,GB3) CALL GAAPP(E12,E22,E32,G11,G12,G22,G13,G23,G33) ELSE C P wave G=G3 IF(G.LE.0.0) THEN C ANISRF-05 CALL ERROR('ANISRF-05: Zero eigenvalue of Christ.matrix') C The eigenvalue of the Christoffel matrix corresponding C to the selected wave IWAVE is zero. This may happen, C e.g., if an S wave is selected in a liquid. END IF WA=G-G1 WB=G-G2 CALL GABP(G,E11,E21,E31,E13,E23,E33,GA1,GA2,GA3) CALL GABP(G,E12,E22,E32,E13,E23,E33,GB1,GB2,GB3) CALL GAAPP(E13,E23,E33,G11,G12,G22,G13,G23,G33) END IF C Slowness vector P1=P1/SQRT(G) P2=P2/SQRT(G) P3=P3/SQRT(G) IF (LANISO) THEN C Recording minimum and maximum of the slowness vector ABSP=SQRT(P1*P1+P2*P2+P3*P3) IF (I.EQ.1) THEN ABSPMI=ABSP ABSPMA=ABSP ELSE IF (ABSP.LT.ABSPMI) ABSPMI=ABSP IF (ABSP.GT.ABSPMA) ABSPMA=ABSP ENDIF ENDIF C Ray-velocity vector V1=0.5*(G11*P1+G12*P2+G13*P3) V2=0.5*(G12*P1+G22*P2+G23*P3) V3=0.5*(G13*P1+G23*P2+G33*P3) C Surface to be triangularized (phase-slownes or ray-velocity) IF(KSURF.EQ.0) THEN C Phase-slownes surface XP(1,I)=P1 XP(2,I)=P2 XP(3,I)=P3 AUX=SQRT(V1*V1+V2*V2+V3*V3) Z1=V1/AUX Z2=V2/AUX Z3=V3/AUX ELSE C Ray-velocity surface XP(1,I)=V1 XP(2,I)=V2 XP(3,I)=V3 AUX=SQRT(P1*P1+P2*P2+P3*P3) Z1=P1/AUX Z2=P2/AUX Z3=P3/AUX END IF ZP(1,I)=Z1 ZP(2,I)=Z2 ZP(3,I)=Z3 IF(WA.NE.0.0.AND.WB.NE.0.0) THEN C Wave-propagation metric tensor WA=G/WA WB=G/WB G11=0.5*G11+WA*GA1*GA1+WB*GB1*GB1 G12=0.5*G12+WA*GA1*GA2+WB*GB1*GB2 G22=0.5*G22+WA*GA2*GA2+WB*GB2*GB2 G13=0.5*G13+WA*GA1*GA3+WB*GB1*GB3 G23=0.5*G23+WA*GA2*GA3+WB*GB2*GB3 G33=0.5*G33+WA*GA3*GA3+WB*GB3*GB3 C Signature of the curvature matrix DG=G11*G22*G33+2.*G12*G13*G23 * -G11*G23*G23-G22*G13*G13-G33*G12*G12 IF(DG.LE.0.0) THEN C Hyperbolic (signature 1,-1) KURV(I)=1 ELSE AUX=Z1*(G11*Z1+G12*Z2+G13*Z3) * +Z2*(G12*Z1+G22*Z2+G23*Z3) * +Z3*(G13*Z1+G23*Z2+G33*Z3) IF(G11+G22+G33.GT.AUX) THEN C Convex (signature 1,1) KURV(I)=0 ELSE C Concave (signature -1,-1) KURV(I)=2 END IF END IF ELSE KURV(I)=-1 END IF 13 CONTINUE ELSE DO 14 I=1,NP ZP(1,I)=XP(1,I) ZP(2,I)=XP(2,I) ZP(3,I)=XP(3,I) KURV(I)=0 14 CONTINUE END IF C C Writing anisotropy strength: IF (LANISO) THEN OPEN(LU,FILE=FILE) WRITE(LU,'(A)') ' /' FORMAT='(A,1(F00.0,A))' CALL FORM1(ABSPMI,ABSPMA,FORMAT(6:13)) FORMAT(13:13)=')' WRITE(LU,FORMAT) 'Minimum slowness: ',ABSPMI WRITE(LU,FORMAT) 'Maximum slowness: ',ABSPMA WRITE(LU,FORMAT) 'Anisotropy strength: ', * 2*(ABSPMA-ABSPMI)/(ABSPMA+ABSPMI) WRITE(LU,'(A)') ' /' CLOSE(LU) ENDIF C C Writing vertices: CALL RSEP3T('VRTX',FILE,'vrtx.out') OPEN(LU,FILE=FILE) WRITE(LU,'(A)') ' /' FORMAT='(A,I0.0,A,3(F00.0,A),3(F7.4,A),I2,A)' I=INT(ALOG10(FLOAT(NP)+0.5))+1 FORMAT(5:5)=CHAR(ICHAR('0')+I) FORMAT(7:7)=CHAR(ICHAR('0')+I) XPMIN=AMIN1(XP(1,1),XP(2,1),XP(3,1)) XPMAX=AMAX1(XP(1,1),XP(2,1),XP(3,1)) DO 21 I=2,NP XPMIN=AMIN1(XP(1,I),XP(2,I),XP(3,I),XPMIN) XPMAX=AMAX1(XP(1,I),XP(2,I),XP(3,I),XPMAX) 21 CONTINUE CALL FORM1(XPMIN,XPMAX,FORMAT(13:20)) FORMAT(20:20)=')' DO 22 I=1,NP WRITE(LU,FORMAT) '''',I,''' ', * XP(1,I),' ',XP(2,I),' ',XP(3,I),' ', * ZP(1,I),' ',ZP(2,I),' ',ZP(3,I),' ',KURV(I),' /' 22 CONTINUE C End of file WRITE(LU,'(A)') ' /' CLOSE(LU) C C Debugging vertices: * DO 23 I=1,NP * AUX=XP(1,I)**2+XP(2,I)**2+XP(3,I)**2 * WRITE(*,*) SQRT(AUX) * 23 CONTINUE * WRITE(*,*) C C Writing triangles: CALL RSEP3T('TRGL',FILE,'trgl.out') OPEN(LU,FILE=FILE) FORMAT='(3(I0,A))' I=INT(ALOG10(FLOAT(NT)+0.5))+1 FORMAT(5:5)=CHAR(ICHAR('0')+I) DO 31 I=1,NT WRITE(LU,FORMAT) KT(1,I),' ',KT(2,I),' ',KT(3,I),' /' 31 CONTINUE CLOSE(LU) C C Debugging triangles: * DO 32 I=1,NT * I1=KT(1,I) * I2=KT(2,I) * I3=KT(3,I) * AUX=(XP(1,I1)-XP(1,I2))**2+(XP(2,I1)-XP(2,I2))**2 * * +(XP(3,I1)-XP(3,I2))**2 * WRITE(*,*) SQRT(AUX) * AUX=(XP(1,I2)-XP(1,I3))**2+(XP(2,I2)-XP(2,I3))**2 * * +(XP(3,I2)-XP(3,I3))**2 * WRITE(*,*) SQRT(AUX) * AUX=(XP(1,I3)-XP(1,I1))**2+(XP(2,I3)-XP(2,I1))**2 * * +(XP(3,I3)-XP(3,I1))**2 * WRITE(*,*) SQRT(AUX) * 32 CONTINUE C RETURN END C C======================================================================= C SUBROUTINE ICOSAH(NDIVTR,KT,XP) INTEGER NDIVTR,KT(3,20*4**NDIVTR) REAL XP(3,10*4**NDIVTR+2) C INTEGER NT,NP,MP,IORDER,IT,JT,JE,J1,J2,J3,I1,I2,I REAL SQRT5,C1,S1,C2,S2,Z,R,X1,X2,X3,XX C C....................................................................... C C Regular icosahedron: C C Calculating the values of the coordinates of vertices: SQRT5=SQRT(5.) C1=0.25*(1.0+SQRT5) S1=1.0/SQRT(2.0+0.4*SQRT5) C2=1.0/(1.0+SQRT5) S2=2.0*C1*S1 Z=C1/(1.0+C1) R=SQRT(1.0+2.0*C1)/(1.0+C1) C1=R*C1 S1=R*S1 C2=R*C2 S2=R*S2 C C Vertices of the icosahedron: XP(1, 1)=0.0 XP(2, 1)=0.0 XP(3, 1)=1.0 XP(1, 2)=R XP(2, 2)=0.0 XP(3, 2)=Z XP(1, 3)=C2 XP(2, 3)=S2 XP(3, 3)=Z XP(1, 4)=-C1 XP(2, 4)=S1 XP(3, 4)=Z XP(1, 5)=-C1 XP(2, 5)=-S1 XP(3, 5)=Z XP(1, 6)=C2 XP(2, 6)=-S2 XP(3, 6)=Z XP(1, 7)=C1 XP(2, 7)=S1 XP(3, 7)=-Z XP(1, 8)=-C2 XP(2, 8)=S2 XP(3, 8)=-Z XP(1, 9)=-R XP(2, 9)=0.0 XP(3, 9)=-Z XP(1,10)=-C2 XP(2,10)=-S2 XP(3,10)=-Z XP(1,11)=C1 XP(2,11)=-S1 XP(3,11)=-Z XP(1,12)=0.0 XP(2,12)=0.0 XP(3,12)=-1.0 C C Triangular faces of the icosahedron: DO 10 I=1,5 KT(1, I)=1 KT(2, I)=1+I KT(3, I)=2+I KT(1, 5+I)=2+I KT(2, 5+I)=1+I KT(3, 5+I)=6+I KT(1,10+I)=1+I KT(2,10+I)=5+I KT(3,10+I)=6+I KT(1,15+I)=6+I KT(2,15+I)=5+I KT(3,15+I)=12 10 CONTINUE KT(3, 5)= 2 KT(1,10)= 2 KT(2,11)=11 KT(2,11)=11 KT(2,16)=11 C C....................................................................... C C Subdividing the regular icosahedron NDIVTR times: DO 90 IORDER=1,NDIVTR C Number of old triangles (number of new triangles is 4*NT) NT=20*4**(IORDER-1) C Number of old vertices NP=10*4**(IORDER-1)+2 C Number of new vertices MP=10*4**IORDER+2 C C Relocating old triangles in the memory: DO 20 IT=NT,1,-1 C Locations for new vertices KT(1,4*IT)=0 KT(2,4*IT)=0 KT(3,4*IT)=0 C Old vertices KT(1,4*IT-3)=KT(1,IT) KT(2,4*IT-3)=KT(2,IT) KT(3,4*IT-3)=KT(3,IT) 20 CONTINUE C C Dividing the edges of the triangles: C Loop over old triangles DO 34 JT=1,NT C Loop over the edges of the triangle DO 33 JE=1,3 IF(JE.EQ.1) THEN J1=2 J2=3 J3=1 ELSE IF(JE.EQ.2) THEN J1=3 J2=1 J3=2 ELSE J1=1 J2=2 J3=3 END IF C IF(KT(J3,4*JT).EQ.0) THEN C Dividing the edge IF(NP.GE.MP) THEN WRITE(*,*) ' Error NP.GE.MP',JE,JT,NP,MP END IF NP=NP+1 KT(J3,4*JT)=NP I1=KT(J1,4*JT-3) I2=KT(J2,4*JT-3) * WRITE(*,*) I1,I2,NP X1=XP(1,I1)+XP(1,I2) X2=XP(2,I1)+XP(2,I2) X3=XP(3,I1)+XP(3,I2) XX=SQRT(X1*X1+X2*X2+X3*X3) XP(1,NP)=X1/XX XP(2,NP)=X2/XX XP(3,NP)=X3/XX C C Searching for the other identical edge DO 31 IT=JT+1,NT IF(KT(1,4*IT-3).EQ.I1.AND.KT(2,4*IT-3).EQ.I2) THEN KT(3,4*IT)=NP GO TO 32 END IF IF(KT(2,4*IT-3).EQ.I1.AND.KT(1,4*IT-3).EQ.I2) THEN KT(3,4*IT)=NP GO TO 32 END IF IF(KT(2,4*IT-3).EQ.I1.AND.KT(3,4*IT-3).EQ.I2) THEN KT(1,4*IT)=NP GO TO 32 END IF IF(KT(3,4*IT-3).EQ.I1.AND.KT(2,4*IT-3).EQ.I2) THEN KT(1,4*IT)=NP GO TO 32 END IF IF(KT(3,4*IT-3).EQ.I1.AND.KT(1,4*IT-3).EQ.I2) THEN KT(2,4*IT)=NP GO TO 32 END IF IF(KT(1,4*IT-3).EQ.I1.AND.KT(3,4*IT-3).EQ.I2) THEN KT(2,4*IT)=NP GO TO 32 END IF 31 CONTINUE 32 CONTINUE END IF 33 CONTINUE 34 CONTINUE C C Creating new triangles: DO 80 IT=NT,1,-1 KT(1,4*IT-1)=KT(2,4*IT-3) KT(2,4*IT-1)=KT(1,4*IT ) KT(3,4*IT-1)=KT(3,4*IT ) KT(1,4*IT-2)=KT(3,4*IT-3) KT(2,4*IT-2)=KT(2,4*IT ) KT(3,4*IT-2)=KT(1,4*IT ) C KT(1,4*IT-3)=KT(1,4*IT-3) KT(2,4*IT-3)=KT(3,4*IT ) KT(3,4*IT-3)=KT(2,4*IT ) 80 CONTINUE C IF(NP.NE.MP) THEN C ANISRF-06 CALL ERROR('ANISRF-06: Wrong number of triangles') C This error should not appear. Contact the author. END IF C 90 CONTINUE RETURN 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 INCLUDE 'hder.for' C hder.for INCLUDE 'eigen.for' C eigen.for C C======================================================================= C