C
C Program SRFSECT to calculate the crossection of a triangulated surface
C with a parallelogram.  Each of the three sides of each triangle of the
C triangulated surface is tested whether it has a crossection with the
C specified parallelogram. The crossection points, if found, are written
C to the output file.
C
C Version: 7.50
C Date: 2018, December 11
C
C Coded by: Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/bulant.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 files with the triangulated surface:
C     VRTX='string'... Name of the file with vertices of the polygons.
C             Description of file VRTX
C             Default: VRTX='vrtx.out'
C     TRGL='string'... Name of the file describing the triangles or
C             polygons.
C             Description of file TRGL
C             Default: TRGL='trgl.out'
C Data specifying the intersection parallelogram:
C     A1=real,A2=real,A3=real,B1=real,B2=real,B3=real,C1=real,C2=real,
C     C3=real... Cartesian coordinates of three points A,B and C
C             defining the intersection plane. The intersection of the
C             given triangulated surface with the parallelogram given
C             by points A,B and C is searched for.
C             Defaults: A1=0.,A2=0.,A3=0.,B1=0.,B2=0.,B3=0.,
C               C1=0.,C2=0.,C3=0.
C             Note that the points A,B and C must be different and thus
C             at least two of the values A1,A2,A3,B1,B2,B3,C1,C2,C3 must
C             be different from the default values.
C Output file:
C     SRFSECT='string'... Name of the file to store the calculated
C             crossection of the surface with the parallelogram.
C             Each of the three sides of each triangle from file TRGL
C             is tested whether it has a crossection with the specified
C             parallelogram. The crossection points, if found, are
C             written to this output file.
C             Description of file SRFSECT.
C             Default: SRFSECT='srfsect.out'
C Optional parameters specifying the form of the real quantities
C written in the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             forms.for.
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 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 (2) / or end of file.
C
C                                                    
C Output formatted file SRFSECT (format LIN):
C (1) / (a slash).
C (2) For each triangle from the input file TRGL whose side(s) have a
C     crossection with the parallelogram given by points A,B and C
C     (2.1), (2.2), and (2.3):
C (2.1) 'wf',/ (text 'wf' and a slash)
C (2.2) Coordinates of the endpoints of the crossection segment:
C (2.2.1) X1,X2,X3,/
C     X1,X2,X3... Coordinates of the first endpoint.
C     /...    A slash.
C (2.2.2) Y1,Y2,Y3,/
C     Y1,Y2,Y3... Coordinates of the second endpoint.
C     /...    A slash.
C (2.3) / (a slash).
C (3) / (a slash).
C General description of format LIN
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
C.......................................................................
C
C     External functions and subroutines:
      EXTERNAL LENGTH,RSEP1,RSEP3T,RSEP3R,ERROR,FORM1
      INTEGER  LENGTH
C
C     Filenames and parameters:
      CHARACTER*80 FSEP,FVRTX,FTRGL,FOUT
      INTEGER LU1,LU2,LU3,IUNDEF
      PARAMETER (LU1=1,LU2=2,LU3=3,IUNDEF=-999999)
C     Other variables:
      CHARACTER*255 TEXT
      REAL U1,U2,U3,V1,V2,V3,A1,A2,A3,B1,B2,B3,C1,C2,C3,
     *     BA1,BA2,BA3,CA1,CA2,CA3,AA,BB,CC,DD
      LOGICAL LCROS
      INTEGER NVRTX,NPLGN,OPLGN,NQ,I1,I
C     LCROS.. Says whether a crossection of the current triangle with
C             the parallelogram was already found.
C
C.......................................................................
C
C     Reading main input data:
      WRITE(*,'(A)') '+SRFSECT: Enter input filename: '
      FSEP=' '
      READ (*,*) FSEP
      IF(FSEP.EQ.' ') THEN
C       SRFSECT-01
        CALL ERROR('SRFSECT-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)') '+SRFSECT: Working...            '
C
C     Reading input and output filenames:
      CALL RSEP3T('VRTX'   ,FVRTX,'vrtx.out')
      CALL RSEP3T('TRGL'   ,FTRGL,'trgl.out')
      CALL RSEP3T('SRFSECT',FOUT,'srfsect.out' )
