C 
C Program TRGL to divide polygons on a curved surface into triangles,
C right-handed with respect to the surface normals
C
C Version: 5.40
C Date: 2000, January 21
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: klimes@seis.karlov.mff.cuni.cz
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 files:
C     VRTX='string'... Name of the file with vertices of the polygons.
C             Description of file VRTX
C             Default: VRTX='vrtx.out'
C     PLGN='string'... Name of the file describing the polygons.
C             Description of file PLGN
C             Default: PLGN='plgn.out'
C Data specifying output file:
C     TRGL='string'... Name of the file describing the triangles.
C             Description of file TRGL
C             Default: TRGL='trgl.out'
C
C                                                    
C Input file VRTX with the vertices:
C (1) None to several strings terminated by / (a slash)
C (2) For each vertex data (2.1):
C (2.1) 'NAME',X1,X2,X3,Z1,Z2,Z3,/
C     'NAME'... Name of the vertex.  Not considered.  May be blank.
C     X1,X2,X3... Coordinates of the vertex.
C     Z1,Z2,Z3... Normal to the surface at the vertex.
C     /...    None to several values terminated by a slash.
C (3) / or end of file.
C
C                                                    
C Input file PLGN with the polygons:
C (1) For each polygon data (1.1):
C (1.1) I1,I2,...,IN,/
C     I1,I2,...,IN... Indices of N vertices of the polygon.
C             The vertices in file VRTX are indexed by positive integers
C             according to their order.
C     /...    List of vertices must be terminated by a slash.
C (2) / or end of file.
C
C                                                    
C Output file TRGL with the triangles:
C (1) For each polygon 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 LENGTH,ERROR,WARN,RSEP1,RSEP3T
      INTEGER  LENGTH
C
C.......................................................................
C
C     Filenames and parameters:
      CHARACTER*80 FSEP,FVRTX,FPLGN,FTRGL
      INTEGER LU,IUNDEF,MVRTX
      PARAMETER (LU=1,IUNDEF=-999999,MVRTX=1000)
C     Input data:
      CHARACTER*9 FORMAT
      CHARACTER*80 TEXT
C     Other variables:
      INTEGER NVRTX,NPLGN,NTRGL,NBIGON,N,I,I1,I2,I3,J1,J2,J3,K1,K2,N1,N2
      REAL X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,A1,A2,A3,D,DMAX,DMIN,HAND,S,C,AUX
C
C     MVRTX...Maximum number of vertices of one polygon.
C     NVRTX...Number of storage locations for the vertices, i.e. 6 times
C             the number of vertices.
C     NPLGN...Last storage location for polygons, i.e. NVRTX + 4 times
C             the number of future triangles.  Each polygon with N
C             vertices occupies 4*(N-2) locations because it will be
C             split into N-2 triangles.  4*(N-2)-N locations following
C             N vertex indices are filled wth zeros.
C     NTRGL...Last storage location with triangles, i.e. NVRTX + 4 times
C             the number of created triangles.
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+TRGL: Enter input filename: '
      FSEP=' '
      READ (*,*) FSEP
      WRITE(*,'(A)') '+TRGL: Working...            '
C
C     Reading all data from the SEP file into the memory:
      IF (FSEP.NE.' ') THEN
        CALL RSEP1(LU,FSEP)
      ELSE
C       TRGL-11
        CALL ERROR('TRGL-11: SEP file not given')
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.
      ENDIF
C
C     Reading input and output filenames:
      CALL RSEP3T('VRTX',FVRTX,'vrtx.out')
      CALL RSEP3T('PLGN',FPLGN,'plgn.out')
      CALL RSEP3T('TRGL',FTRGL,'trgl.out')
C
C     Reading vertices:
      OPEN(LU,FILE=FVRTX)
      READ(LU,*) (TEXT,I=1,20)
      NVRTX=0
   10 CONTINUE
        IF(NVRTX+6.GT.MRAM) THEN
