      PROGRAM PLOTRAY

C     Plots central rays of beams generated by BELLHOP
C     Basically, all this does is change the file format to
C       one useable by the xmgr/xvgr plot package

      INTEGER RAYFIL, PRTFIL
      PARAMETER ( RAYFIL = 5, PRTFIL = 6, MAXNR = 100000 )

      LOGICAL     SURBNC, BOTBNC, NEWSET, SURBNCOLD, BOTBNCOLD
      REAL*4      R( MAXNR ), Z( MAXNR ), RV( MAXNR ), ZV( MAXNR )
      CHARACTER   TITLE*80

C      OPEN ( FILE = 'RAYFIL', UNIT = RAYFIL, STATUS = 'OLD',
C     &       FORM = 'FORMATTED' )

      READ( RAYFIL, * ) TITLE
      READ( RAYFIL, * ) FREQ
      READ( RAYFIL, * ) NBEAMS
      READ( RAYFIL, * ) DEPTHT
      READ( RAYFIL, * ) DEPTHB

C     *** Write data ***

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

C    *** Read rays ***

      NEWSET = .FALSE.
      SURBNCOLD = .FALSE.
      BOTBNCOLD = .FALSE.

      DO I = 1, NBEAMS
         SURBNC = .FALSE.
         BOTBNC = .FALSE.

C        --- read in a ray trajectory
         READ( RAYFIL, *, END = 3000 ) NSTEPS
         DO ISTEP = 1, NSTEPS
            READ( RAYFIL, * ) RV( ISTEP ), ZV( ISTEP )
         ENDDO

         R( 1 ) = RV( 1 ) / 1000.0
         Z( 1 ) = ZV( 1 )

         DO ISTEP = 2, NSTEPS
            R( ISTEP ) = RV( ISTEP ) / 1000.0
            Z( ISTEP ) = ZV( ISTEP )

C           --- Check ray lies within plot window
C            IF ( R( ISTEP - 1 ) .LT. RMAX .AND.
C     &           Z( ISTEP - 1 ) .LT. ZMAX ) N = ISTEP

C           --- Make note of ray type
            IF      ( ZV( ISTEP ) .LE. DEPTHT ) THEN
               SURBNC = .TRUE.
            ELSE IF ( ZV( ISTEP ) .GE. DEPTHB ) THEN
               BOTBNC = .TRUE.
            ENDIF
         END DO

C        --- is this a new ray type?
         IF ( SURBNC .NEQV. SURBNCOLD ) NEWSET = .TRUE.
         IF ( BOTBNC .NEQV. BOTBNCOLD ) NEWSET = .TRUE.

C        --- store ray type
         SURBNCOLD = SURBNC
         BOTBNCOLD = BOTBNC

C        *** Plot ***

C        --- separate ray families by '&' to mark a new 'set'
         IF ( NEWSET .AND. I .NE. 1 ) THEN
            WRITE( PRTFIL, * ) '&'
            NEWSET = .FALSE.
         END IF

         DO ISTEP = 1, NSTEPS
            WRITE( PRTFIL, * ) R( ISTEP ), Z( ISTEP )
         END DO

         WRITE( PRTFIL, * ) R( NSTEPS ), Z( 1 )
         WRITE( PRTFIL, * ) R( 1      ), Z( 1 )

      END DO

 3000 CONTINUE
      WRITE( PRTFIL, 2000 ) NINT( Z( 1 ) )
 2000 FORMAT(  '@ subtitle "',
     &         'Sd =', I4, ' m"' )

      RMIN = 0.0
      RMAX = R( NSTEPS )
      WRITE( PRTFIL, 4000 ) INT( RMIN   ), INT( RMAX   ),
     &                      INT( DEPTHT ), INT( DEPTHB )
 4000 FORMAT( '@ world xmin ', I5, /,
     &        '@ world xmax ', I5, /,
     &        '@ world ymin ', I5, /,
     &        '@ world ymax ', I5 )

      STOP
      END