C
C     Reading input parameters for the intersection parallelogram:
      CALL RSEP3R('A1',A1,0.)
      CALL RSEP3R('A2',A2,0.)
      CALL RSEP3R('A3',A3,0.)
      CALL RSEP3R('B1',B1,0.)
      CALL RSEP3R('B2',B2,0.)
      CALL RSEP3R('B3',B3,0.)
      CALL RSEP3R('C1',C1,0.)
      CALL RSEP3R('C2',C2,0.)
      CALL RSEP3R('C3',C3,0.)
C     The three points must be mutually different:
      IF (((A1.EQ.B1).AND.(A2.EQ.B2).AND.(A3.EQ.B3)).OR.
     *    ((A1.EQ.C1).AND.(A2.EQ.C2).AND.(A3.EQ.C3)).OR.
     *    ((B1.EQ.C1).AND.(B2.EQ.C2).AND.(B3.EQ.C3))) THEN
C       SRFSECT-02
        CALL ERROR('SRFSECT-02: Wrong input parallelogram')
C       The points A,B and C must be mutually different.
      ENDIF
C
C     Reading vertices:
      NQ=3
      OPEN(LU1,FILE=FVRTX,STATUS='OLD')
      READ(LU1,*) (TEXT,I=1,20)
      NVRTX=0
   20 CONTINUE
        IF(NVRTX+NQ.GT.MRAM) THEN
C         SRFSECT-03
          CALL ERROR('SRFSECT-03: Too small array RAM')
        END IF
        TEXT='$'
        DO 21 I=NVRTX+1,NVRTX+NQ
          RAM(I)=0.
   21   CONTINUE
        READ(LU1,*,END=29) TEXT,(RAM(I),I=NVRTX+1,NVRTX+NQ)
        IF(TEXT.EQ.'$') THEN
          GO TO 29
        END IF
C       Number of storage locations in RAM used for the vertices
        NVRTX=NVRTX+NQ
      GO TO 20
   29 CONTINUE
      CLOSE(LU1)
C
C     Reading the triangles:
      DO 81 I=NVRTX+1,MRAM
        IRAM(I)=0
   81 CONTINUE
      OPEN(LU1,FILE=FTRGL)
      NPLGN=0
      OPLGN=NVRTX
   82 CONTINUE
        IF(OPLGN+NPLGN+3.GT.MRAM) THEN
C         SRFSECT-04
          CALL ERROR('SRFSECT-04: Too small array RAM')
        END IF
        IRAM(OPLGN+NPLGN+1)=IUNDEF
        IRAM(OPLGN+NPLGN+2)=IUNDEF
        IRAM(OPLGN+NPLGN+3)=IUNDEF
        READ(LU1,*,END=89) (IRAM(I),I=OPLGN+NPLGN+1,OPLGN+NPLGN+3)
        IF(IRAM(OPLGN+NPLGN+1).EQ.IUNDEF) THEN
C         End of triangles:
          GO TO 89
        END IF
        IF ((IRAM(OPLGN+NPLGN+2).EQ.IUNDEF).OR.
     *      (IRAM(OPLGN+NPLGN+3).EQ.IUNDEF)) THEN
C         SRFSECT-05
          WRITE(TEXT,'(A,I6)')
     *      'SRFSECT-05: Missing vertices in triangle',NPLGN/3+1
          CALL ERROR(TEXT(1:LENGTH(TEXT)))
        END IF
        IF ((IRAM(OPLGN+NPLGN+1).EQ.IRAM(OPLGN+NPLGN+2)).OR.
     *      (IRAM(OPLGN+NPLGN+1).EQ.IRAM(OPLGN+NPLGN+3)).OR.
     *      (IRAM(OPLGN+NPLGN+2).EQ.IRAM(OPLGN+NPLGN+3))) THEN
C         SRFSECT-06
          WRITE(TEXT,'(A,I6)')
     *      'SRFSECT-06: Equal vertices in triangle',NPLGN/3+1
          CALL ERROR(TEXT(1:LENGTH(TEXT)))
        END IF
        NPLGN=NPLGN+3
      GO TO 82
   89 CONTINUE
      CLOSE(LU1)