C         TRGL-01
          CALL ERROR('TRGL-01: Too small array RAM')
        END IF
        TEXT='$'
        READ(LU,*,END=18) TEXT,(RAM(I),I=NVRTX+1,NVRTX+6)
        IF(TEXT.EQ.'$') THEN
          GO TO 18
        END IF
C       Normalizing the normal
        AUX=SQRT(RAM(NVRTX+4)**2+RAM(NVRTX+5)**2+RAM(NVRTX+6)**2)
        RAM(NVRTX+4)=RAM(NVRTX+4)/AUX
        RAM(NVRTX+5)=RAM(NVRTX+5)/AUX
        RAM(NVRTX+6)=RAM(NVRTX+6)/AUX
        NVRTX=NVRTX+6
      GO TO 10
   18 CONTINUE
      CLOSE(LU)
C
C     Reading polygons:
      DO 19 I=NVRTX+1,MRAM
        IRAM(I)=0
   19 CONTINUE
      OPEN(LU,FILE=FPLGN)
      NPLGN=NVRTX
      NBIGON=0
   20 CONTINUE
        IRAM(NPLGN+1)=IUNDEF
        READ(LU,*,END=29) (IRAM(I),I=NPLGN+1,MIN0(NPLGN+MVRTX+1,MRAM))
        IF(IRAM(NPLGN+1).EQ.IUNDEF) THEN
          GO TO 29
        END IF
        DO 21 I=NPLGN+1,MIN0(NPLGN+MVRTX+1,MRAM)
          IF(IRAM(I).LE.0) THEN
C           Number of polygon vertices
            N=I-1-NPLGN
            GO TO 22
          ELSE IF(IRAM(I).GT.NVRTX/6) THEN
C           TRGL-02
            WRITE(TEXT,'(A,I6)') 'TRGL-02: Wrong vertex index:',IRAM(I)
            CALL ERROR(TEXT(1:LENGTH(TEXT)))
          END IF
   21   CONTINUE
C         TRGL-03
          CALL ERROR('TRGL-03: Too many vertices in polygons')
   22   CONTINUE
        IF(N.LT.3) THEN
C         TRGL-04
          CALL ERROR('TRGL-04: Polygon of less than 3 vertices')
        END IF
        IF(NPLGN+4*(N-2).GT.MRAM) THEN
C         TRGL-05
          CALL ERROR('TRGL-05: Too small array RAM')
        END IF
C       Checking vertex indices:
        DO 24 I2=NPLGN+1,NPLGN+N
          DO 23 I1=I2+1,NPLGN+N
          IF(IRAM(I2).EQ.IRAM(I1)) THEN
C           TRGL-06
            WRITE(TEXT,'(A,I6)')
     *        'TRGL-06: The same vertex twice in a polygon:',IRAM(I2)
            CALL ERROR(TEXT(1:LENGTH(TEXT)))
C           All vertices of a polygon must have different indices.
          END IF
   23     CONTINUE
   24   CONTINUE
C       Check for identical indices:
   25   CONTINUE
        DO 27 I2=NPLGN+1,NPLGN+N
          J1=6*(IRAM(I2)-1)
          IF(I2.LT.NPLGN+N) THEN
            J2=6*(IRAM(I2+1)-1)
          ELSE
            J2=6*(IRAM(NPLGN+1)-1)
          END IF
          IF(RAM(J1+1).EQ.RAM(J2+1).AND.
     *       RAM(J1+2).EQ.RAM(J2+2).AND.
     *       RAM(J1+3).EQ.RAM(J2+3)) THEN
            IF(RAM(J1+1).NE.RAM(J2+1).OR.
     *         RAM(J1+2).NE.RAM(J2+2).OR.
     *         RAM(J1+3).NE.RAM(J2+3)) THEN
C             TRGL-10
              WRITE(TEXT,'(A,3I6)')
     *          'TRGL-10: Different normals at identical points:',
     *          IRAM(I2),IRAM(I2+1)
              CALL ERROR(TEXT(1:LENGTH(TEXT)))
