C 
C Program TRGLNORM to compute normals to given triangles
C
C Version: 5.50
C Date: 2000, September 4
C
C Coded by: Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: bulant@seis.karlov.mff.cuni.cz
C
C.......................................................................
C This program reads file TRGL with triangles (or polygons) and the
C corresponding file VRTX with coordinates of the vertices of the
C triangles (polygons).  Then it computes vector product W of the vector
C [first vertex , second vertex] with the vector
C [second vertex , third vertex].  If the scalar product of vector W
C with the input vector given by VECT1, VECT2 and VECT3 is nonnegative,
C the program writes the triangle (polygon) to the output file TRGLN,
C and the coordinates of vertices together with the vector W to the
C output file VRTXN.  The vertex being used in several triangles
C in file TRGL is thus written to file VRTXN several times with
C different normals W.
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 triangles.
C             Description of file VRTX
C             Default: VRTX='vrtx.out'
C     TRGL='string'... Name of the file with the triangles or polygons.
C             Description of file TRGL
C             Default: TRGL='trgl.out'
C Selection vector:
C     VECT1=real,VECT2=real,VECT3=real ... Three components of the
C             selection vector. The triangles (polygons) and the
C             corresponding vertices are written to the output files
C             only if the scalar product of the selection vector with
C             the normal to the triangle is nonnegative.
C             Default: VECT1=0.,VECT2=0.,VECT3=0. means that all the
C               triangles are written to the output files.
C Data specifying output files:
C     VRTXN='string'... Name of the file with vertices of the triangles.
C             Description of file VRTXN
C             Default: VRTXN='vrtxn.out'
C     TRGLN='string'... Name of the file with the triangles or polygons.
C             Description of file TRGLN
C             Default: TRGLN='trgln.out'
C Data specifying form of output files:
C     KOLUM1=integer,KOLUM2=integer,KOLUM3=integer ... Indices of the
C             columns, where to write the components of the computed
C             normal vector.
C             Default: KOLUM1=4,KOLUM2=5,KOLUM3=6
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,/
C     'NAME'... Name of the vertex.  Not considered.  May be blank.
C     X1,X2,X3... Coordinates of the vertex.
C     /...    None to several values terminated by a slash.
C (3) / or end of file.
C
C                                                    
C Input 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 (or more in case of polygon) vertices
C             of the triangle.
C             The vertices in file VRTX are indexed by positive integers
C             according to their order.
C     /...    List of vertices of the triangle is terminated by a slash.
C
C                                                   
C Output file VRTXN with the vertices:
C (1) / (a slash)
C (2) For each vertex data (2.1):
C (2.1) 'NAME',X1,X2,X3,W1,W2,W3,/
C     'NAME'..Name of the vertex.  String in apostrophes containing
C             the index of the vertex corresponding to file TRGLN.
C     X1,X2,X3,W1,W2,W3... KOLUMi-th column of the input file is
C             replaced by the value of i-th component of the normal
C             to the corresponding triangle in file TRGLN.
C             If the input file VRTX contains less than
C             MAX(KOLUM1,KOLUM2,KOLUM3) values, the missing values
C             are replaced by the value of parameter UNDEF, see below.
C     /  ...    A slash.
C (3) / (a slash)
C
C                                                   
C Output file TRGLN with the triangles:
C (1) For each triangle data (1.1):
C (1.1) I1,I2,I3,/
C     I1,I2,I3... Indices of 3 (or more in case of polygon) vertices
C             of the triangle.
C             The vertices in file VRTXN are indexed by positive
C             integers according to their order.
C     /...    List of vertices of the triangle 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,FORM1,WARN,RSEP1,RSEP3T,RSEP3R,RSEP3I
      INTEGER  LENGTH
C
C.......................................................................
C
C     Filenames and parameters:
      CHARACTER*80 FSEP,FVRTX,FTRGL,FVRTXN,FTRGLN
      INTEGER LU,LU1,LU2,IUNDEF,MVRTX
      PARAMETER (LU=1,LU1=2,LU2=3,IUNDEF=-999999,MVRTX=1000)
      REAL UNDEF
      PARAMETER (UNDEF=9.9E9)
C     Input data:
      REAL VECT1,VECT2,VECT3
      CHARACTER*10 FORMA1
      CHARACTER*26 FORMA2
      CHARACTER*80 TEXT