C
C     Calculating the crossections:
      OPEN(LU3,FILE=FOUT)
      WRITE(LU3,'(A)') '/'
C     Parametric equation of the parallelogram:
C     X = A + s*[B-A] + t*[C-A]
      BA1=B1-A1
      BA2=B2-A2
      BA3=B3-A3
      CA1=C1-A1
      CA2=C2-A2
      CA3=C3-A3
C     General equation of the plane with the parallelogram:
C     AA*X1 + BB*X2 + CC*X3 + DD = 0
      AA=BA2*CA3-BA3*CA2
      BB=BA3*CA1-BA1*CA3
      CC=BA1*CA2-BA2*CA1
      DD=-1.*(AA*A1+BB*A2+CC*A3)
C     Loop over the triangles:
      DO 200, I1=OPLGN,OPLGN+NPLGN-3,3
        LCROS=.FALSE.
C       First side of the triangle:
        U1=RAM((IRAM(I1+1)-1)*NQ+1)
        U2=RAM((IRAM(I1+1)-1)*NQ+2)
        U3=RAM((IRAM(I1+1)-1)*NQ+3)
        V1=RAM((IRAM(I1+2)-1)*NQ+1)
        V2=RAM((IRAM(I1+2)-1)*NQ+2)
        V3=RAM((IRAM(I1+2)-1)*NQ+3)
        CALL SRFCRO(U1,U2,U3,V1,V2,V3,
     *       A1,A2,A3,BA1,BA2,BA3,CA1,CA2,CA3,AA,BB,CC,DD,LCROS,LU3)
C       Second side of the triangle:
        U1=RAM((IRAM(I1+1)-1)*NQ+1)
        U2=RAM((IRAM(I1+1)-1)*NQ+2)
        U3=RAM((IRAM(I1+1)-1)*NQ+3)
        V1=RAM((IRAM(I1+3)-1)*NQ+1)
        V2=RAM((IRAM(I1+3)-1)*NQ+2)
        V3=RAM((IRAM(I1+3)-1)*NQ+3)
        CALL SRFCRO(U1,U2,U3,V1,V2,V3,
     *       A1,A2,A3,BA1,BA2,BA3,CA1,CA2,CA3,AA,BB,CC,DD,LCROS,LU3)
C       Third side of the triangle:
        U1=RAM((IRAM(I1+2)-1)*NQ+1)
        U2=RAM((IRAM(I1+2)-1)*NQ+2)
        U3=RAM((IRAM(I1+2)-1)*NQ+3)
        V1=RAM((IRAM(I1+3)-1)*NQ+1)
        V2=RAM((IRAM(I1+3)-1)*NQ+2)
        V3=RAM((IRAM(I1+3)-1)*NQ+3)
        CALL SRFCRO(U1,U2,U3,V1,V2,V3,
     *       A1,A2,A3,BA1,BA2,BA3,CA1,CA2,CA3,AA,BB,CC,DD,LCROS,LU3)
        IF (LCROS) THEN
          WRITE(LU3,'(A)') '/'
        ENDIF
  200 CONTINUE
      WRITE(LU3,'(A)') '/'
C
      WRITE(*,'(A)') '+SRFSECT: Done.                 '
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE SRFCRO(U1,U2,U3,V1,V2,V3,
     *     A1,A2,A3,BA1,BA2,BA3,CA1,CA2,CA3,AA,BB,CC,DD,LCROS,LU3)
      INTEGER LU3
      REAL U1,U2,U3,V1,V2,V3,
     *     A1,A2,A3,BA1,BA2,BA3,CA1,CA2,CA3,AA,BB,CC,DD
      LOGICAL LCROS
C
C Subroutine designed to determine whether the abscissa formed by points
C U and V intersects the parallelogram defined by points A, B and C.
C
C Input:
C     LU3 ... Number of logical unit connected to the output file.
C     U1,U2,U3,V1,V2,V3... Coordinates of the points U and V.
C     A1,A2,A3... Coordinates of the point A.
C     BA1,BA2,BA3,CA1,CA2,CA3... Coordinates of the vectors B-A and C-A.
C     AA,BB,CC,DD... Coefficients of the general equation of the plane
C              formed by points A, B and C.
C     LCROS... Identifies whether any crossection was already found.
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
C.......................................................................
C     Local storage locations:
      INTEGER J1,J2
      CHARACTER*(25) FORMAT
      REAL VU1,VU2,VU3,W,SS,TT,X1,X2,X3,X(3),A(3),BA(3),CA(3),AUX1,AUX2