C             Two subsequent points of a polygon have the same
C             coordinates but different normals.
            END IF
            DO 26 I1=I2+1,NPLGN+N-1
              IRAM(I1)=IRAM(I1+1)
   26       CONTINUE
            IRAM(NPLGN+N)=0
            N=N-1
            GO TO 25
          END IF
   27   CONTINUE
        IF(N.LT.3) THEN
          NBIGON=NBIGON+1
          GO TO 20
        END IF
C       Leaving the space for N-2 triangles, each terminated by zero:
        DO 28 I=NPLGN+N+1,NPLGN+4*(N-2)
          IRAM(I)=0
   28   CONTINUE
        NPLGN=NPLGN+4*(N-2)
      GO TO 20
   29 CONTINUE
      IF(NBIGON.GT.0) THEN
C       TRGL-53
        WRITE(TEXT,'(A,3I6)')
     *    'TRGL-53: Number of polygons with less than 3 vertices:',
     *    NBIGON
        CALL WARN(TEXT(1:LENGTH(TEXT)))
C       Polygons with less than 3 vertices are created from polygons
C       with subsequent vertices of identical coordinates.
      END IF
      CLOSE(LU)
C
C     Dividing polygons into triangles:
      IF(NVRTX.EQ.NPLGN) THEN
C       There are no polygons
        GO TO 49
      END IF
      NTRGL=NVRTX
C     Loop over polygons
   30 CONTINUE
C       Determining number N of polygon vertices
        DO 31 I=NTRGL+1,NPLGN
          IF(IRAM(I).LE.0) THEN
            N=I-1-NTRGL
            GO TO 32
          END IF
   31   CONTINUE
   32   CONTINUE
        IF(N.EQ.3) THEN
C         Current polygon is a triangle, proceeding to the next polygon
          NTRGL=NTRGL+4
          IF(NTRGL.LT.NPLGN) THEN
            GO TO 30
          ELSE
C           All polygons are divided into triangles
            GO TO 49
          END IF
        END IF
C
C       Dividing the polygon into 2 polygons:
        K1=0
        DMIN=1.0E20
C       Loop over polygon vertices
C       (a) Determination of the handedness of the polygon
        HAND=0.
        DO 35 I2=1,N
          I3=MOD(I2,N)+1
          IF(I2.GT.1) THEN
            I1=I2-1
          ELSE
            I1=N
          END IF
          J1=6*(IRAM(NTRGL+I1)-1)
          J2=6*(IRAM(NTRGL+I2)-1)
          J3=6*(IRAM(NTRGL+I3)-1)
          X1=RAM(J2+1)-RAM(J1+1)
          X2=RAM(J2+2)-RAM(J1+2)
          X3=RAM(J2+3)-RAM(J1+3)
          Y1=RAM(J3+1)-RAM(J2+1)
          Y2=RAM(J3+2)-RAM(J2+2)
          Y3=RAM(J3+3)-RAM(J2+3)
          Z1=X2*Y3-X3*Y2
          Z2=X3*Y1-X1*Y3
          Z3=X1*Y2-X2*Y1
          S=Z1*RAM(J2+4)+Z2*RAM(J2+5)+Z3*RAM(J2+6)
          C=X1*Y1+X2*Y2+X3*Y3
          HAND=HAND+ATAN2(S,C)
   35   CONTINUE
C       (b) Selection of the triangle to be separated
        DO 37 I2=1,N
          I3=MOD(I2,N)+1
          IF(I2.GT.1) THEN
            I1=I2-1
          ELSE
            I1=N
          END IF
          J1=6*(IRAM(NTRGL+I1)-1)
          J2=6*(IRAM(NTRGL+I2)-1)
          J3=6*(IRAM(NTRGL+I3)-1)
          X1=RAM(J2+1)-RAM(J1+1)
          X2=RAM(J2+2)-RAM(J1+2)
          X3=RAM(J2+3)-RAM(J1+3)
          Y1=RAM(J3+1)-RAM(J2+1)
          Y2=RAM(J3+2)-RAM(J2+2)
          Y3=RAM(J3+3)-RAM(J2+3)
          Z1=X2*Y3-X3*Y2
          Z2=X3*Y1-X1*Y3
          Z3=X1*Y2-X2*Y1
          S=Z1*RAM(J2+4)+Z2*RAM(J2+5)+Z3*RAM(J2+6)
          IF(S*HAND.GT.0.) THEN