C     Other variables:
      INTEGER NVRTX,NBIGON,N,I,I1,I2,J1,J2,J3,NQ
      INTEGER NPTS,KOLUM1,KOLUM2,KOLUM3
      REAL A1,A2,A3,B1,B2,B3,W1,W2,W3,AUX,OUTMIN,OUTMAX
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     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)') '+TRGLNORM: Enter input filename: '
      FSEP=' '
      READ (*,*) FSEP
      WRITE(*,'(A)') '+TRGLNORM: Working...            '
C
C     Reading all data from the SEP file into the memory:
      IF (FSEP.NE.' ') THEN
        CALL RSEP1(LU,FSEP)
      ELSE
C       TRGLNORM-01
        CALL ERROR('TRGLNORM-01: 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('TRGL',FTRGL,'trgl.out')
      CALL RSEP3T('VRTXN',FVRTXN,'vrtxn.out')
      CALL RSEP3T('TRGLN',FTRGLN,'trgln.out')
C
C     Reading the selection vector:
      CALL RSEP3R('VECT1',VECT1,0.)
      CALL RSEP3R('VECT2',VECT2,0.)
      CALL RSEP3R('VECT3',VECT3,0.)
C
C     Reading the columns where to write the normal:
      CALL RSEP3I('KOLUM1',KOLUM1,4)
      CALL RSEP3I('KOLUM2',KOLUM2,5)
      CALL RSEP3I('KOLUM3',KOLUM3,6)
C
C     Reading vertices:
      OPEN(LU,FILE=FVRTX,STATUS='OLD')
      READ(LU,*) (TEXT,I=1,20)
      NVRTX=0
      NQ=MAX0(3,KOLUM1,KOLUM2,KOLUM3)
   10 CONTINUE
        IF(NVRTX+NQ.GT.MRAM) THEN
C         TRGLNORM-02
          CALL ERROR('TRGLNORM-02: Too small array RAM')
        END IF
        DO 12, I=1,NQ
          RAM(NVRTX+I)=UNDEF
  12    CONTINUE
        TEXT='$'
        READ(LU,*,END=18) TEXT,(RAM(I),I=NVRTX+1,NVRTX+NQ)
        IF(TEXT.EQ.'$') THEN
          GO TO 18
        END IF
        NVRTX=NVRTX+NQ
      GO TO 10
   18 CONTINUE
      CLOSE(LU)
C
C     Output format for the file with polygons
      FORMA1='(99(I0,A))'
      FORMA2='(A,1I0.0,A,00(F00.0,1X),A)'
      I=INT(ALOG10(FLOAT(NVRTX/NQ)))+1
      IF (I.GT.9) THEN
C       TRGLNORM-03
        CALL ERROR('TRGLNORM-03: Too many vertices in file VRTX')
C       This format specification allows for maximum of 100 000 000
C       of vertices in file VRTX
      ENDIF
      FORMA1(6:6)=CHAR(ICHAR('0')+I)
      FORMA2(6:6)=FORMA1(6:6)
      FORMA2(8:8)=FORMA1(6:6)
      FORMA2(13:13)=CHAR(ICHAR('0')+MOD(NQ/1,10))
      FORMA2(12:12)=CHAR(ICHAR('0')+MOD(NQ/10,10))
C
C     Reading polygons:
      DO 19 I=NVRTX+1,MRAM
        IRAM(I)=0
   19 CONTINUE
      OPEN(LU,FILE=FTRGL,STATUS='OLD')
      IF (FVRTXN.NE.' ') OPEN(LU1,FILE=FVRTXN)
      IF (FTRGLN.NE.' ') OPEN(LU2,FILE=FTRGLN)
      IF (FVRTXN.NE.' ') WRITE(LU1,'(A)') '/'
      NPTS=0
      NBIGON=0
C     Loop over the polygons:
   20 CONTINUE
        IRAM(NVRTX+1)=IUNDEF
        READ(LU,*,END=29) (IRAM(I),I=NVRTX+1,MIN0(NVRTX+MVRTX+1,MRAM))
        IF(IRAM(NVRTX+1).EQ.IUNDEF) THEN
          GO TO 29
        END IF
        DO 21 I=NVRTX+1,MIN0(NVRTX+MVRTX+1,MRAM)
          IF(IRAM(I).LE.0) THEN
C           Number of polygon vertices
            N=I-1-NVRTX
            GO TO 22
          ELSE IF(IRAM(I).GT.NVRTX/NQ) THEN
C           TRGLNORM-04
            WRITE(TEXT,'(A,I6)')
     *      'TRGLNORM-04: Wrong vertex index:',IRAM(I)
            CALL ERROR(TEXT(1:LENGTH(TEXT)))
          END IF
   21   CONTINUE
C         TRGLNORM-05
          CALL ERROR('TRGLNORM-05: Too many vertices in polygons')
   22   CONTINUE
        IF(N.LT.3) THEN
C         TRGLNORM-06
          CALL ERROR('TRGLNORM-06: Polygon of less than 3 vertices')
        END IF
C       Check for identical indices of vertices:
        DO 24 I2=NVRTX+1,NVRTX+N
          DO 23 I1=I2+1,NVRTX+N
            IF(IRAM(I2).EQ.IRAM(I1)) THEN
C             TRGLNORM-07
              WRITE(TEXT,'(A,I6)')
     *       'TRGLNORM-07: 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 coordinates of vertices:
   25   CONTINUE
        DO 27 I2=NVRTX+1,NVRTX+N
          J1=NQ*(IRAM(I2)-1)
          IF(I2.LT.NVRTX+N) THEN
            J2=NQ*(IRAM(I2+1)-1)
          ELSE
            J2=NQ*(IRAM(NVRTX+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
C           Two subsequent points of a polygon have the same coordinates
            DO 26 I1=I2+1,NVRTX+N-1
              IRAM(I1)=IRAM(I1+1)
   26       CONTINUE
            IRAM(NVRTX+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
        J1=NQ*(IRAM(NVRTX+1)-1)
        J2=NQ*(IRAM(NVRTX+2)-1)
        J3=NQ*(IRAM(NVRTX+3)-1)
        A1=RAM(J2+1)-RAM(J1+1)
        A2=RAM(J2+2)-RAM(J1+2)
        A3=RAM(J2+3)-RAM(J1+3)
        B1=RAM(J3+1)-RAM(J1+1)
        B2=RAM(J3+2)-RAM(J1+2)
        B3=RAM(J3+3)-RAM(J1+3)
        W1=A2*B3-B2*A3
        W2=A3*B1-B3*A1
        W3=A1*B2-B1*A2
C       Normalizing the normal:
        AUX=SQRT(W1**2+W2**2+W3**2)
        W1=W1/AUX
        W2=W2/AUX
        W3=W3/AUX
C       Scalar product with selection vector:
        AUX=VECT1*W1+VECT2*W2+VECT3*W3
        IF (AUX.GE.0.) THEN
          IF (FTRGLN.NE.' ') THEN
C           Writing the polygon:
            WRITE(LU2,FORMA1)   (NPTS+I,' ',I=1,N-1),NPTS+N,' /'
          ENDIF
          IF (FVRTXN.NE.' ') THEN
C           Writing the vertices:
            DO 28, I1=NVRTX+1,NVRTX+N
              J1=NQ*(IRAM(I1)-1)
              RAM(J1+KOLUM1)=W1
              RAM(J1+KOLUM2)=W2
              RAM(J1+KOLUM3)=W3
              OUTMIN=0.
              OUTMAX=0.
              DO 16, I=J1+1,J1+NQ
                IF(RAM(I).LT.OUTMIN) OUTMIN=RAM(I)
                IF(RAM(I).GT.OUTMAX) OUTMAX=RAM(I)
  16          CONTINUE
              CALL FORM1(OUTMIN,OUTMAX,FORMA2(15:22))
              FORMA2(21:24)= '1X),'
              NPTS=NPTS+1
              WRITE(LU1,FORMA2)
     *              ' ''',NPTS,''' ',(RAM(J1+I),I=1,NQ),'/'
  28        CONTINUE
          ENDIF
        ENDIF
      GO TO 20
C     End of the loop over the polygons.
   29 CONTINUE
      IF(NBIGON.GT.0) THEN
C       TRGLNORM-08
        WRITE(TEXT,'(A,3I6)')
     *    'TRGLNORM-08: 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)
      IF (FVRTXN.NE.' ') WRITE(LU1,'(A)') '/'
      IF (FVRTXN.NE.' ') CLOSE(LU1)
      IF (FTRGLN.NE.' ') CLOSE(LU2)
      WRITE(*,'(A)') '+TRGLNORM: Done.                 '
C
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C