C
C.......................................................................
C
C     Parametric equation of the side:
C     X = U + w*[V-U]
      VU1=V1-U1
      VU2=V2-U2
      VU3=V3-U3
C     Value of w for the crossection point:
      AUX1=(VU1*AA+VU2*BB+VU3*CC)
      IF (AUX1.EQ.0.) THEN
C       Vector VU is parallel to the parallelogram:
        AUX2=AA*X1+BB*X2+CC*X3+DD
        IF (AUX2.EQ.0.) THEN
C         Abscissa VU lies in the plane of the parallelogram.
C         Such abscissa is excluded in the current version of the code.
          W=-999.
        ELSE
C         Abscissa VU lies away of the plane of the parallelogram.
          W=-999.
        ENDIF
      ELSE
        W=-1.*(U1*AA+U2*BB+U3*CC+DD)/AUX1
      ENDIF
      IF ((W.LE.1.).AND.(W.GE.0.)) THEN
C       Crossection of the triangle side with the parallelogram plane:
        X1=U1+W*VU1
        X2=U2+W*VU2
        X3=U3+W*VU3
C       Computing values of parameters s and t for the parallelogram:
C       Selecting the two of the three parametric equations of the
C       parallelogram with at least one nonzero coefficient for s or t:
        J1=0
        J2=0
        IF ((BA1.NE.0.).OR.(CA1.NE.0.)) THEN
          J1=1
        ENDIF
        IF ((BA2.NE.0.).OR.(CA2.NE.0.)) THEN
          IF (J1.EQ.0) THEN
            J1=2
            IF ((BA3.NE.0.).OR.(CA3.NE.0.)) THEN
              J2=3
            ELSE
C             SRFSECT-07
              CALL ERROR('SRFSECT-07: Error in parametric equation')
            ENDIF
          ELSE
            J2=2
          ENDIF
        ELSE
          IF (((BA3.NE.0.).OR.(CA3.NE.0.)).AND.(J1.NE.0)) THEN
            J2=3
          ELSE
C           SRFSECT-08
            CALL ERROR('SRFSECT-08: Error in parametric equation')
          ENDIF
        ENDIF
        A(1)=A1
        A(2)=A2
        A(3)=A3
        BA(1)=BA1
        BA(2)=BA2
        BA(3)=BA3
        CA(1)=CA1
        CA(2)=CA2
        CA(3)=CA3
        X(1)=X1
        X(2)=X2
        X(3)=X3
C       Calculating parameters s and t:
        SS=(X(J1)*CA(J2)-X(J2)*CA(J1)-A(J1)*CA(J2)+A(J2)*CA(J1))/
     *     (BA(J1)*CA(J2)-BA(J2)*CA(J1))
        TT=(X(J1)*BA(J2)-X(J2)*BA(J1)-A(J1)*BA(J2)+A(J2)*BA(J1))/
     *     (CA(J1)*BA(J2)-CA(J2)*BA(J1))
C       Values of s and t for the parallelogram must be within [0,1]:
        IF ((SS.GE.0.).AND.(SS.LT.1.)) THEN
        IF ((TT.GE.0.).AND.(TT.LT.1.)) THEN
          IF (.NOT.LCROS) THEN
C           First crossection of this triangle:
            WRITE(LU3,'(2A)') '''wf''',' /'
            LCROS=.TRUE.
          ENDIF
          FORMAT='(3('
          CALL FORM1(AMIN1(X1,X2,X3),AMAX1(X1,X2,X3),FORMAT(4:11))
          FORMAT(11:12)='))'
          WRITE(LU3,FORMAT) X1,' ',X2,' ',X3,' /'
        ENDIF
        ENDIF
      ENDIF
      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
C
C=======================================================================
C