C           Normal Z to the separation plane
            X1=RAM(J3+1)-RAM(J1+1)
            X2=RAM(J3+2)-RAM(J1+2)
            X3=RAM(J3+3)-RAM(J1+3)
            Y1=RAM(J3+4)+RAM(J1+4)
            Y2=RAM(J3+5)+RAM(J1+5)
            Y3=RAM(J3+6)+RAM(J1+6)
            Z1=X2*Y3-X3*Y2
            Z2=X3*Y1-X1*Y3
            Z3=X1*Y2-X2*Y1
C           Point on the diagonal
            Y1=RAM(J1+1)
            Y2=RAM(J1+2)
            Y3=RAM(J1+3)
C           Vertex to be separated
            D=Z1*(RAM(J2+1)-Y1)+Z2*(RAM(J2+2)-Y2)+Z3*(RAM(J2+3)-Y3)
            DMAX=0.
            IF(D*HAND.GT.0) THEN
C             Loop over opposite vertices
              DO 36 I=1,N
                IF(I.NE.I1.AND.I.NE.I2.AND.I.NE.I3) THEN
                  J2=6*(IRAM(NTRGL+I)-1)
                  AUX=Z1*(RAM(J2+1)-Y1)+Z2*(RAM(J2+2)-Y2)
     *                                 +Z3*(RAM(J2+3)-Y3)
                  DMAX=AMAX1(-AUX*SIGN(1.,HAND),DMAX)
                END IF
   36         CONTINUE
              AUX=AMIN1(ABS(D),DMAX)
              IF(AUX.GT.0.) THEN
                AUX=ABS((RAM(J3+4)-RAM(J1+4))*X1
     *                 +(RAM(J3+5)-RAM(J1+5))*X2
     *                 +(RAM(J3+6)-RAM(J1+6))*X3)
     *              *SQRT(X1*X1+X2*X2+X3*X3)/AUX
                IF(K1.EQ.0.OR.AUX.LT.DMIN) THEN
                  K1=MIN0(I1,I3)
                  K2=MAX0(I1,I3)
                  DMIN=AUX
                END IF
              END IF
            END IF
          END IF
   37   CONTINUE
