      PROGRAM PLTTRI

C     Plot the triangulation

      INTEGER PRTFIL, FLPFIL
      PARAMETER ( PRTFIL = 6, FLPFIL = 5,
     &   MAXNSD = 21, MAXNRD = 51, MAXNT = 101,
     &   MXNODE = 10000, PI = 3.141592, DEGRAD = PI / 180.0 )

      REAL THETA( MAXNT ), SD( MAXNSD ), RD( MAXNRD ),
     &   X( MXNODE ), Y( MXNODE ), XV( 4 ), YV( 4 )
      CHARACTER OPT*50, TITLE*80, FILNAM( MXNODE )*50

C     *** Read FLPFIL header information ***

      READ( FLPFIL, * ) TITLE
      READ( FLPFIL, * ) OPT
      READ( FLPFIL, * ) M

C     --- Read source/receiver coordinates
      READ( FLPFIL, * ) XS, YS

C     --- Read source/rcvr depths
      CALL SDRD( FLPFIL, PRTFIL, 0.0, 1.0E6,
     &   MAXNSD, NSD, SD,
     &   MAXNRD, NRD, RD )

C     --- Radii of receivers
      READ( FLPFIL, * ) RMIN, RMAX, NR

C     --- Read angles for radials
      THETA( 3 ) = -999.9
      READ( FLPFIL, * ) NTHETA, ( THETA( ITH ), ITH = 1, NTHETA )
      CALL SUBTAB( THETA, NTHETA )
      IF ( NRD .GT. 1 ) NTHETA = 1

C     *** Nodal coordinates ***

      XMIN =  1.0E9
      XMAX = -1.0E9
      YMIN =  1.0E9
      YMAX = -1.0E9

      READ( FLPFIL, * ) NNODES
C     WRITE( *, * ) ' Number of Nodes', NNODES

      DO I = 1, NNODES
         READ( FLPFIL, * ) X( I ), Y( I ), FILNAM( I )
C        WRITE( *, * ) X(I), Y(I), FILNAM(I)
         X(I) = X( I )
         Y(I) = Y( I )
         XMIN = MIN( X( I ), XMIN )
         XMAX = MAX( X( I ), XMAX )
         YMIN = MIN( Y( I ), YMIN )
         YMAX = MAX( Y( I ), YMAX )
      END DO

      LAST = INDEX( TITLE, '    ' ) - 1
      WRITE( PRTFIL, 1000 ) TITLE(1:LAST)
 1000 FORMAT( '@ title "', A, '"' )

      WRITE( PRTFIL, 2000 ) INT( XMIN ), INT( XMAX ),
     &                      INT( YMIN ), INT( YMAX )
 2000 FORMAT( '@ world xmin ', I5, /,
     &        '@ world xmax ', I5, /,
     &        '@ world ymin ', I5, /,
     &        '@ world ymax ', I5 )

C     ------ Element definition
      READ( FLPFIL, * ) NELTS
C      WRITE( *, * ) 'NELTS = ', NELTS

      DO IELT = 1, NELTS
         READ( FLPFIL, * ) NODE1, NODE2, NODE3
C        WRITE( *, * ) 'NODE1, NODE2, NODE3 = ', NODE1, NODE2, NODE3
         XV( 1 ) = X( NODE1 )
         YV( 1 ) = Y( NODE1 )
         XV( 2 ) = X( NODE2 )
         YV( 2 ) = Y( NODE2 )
         XV( 3 ) = X( NODE3 )
         YV( 3 ) = Y( NODE3 )
         XV( 4 ) = X( NODE1 )
         YV( 4 ) = Y( NODE1 )

         DO INODE = 1, 4
            WRITE( *, * ) XV( INODE ), YV( INODE )
         END DO

C         IF ( FILNAM( NODE1 ) .EQ. FILNAM( NODE2 ) ) THEN
C            CALL SETCLR( 'BLACK' )
C         ELSE
C            CALL SETCLR( 'RED' )
C         ENDIF
C         CALL CURVE( XV( 1 ), YV( 1 ), 2, 1 )

C         IF ( FILNAM( NODE2 ) .EQ. FILNAM( NODE3 ) ) THEN
C            CALL SETCLR( 'BLACK' )
C         ELSE
C            CALL SETCLR( 'RED' )
C         ENDIF
C         CALL CURVE( XV( 2 ), YV( 2 ), 2, 1 )

C         IF ( FILNAM( NODE3 ) .EQ. FILNAM( NODE1 ) ) THEN
C            CALL SETCLR( 'BLACK' )
C         ELSE
C            CALL SETCLR( 'RED' )
C         ENDIF
C         CALL CURVE( XV( 3 ), YV( 3 ), 2, 1 )

C        ------ Label the element
         XMID = ( X( NODE1 ) + X( NODE2 ) + X( NODE3 ) ) / 3.0
         YMID = ( Y( NODE1 ) + Y( NODE2 ) + Y( NODE3 ) ) / 3.0
C         CALL RLMESS( '.$', 100, XMID, YMID )
C         CALL RLINT( IELT, XMID, YMID )
      END DO   ! next element

      WRITE( *, * ) '&'

C     Label the nodes

C      DO INODE = 1, NNODES
C         CALL RLINT( INODE, X( INODE ), Y( INODE ) )
C      END DO

C     Plot the bearing lines

      DO ITH = 1, NTHETA
         XV( 1 ) = XS + RMIN * COS( DEGRAD * THETA( ITH ) )
         YV( 1 ) = YS + RMIN * SIN( DEGRAD * THETA( ITH ) )

         XV( 2 ) = XS + RMAX * COS( DEGRAD * THETA( ITH ) )
         YV( 2 ) = YS + RMAX * SIN( DEGRAD * THETA( ITH ) )

         WRITE( *, * ) XV( 1 ), YV( 1)
         WRITE( *, * ) XV( 2 ), YV( 2)
      END DO
      WRITE( *, * ) '&'

      STOP
      END
