C
C Program ANISRF to triangularize the phase-slowness surface
C or the ray-velocity surface
C
C Version: 7.10
C Date: 2014, June 18
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 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
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
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,'animod.dat')
      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       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)
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 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