C       ^^^^^^^^
C       Loop over polygon diagonals to find the best one for division
*521    DO 42 I1=1,N-2
*       DO 42 I1=1,N
*         J1=6*(IRAM(NTRGL+I1)-1)
C521      Loop over all polygon diagonals
*521      DO 41 I2=I1+2,N+MIN0(I1-2,0)
C-NEW       ***
C         Limiting the loop to a single diagonal per vertex
*         I2=MOD(I1+1,N)+1
*           J2=6*(IRAM(NTRGL+I2)-1)
C           Normal Z to the separation plane
*           X1=RAM(J2+1)-RAM(J1+1)
*           X2=RAM(J2+2)-RAM(J1+2)
*           X3=RAM(J2+3)-RAM(J1+3)
*           Y1=RAM(J2+4)+RAM(J1+4)
*           Y2=RAM(J2+5)+RAM(J1+5)
*           Y3=RAM(J2+6)+RAM(J1+6)
*           Z1=X2*Y3-X3*Y2
*           Z2=X3*Y1-X1*Y3
*           Z3=X1*Y2-X2*Y1
C           Point on the diagonal
*           Y1=RAM(J1+1)
*           Y2=RAM(J1+2)
*           Y3=RAM(J1+3)
C           Vertex to be separated
*           I3=MOD(I1,N)+1
*           J3=6*(IRAM(NTRGL+I3)-1)
*           D=Z1*(RAM(J3+1)-Y1)+Z2*(RAM(J3+2)-Y2)+Z3*(RAM(J3+3)-Y3)
*           DMAX=0.
*           IF(D.NE.0) THEN
C             Loop over opposite vertices
*             DO 35 I=1,N
*               IF(I.NE.I1.AND.I.NE.I2.AND.I.NE.I3) THEN
*                 J3=6*(IRAM(NTRGL+I)-1)
*                 AUX=Z1*(RAM(J3+1)-Y1)+Z2*(RAM(J3+2)-Y2)
*    *                                 +Z3*(RAM(J3+3)-Y3)
*                 IF(AUX*D.GE.0.) THEN
C                   Vertex I3 is not separated from opposite vertices
*                   GO TO 41
*                 END IF
*                 DMAX=AMAX1(ABS(AUX),DMAX)
*               END IF
*  35         CONTINUE
*             AUX=SQRT((RAM(J2+4)-RAM(J1+4))**2
*    *                +(RAM(J2+5)-RAM(J1+5))**2
*    *                +(RAM(J2+6)-RAM(J1+6))**2)
*    *            *(X1*X1+X2*X2+X3*X3)/AMIN1(ABS(D),DMAX)
*             IF(K1.EQ.0.OR.AUX.LT.DMIN) THEN
*               K1=MIN0(I1,I2)
*               K2=MAX0(I1,I2)
*               DMIN=AUX
*             END IF
*           END IF
*  41     CONTINUE
*  42   CONTINUE
C       ^^^^^^^^
        IF(K1.EQ.0) THEN
C         TRGL-09
          WRITE(TEXT,'(A,7I6)')
     *      'TRGL-09: Polygon cannot be divided:',
     *      (IRAM(I),I=NTRGL+1,NTRGL+MIN0(7,N))
          IF(N.GT.7) THEN
            TEXT(78:80)='...'
          END IF
          CALL ERROR(TEXT(1:LENGTH(TEXT)))
C         No polygon vertex can be separated from all opposite polygon
C         vertices by the plane determined by the diagonal connecting
C         the two adjacent vertices and the sum of the normals in the
C         adjacent vertices.
        END IF
C       v---521---v
C       Dividing polygon along diagonal between K1-th and K2-th vertices
        N1=K1+N-K2+1
        N2=K2-K1+1
C       Doubling K2-th and K1-th vertices
        DO 43 I=NTRGL+N,NTRGL+K2,-1
          IRAM(I+1)=IRAM(I)
   43   CONTINUE
        DO 44 I=NTRGL+N+1,NTRGL+K1,-1
          IRAM(I+1)=IRAM(I)
   44   CONTINUE
C       Moving vertices K2,K2+1,...,N
C       between vertices 1,2,...,K1 and K1,K1+1,...,K2
        DO 46 I2=K2,N
          I=IRAM(NTRGL+N+2)
          DO 45 I1=NTRGL+N+1,NTRGL+K1+1,-1
            IRAM(I1+1)=IRAM(I1)
   45     CONTINUE
          IRAM(NTRGL+K1+1)=I
   46   CONTINUE
C       Moving the second polygon to the proper location
        DO 47 I=NTRGL+N2,NTRGL+1,-1
          IRAM(I+4*(N1-2))=IRAM(I+N1)
   47   CONTINUE
C       Storing zeros between the polygons
        DO 48 I=NTRGL+4*(N1-2),NTRGL+N1+1,-1
          IRAM(I)=0
   48   CONTINUE
      GO TO 30
   49 CONTINUE
C
C     Making triangles right-handed with respect to the normals:
      DO 51 I=NVRTX+1,NPLGN,4
        IF(IRAM(I+3).NE.0) THEN
C         TRGL-07
          CALL ERROR('TRGL-07: Triangle not terminated by 0')
C         This errror should not appear.  Contact the author.
        END IF
        J1=6*(IRAM(I  )-1)
        J2=6*(IRAM(I+1)-1)
        J3=6*(IRAM(I+2)-1)
        X1=RAM(J2+1)-RAM(J1+1)
        X2=RAM(J2+2)-RAM(J1+2)
        X3=RAM(J2+3)-RAM(J1+3)
        Y1=RAM(J3+1)-RAM(J2+1)
        Y2=RAM(J3+2)-RAM(J2+2)
        Y3=RAM(J3+3)-RAM(J2+3)
        Z1=X2*Y3-X3*Y2
        Z2=X3*Y1-X1*Y3
        Z3=X1*Y2-X2*Y1
        IF(Z1.EQ.0..AND.Z2.EQ.0..AND.Z3.EQ.0.) THEN
C         TRGL-51
          WRITE(TEXT,'(A,3I6)')
     *      'TRGL-51: Triangle has shrunk into a line:',
     *      IRAM(I),IRAM(I+1),IRAM(I+2)
          CALL WARN(TEXT(1:LENGTH(TEXT)))
C         The sides of a triangle are parallel.
C         Marking the triangle linear in order not to write it into the
C         output file:
          IRAM(I+3)=-1
        ELSE
          A1=Z1*RAM(J1+4)+Z2*RAM(J1+5)+Z3*RAM(J1+6)
          A2=Z1*RAM(J2+4)+Z2*RAM(J2+5)+Z3*RAM(J2+6)
          A3=Z1*RAM(J3+4)+Z2*RAM(J3+5)+Z3*RAM(J3+6)
          IF(A1.LT.0..AND.A2.LT.0..AND.A3.LT.0.) THEN
C           Changing left-handed triangle into right-handed one
            J2=IRAM(I+1)
            J3=IRAM(I+2)
            IRAM(I+1)=J3
            IRAM(I+2)=J2
          ELSE IF(A1.LE.0..OR.A2.LE.0..OR.A3.LE.0.) THEN
            AUX=X1*(X1+Y1)+Y1*Y1+X2*(X2+Y2)+Y2*Y2+X3*(X3+Y3)+Y3*Y3
            AUX=SQRT(Z1*Z1+Z2*Z2+Z3*Z3)/AUX
            IF(AUX.LE.0.000010) THEN
C             TRGL-52
              WRITE(TEXT,'(A,3I6)')
     *          'TRGL-52: Triangle is too narrow:',
     *          IRAM(I),IRAM(I+1),IRAM(I+2)
              CALL WARN(TEXT(1:LENGTH(TEXT)))
C             Triangle is too narrow, i.e., the area of the triangle
C             is too small compared with the sum of the squares of its
C             sides.
C             Marking the triangle linear in order not to write it into
C             the output file:
              IRAM(I+3)=-1
            ELSE
C             TRGL-08
              WRITE(TEXT,'(A,3I6)')
     *          'TRGL-08: Wrong normals at the triangle vertices:',
     *          IRAM(I),IRAM(I+1),IRAM(I+2)
              CALL ERROR(TEXT(1:LENGTH(TEXT)))
C             Normals at the triangle vertices do not point at the same
C             side of the triangle.
            END IF
          END IF
        END IF
   51 CONTINUE
C
C     Writing triangles:
      OPEN(LU,FILE=FTRGL)
      FORMAT='(3(I0,A))'
      I=INT(ALOG10(FLOAT(NVRTX)+0.5))+1
      FORMAT(5:5)=CHAR(ICHAR('0')+I)
      DO 61 I=NVRTX+1,NPLGN,4
        IF(IRAM(I+3).EQ.0) THEN
          WRITE(LU,FORMAT) IRAM(I),' ',IRAM(I+1),' ',IRAM(I+2),' /'
        END IF
   61 CONTINUE
      CLOSE(LU)
C
      IF(NVRTX.EQ.NPLGN) THEN
C       There are no polygons
        WRITE(*,'(A)') '+TRGL: No triangles.         '
      ELSE
        WRITE(*,'(A)') '+TRGL: Done.                 '
      END IF
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C