C
C=======================================================================
C
      PROGRAM RPPLOT
C
C-----------------------------------------------------------------------
C
C Version: 5.20
C Date: 1998, November 10
C
C Program to create simple PostScript plot of the distribution of the
C rays and homogeneous triangles either on the normalized ray domain or
C on the reference surface. The program is able to plot any of following
C objects: basic rays, two-point rays, auxiliary rays, homogeneous
C triangles, receivers. Objects are colour-coded according to the
C ray history, colours may be chosen. Rays are symbol-coded according
C to the ray history, symbols and their heights may be chosen. The
C limits of displayed part of the ray domain (or reference surface)
C may also be given by input data.
C
C.......................................................................
C
C                                                    
C Main input data read from external interactive device (*):
C     The data consist of character string, read by list
C     directed (free format) input.  The string has thus to be
C     enclosed in apostrophes.
C 'INPUTFILE'
C      Name of file with input data.
C Default: 'INPUTFILE'='rpplot.dat'
C
C                                                  
C Structure of the input data file 'INPUTFILE':
C     The data consist of following lines of integers, reals and
C     character strings, read by list directed (free format) input:
C (1) IRBAS,IRTWO,IRAUX,ITHOM,ISANG,ISHP,ISUC
C    Integer quantities describing, what is to be plotted:
C     IRBAS ... 1 if basic rays are to be plotted.
C     IRTWO ... 1 if two-point rays are to be plotted.
C     IRAUX ... 1 if auxiliary rays are to be plotted.
C     ITHOM ... 1 if homogeneous triangles are to be plotted.
C     ISANG ... 1 if the distribution of rays and (or) triangles
C               on the normalized ray domain is to be plotted. Otherwise
C               the distribution of (end)points of the successful rays
C               on the receiver surface will be plotted.
C     ISHP  ... 0 if objects of all histories are to be plotted,
C               otherwise only the objects of history ISHP will
C                 be plotted.
C     ISUC  ... 0 if objects of both positive and negative histories
C                 are to be plotted in colors according to 'COLORS',
C               1 if objects of positive histories are to be plotted
C                 in colors, and objects of negative histories are to
C                 be ploted in black.
C Defaults: IRBAS=1,IRTWO=1,IRAUX=1,ITHOM=1,ISANG=1,ISHP=0,ISUC=0
C
C (2) PLIMIT
C     PLIMIT... Limits of the plot, either in normalized ray
C               parameters (for ISANG=1) or in receiver surface
C               coordinates.
C     PLIMIT(1)... Minimum of first parameter (coordinate).
C     PLIMIT(2)... Maximum of first parameter (coordinate).
C     PLIMIT(3)... Minimum of second parameter (coordinate).
C     PLIMIT(4)... Maximum of second parameter (coordinate).
C Defaults: PLIMIT(1)=0,PLIMIT(2)=1,PLIMIT(3)=0,PLIMIT(4)=1
C
C (3) HRBAS,HRTWO,HRAUX,HOR,VER,HTEXT,'SYMBOLS','COLORS'
C    Real quantities describing height of the symbols on the plot:
C     HRBAS...  Height of symbols for basic rays.
C     HRTWO...  Height of symbols for two-point rays.
C     HRAUX...  Height of symbols for auxiliary rays.
C     HOR...    Horizontal dimension of the plot.
C     VER...    Vertical dimension of the plot.
C     HTEXT...  Height of text.
C     'SYMBOLS'... Name of the file with numbers of symbols.
C               This file has on each line just one integer number,
C               which is the number of symbol to be used for rays
C               and triangles of the history equal to number of the
C               line. i.e. on the first line is a number of symbol
C               to be used for plotting of rays and triangles with
C               history 1 or -1.
C               Description of file SYMBOLS
C     'COLORS'..Name of the file with numbers of colors.
C               Description of file COLORS
C     'RPPLOT'..Name of the output PostScript file with the figure.
C               If 'RPPLOT'=' ' (default), files 'plot00.ps',
C               'plot01.ps', ... are generated.
C Defaults: HRBAS=0.2,HRTWO=0.2,HRAUX=0.2,HOR=16.,VER=16.,HTEXT=0.5,
C           'SYMBOLS'=' ','COLORS'=' ','RPPLOT'=' '
C
C (4) 'CRT-I','CRT-T',NWAVE
C    Names of input files:
C     'CRT-I'...Name of the file with the quantities at the
C               initial points of rays (see C.R.T.6.1).
C               Description of file CRT-I
C     'CRT-T'...Name of the file with the indices of rays
C               forming the homogeneous triangles in the
C               ray-parameter domain.
C               Description of file CRT-T
C     NWAVE ... Index of the wave to be plotted. File with initial
C               points may contain several waves, but only one wave may
C               be plotted by RPPLOT.
C Defaults: 'CRT-I'='r01i.out', 'CRT-T'='t01.out', NWAVE=1.
C
C (5) 'RECEIVERS',ISREC,ICREC,HREC
C    Quantities and filename used for plotting receivers. Used only
C    when ISANG=0.
C     'RECEIVERS'... Name of the file with receiver coordinates. This
C               file should be of the structure of CRT input file.
C               Description of file RECEIVERS
C     ISREC ... Index of symbol to be used for plotting the receivers.
C     ICREC ... Color to be used for plotting the receivers.
C     HREC  ... Height of the symbols used for plotting the receivers.
C Defaults: 'RECEIVERS'='rec.dat',ISREC=3,ICREC=1,HREC=0.3.
C
C Example of data RPPLOT: len-rp2.dat
C Example of data RPPLOT: len-rp4.dat
C Example of data RPPLOT: len-rp5.dat
C Example of data RPPLOT: len-rp6.dat
C Example of data RPPLOT: len-rp7.dat
C
C
C                                                 
C Input formatted file 'SYMBOLS':
C     The file contains in I-th line just one integer telling the index
C     of symbols to be used to plot the rays or triangles of the history
C     I or -I. Only MSYMB different symbols are available, all the rays
C     with value of ray history greater or equal (in absolut value)
C     MSYMB will be ploted with the same symbols.
C For MSYMB refer to the file rpplot.inc.
C Example of data SYMBOL: symbol.dat
C Example of data SYMBOL: symbol7.dat
C
C                                                  
C Input formatted file 'COLORS':
C     The file contains in I-th line just one integer telling the index
C     of colour to be used to plot the rays or triangles of the history
C     I or -I. Only MCOL different colors are available, all the rays
C     and triangles with value of ray history greater or equal (in
C     absolut value) MCOL will be ploted in the same colors.
C For MCOL refer to the file rpplot.inc.
C Example of data COLORS: color.dat
C
C
C Coded by Petr Bulant
C     bulant@seis.karlov.mff.cuni.cz
C
C ...........................
C Common block /RAM/ to store the information about triangles and rays:
      INCLUDE 'rpplot.inc'
C
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER IRAY1,IRAY2,IRAY3,IR0
      INTEGER I1,I2,IT1
      INTEGER I,J,K,L,M,N
      INTEGER IRBAS,IRTWO,IRAUX,ITHOM,ISANG,ISREC,ICREC,NWAVE
      REAL HREC,P1,P2
      INTEGER LU0,LU2,LU3,LU4,LU5,LU6
      LOGICAL LEND
      PARAMETER (LU0=10,LU2=2,LU3=3,LU4=4,LU5=5,LU6=6)
      CHARACTER*80 FILE0,FILE2,FILE3,FILE4,FILE5,FILE6,CH,FILEPS
C-----------------------------------------------------------------------
C
      LEND=.FALSE.
C     Reading in a name of the main input data file:
      FILE0='rpplot.dat'
      WRITE(*,'(A)')
     *  ' Enter the name of the main input data file: '
      READ(*,*) FILE0
      OPEN(LU0,FILE=FILE0,FORM='FORMATTED',STATUS='OLD')
C     Reading in quantities telling what is to be plotted:
      IRBAS =1
      IRTWO =1
      IRAUX =1
      ITHOM =1
      ISANG =1
      ISHP  =0
      ISUC  =0
      READ(LU0,*) IRBAS,IRTWO,IRAUX,ITHOM,ISANG,ISHP,ISUC
      LRBAS=.FALSE.
      LRTWO=.FALSE.
      LRAUX=.FALSE.
      LTHOM=.FALSE.
      LSANG=.FALSE.
      IF (IRBAS.EQ.1) LRBAS=.TRUE.
      IF (IRTWO.EQ.1) LRTWO=.TRUE.
      IF (IRAUX.EQ.1) LRAUX=.TRUE.
      IF (ITHOM.EQ.1) LTHOM=.TRUE.
      IF (ISANG.EQ.1) LSANG=.TRUE.
C
C     Reading in limits of the plot:
      PLIMIT(1)=0.
      PLIMIT(2)=1.
      PLIMIT(3)=0.
      PLIMIT(4)=1.
      READ(LU0,*) PLIMIT(1),PLIMIT(2),PLIMIT(3),PLIMIT(4)
C
C     Reading in heights of symbols and names of files:
      HRBAS=0.2
      HRTWO=0.2
      HRAUX=0.2
      HOR=16.
      VER=16.
      HTEXT=0.5
      FILE4=' '
      FILE5=' '
      FILEPS=' '
      READ(LU0,*) HRBAS,HRTWO,HRAUX,HOR,VER,HTEXT,FILE4,FILE5,FILEPS
C
C     Reading in filenames of the files with computed triangles and rays
C     and index of the wave:
      FILE2='r01i.out'
      FILE3='t01.out'
      NWAVE=1
      READ(LU0,*) FILE2,FILE3,NWAVE
      OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU3,FILE=FILE3,FORM='FORMATTED',STATUS='OLD')
C
C     Reading in quantities and filename used for plotting receivers:
      FILE6='rec.dat'
      ISREC=3
      ICREC=1
      HREC=0.3
      READ(LU0,*) FILE6,ISREC,ICREC,HREC
C
C
      IF (LTHOM) THEN
C       Reading in file with computed triangles,
C       sorting the rays in triangles:
        NT=0
        NRAMP=0
        IRAYMI=1
        IRAYMA=0
  10    CONTINUE
          READ(LU3,*,END=20) IRAY1,IRAY2,IRAY3
          IF (IRAY1.GT.IRAYMA) IRAYMA=IRAY1
          IF (IRAY2.GT.IRAYMA) IRAYMA=IRAY2
          IF (IRAY3.GT.IRAYMA) IRAYMA=IRAY3
          IF (IRAY1.LT.IRAY2) THEN
            I=IRAY1
            IRAY1=IRAY2
            IRAY2=I
          ENDIF
          IF (IRAY2.LT.IRAY3) THEN
            I=IRAY2
            IRAY2=IRAY3
            IRAY3=I
          ENDIF
          IF (IRAY1.LT.IRAY2) THEN
            I=IRAY1
            IRAY1=IRAY2
            IRAY2=I
          ENDIF
          NRAMP=NRAMP+1
          IF (NRAMP.GT.MRAM) CALL RPPLER
          IRAM(NRAMP)=IRAY1
          NRAMP=NRAMP+1
          IF (NRAMP.GT.MRAM) CALL RPPLER
          IRAM(NRAMP)=IRAY2
          NRAMP=NRAMP+1
          IF (NRAMP.GT.MRAM) CALL RPPLER
          IRAM(NRAMP)=IRAY3
          NT=NT+1
        GOTO 10
  20    CONTINUE
        CLOSE(LU3)
        NR=IRAYMA
C
C       Sorting the triangles:
        DO 40, I1=NRAMP-5,1,-3
          DO 30, I2=1,I1,3
            L=I2+3
            IF (IRAM(I2).GT.IRAM(L)) THEN
              J=I2+1
              K=I2+2
              M=I2+4
              N=I2+5
              I       =IRAM(I2)
              IRAM(I2)=IRAM(L)
              IRAM(L) =I
              I       =IRAM(J)
              IRAM(J) =IRAM(M)
              IRAM(M) =I
              I       =IRAM(K)
              IRAM(K) =IRAM(N)
              IRAM(N) =I
            ENDIF
  30      CONTINUE
  40    CONTINUE
C
C
C       Forming an auxiliary array with information about when can be
C       rays erased from the memory ("deleting array"):
        IF (NRAMP+NR.GT.MRAM) CALL RPPLER
        DO 50, I1=NRAMP+1,NRAMP+NR
          IRAM(I1)=I1-NRAMP
  50    CONTINUE
        NRAMP=NRAMP+NR
        ORAYE=-3*NT
        DO 60, I1=1,3*NT,3
          IRAM(IRAM(I1  )-ORAYE)=IRAM(I1)
          IRAM(IRAM(I1+1)-ORAYE)=IRAM(I1)
          IRAM(IRAM(I1+2)-ORAYE)=IRAM(I1)
  60    CONTINUE
      ELSE
        NR=0
        NT=0
        ORAYE=0
        NRAMP=0
        IRAYMI=0
      ENDIF
C
C
C     Forming an auxiliary array with information about addresses
C     of the ends of records for rays in array RAM ("addressing array"):
C     "Ray" IRAYMI-1:
      NRAMP=NRAMP+1
      IF (NRAMP.GT.MRAM) CALL RPPLER
      IRAM(NRAMP)=NRAMP+NR
C     All other rays:
      IF (NRAMP+NR.GT.MRAM) CALL RPPLER
      DO 70, I1=NRAMP+1,NRAMP+NR
        IRAM(I1)=0
  70  CONTINUE
      NRAMP=NRAMP+NR
      ORAYA=-3*NT-NR-1
C
C
C     Choosing symbols:
      DO 71, I1=1,MSYMB
        ISYMB(I1)=I1
  71  CONTINUE
C
      OPEN(LU4,FILE=FILE4,STATUS='OLD',ERR=74)
      DO 72, I1=1,MSYMB
        READ (LU4,*,END=74) ISYMB(I1)
 72   CONTINUE
 74   CONTINUE
      CLOSE (LU4)
C
C     Choosing colours:
      DO 75, I1=1,MCOL
        ICOL(I1)=I1
  75  CONTINUE
C
      OPEN(LU5,FILE=FILE5,STATUS='OLD',ERR=78)
      DO 76, I1=1,MCOL
        READ (LU5,*,END=78) ICOL(I1)
 76   CONTINUE
 78   CONTINUE
      CLOSE (LU5)
C
      P1DIF=PLIMIT(2)-PLIMIT(1)
      P2DIF=PLIMIT(4)-PLIMIT(3)
      DO=(HOR+VER)/13.
C
      IF(FILEPS.NE.' ') THEN
        CALL PLOTN(FILEPS,0)
      END IF
      CALL PLOTS(0,0,0)
      CALL RPSYMB(0,0.,0.,0.)
C
C     Contures:
      CALL NEWPEN(1)
      CALL PLOT(HOR*((PLIMIT(1)-PLIMIT(1))/P1DIF)+DO,
     *          VER*((PLIMIT(3)-PLIMIT(3))/P2DIF)+DO,3)
      CALL PLOT(HOR*((PLIMIT(2)-PLIMIT(1))/P1DIF)+DO,
     *          VER*((PLIMIT(3)-PLIMIT(3))/P2DIF)+DO,2)
      CALL PLOT(HOR*((PLIMIT(2)-PLIMIT(1))/P1DIF)+DO,
     *          VER*((PLIMIT(4)-PLIMIT(3))/P2DIF)+DO,2)
      CALL PLOT(HOR*((PLIMIT(1)-PLIMIT(1))/P1DIF)+DO,
     *          VER*((PLIMIT(4)-PLIMIT(3))/P2DIF)+DO,2)
      CALL PLOT(HOR*((PLIMIT(1)-PLIMIT(1))/P1DIF)+DO,
     *          VER*((PLIMIT(3)-PLIMIT(3))/P2DIF)+DO,2)
      CALL PLOT(HOR*((PLIMIT(1)-PLIMIT(1))/P1DIF)+DO,
     *          VER*((PLIMIT(3)-PLIMIT(3))/P2DIF)+DO,3)
C
C     Labels along the figure:
      IF (LSANG) THEN
        CALL SYMBOL(HOR*( 0.1                     )+DO,
     *              VER*( 0.                      )+DO-1.5*HTEXT,
     *              HTEXT,'FIRST SHOOTING PARAMETER',0.,24)
        CALL SYMBOL(HOR*(0.                       )+DO-0.5*HTEXT,
     *              VER*(0.1                      )+DO       ,
     *              HTEXT,'SECOND SHOOTING PARAMETER',90.,25)
      ELSE
        CALL SYMBOL(HOR*( 0.1                     )+DO,
     *              VER*( 0.                      )+DO-1.5*HTEXT,
     *              HTEXT,'FIRST SURFACE COORDINATE',0.,24)
        CALL SYMBOL(HOR*(0.                       )+DO-0.5*HTEXT,
     *              VER*(0.1                      )+DO       ,
     *              HTEXT,'SECOND SURFACE COORDINATE',90.,25)
      ENDIF
C
C     Plotting receivers:
      IF (.NOT.LSANG) THEN
        OPEN(LU6,FILE=FILE6,FORM='FORMATTED',STATUS='OLD',ERR=80)
        READ(LU6,*,END=80) CH
        CALL NEWPEN(ICREC)
  81    CONTINUE
          CH='$'
          READ(LU6,*,END=80) CH,P1,P2
          IF (CH.EQ.'$') GOTO 81
          IF ((P1.GE.PLIMIT(1)).AND.(P1.LE.PLIMIT(2)).AND.
     *        (P2.GE.PLIMIT(3)).AND.(P2.LE.PLIMIT(4)))
     *         CALL RPSYMB(ISREC,
     *              HOR*((P1-PLIMIT(1))/P1DIF)+DO,
     *              VER*((P2-PLIMIT(3))/P2DIF)+DO,HREC)
        GOTO 81
      ENDIF
C
C
  80  CLOSE(LU6)
      IR0=0
      IF (LTHOM) THEN
C       Loop for all the triangles:
        IT1=1
  100   CONTINUE
C
C         If necessary, reading in new rays:
          IF ((IRAM(IRAM(IT1)-ORAYA+1).EQ.0).AND.(.NOT.LEND)) THEN
            CALL CIREAD(LU2,IRAM(IT1),NWAVE,LEND)
          ENDIF
C
C         Empting the array RAM:
          IF ((MRAM-NRAMP).LT.(MRAM/10.)) CALL CIERAS
C
C         Plotting the triangle:
          CALL RPTPL(IT1)
C
C         Plotting the rays:
          DO 102, I1=IR0+1,IRAM(IT1)
            CALL RPRPL(I1)
 102      CONTINUE
          IR0=IRAM(IT1)
C
        IT1=IT1+3
        IF (IT1.LT.3*NT) GOTO 100
C       End of the loop for all the triangles.
      ENDIF
C     Plotting remaining rays:
  110 CONTINUE
      IF (.NOT.LEND) THEN
        IR0=IR0+1
        CALL CIREAD(LU2,IR0,NWAVE,LEND)
        CALL RPRPL(IR0)
        GOTO 110
      ENDIF
C
      CALL PLOT(0.,0.,999)
      STOP
      END
C
C
C=======================================================================
C
      SUBROUTINE CIREAD(LU2,IR1,NWAVE,LEND)
C
C----------------------------------------------------------------------
C Subroutine to read the unformatted output of program CRT and
C to write it into array (I)RAM.
C Reading the output files is completed by a simple invocation of
C subroutine AP00 of file 'ap.for'.
C
      INTEGER LU2,IR1,NWAVE
      LOGICAL LEND
C Input:
C     LU2 ... Number of logical unit corresponding to the file with
C             the quantities at the initial points of rays.
C     IR1 ... Index of the first ray of the actually processed
C             triangle.
C     NWAVE.. Index of the elementary wave to be plotted.
C Output:
C     LEND .. .TRUE. when the end of file with rays is reached,
C             othervise .FALSE..
C
C Coded by Petr Bulant
C
C ...........................
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C None of the storage locations of the common block are altered.
C ...........................
C Common block /RAM/ to store the information about triangles and rays:
      INCLUDE 'rpplot.inc'
C
C-----------------------------------------------------------------------
C     Loop for the points of rays:
   10 CONTINUE
        IF ((NRAMP+2*NQ).GT.MRAM) THEN
C         Freeing the memory:
          CALL CIERAS
          IF ((NRAMP+2*NQ).GT.MRAM) CALL RPPLER
        ENDIF
C       Reading the results of the complete ray tracing:
        CALL AP00(0,0,LU2)
        IF (IWAVE.LT.1) THEN
C         End of rays:
          CLOSE(LU2)
          LEND=.TRUE.
          RETURN
        ENDIF
        IF (NWAVE.NE.IWAVE)
C         Skipping this elementary wave:
     *    GOTO 10
        IF (IRAY.LT.IRAYMI) GOTO 10
        IF (IRAM(IRAY-ORAYE).NE.0) THEN
C         Writing the results of the complete ray tracing - recording
C         new initial point on a ray:
          IF (LSANG) THEN
C           Normalized ray parameters:
            RAM(NRAMP+1)=YI(26)
            RAM(NRAMP+2)=YI(27)
          ELSE
C           Reference surface coordinates:
            RAM(NRAMP+1)=YI(28)
            RAM(NRAMP+2)=YI(29)
          ENDIF
C         Index of the receiver:
          IRAM(NRAMP+3)=IREC
C         History:
          IF (ISHEET.EQ.0) ISHEET=1
          IRAM(NRAMP+4)=ISHEET
          NRAMP=NRAMP+NQ
        ENDIF
        IRAM(IRAY-ORAYA)=NRAMP
        IF (IRAY.GT.IR1) RETURN
      GOTO 10
C
      END
C
C
C=======================================================================
C
      SUBROUTINE CIERAS
C
C----------------------------------------------------------------------
C Subroutine for empting the array (I)RAM. All the parameters
C of all the rays, which will no more be used, are erased.
C
C No input.
C No output.
C
C Subroutines and external functions required:
C
C Coded by Petr Bulant
C
C ...........................
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C     IRAY .. Index of the ray being actually read in by CIREAD.
C             This procedure supposes, that any ray with higher
C             index than IRAY was not read in.
C None of the storage locations of the common block are altered.
C ...........................
C Common block /RAM/ to store the information about triangles and rays:
      INCLUDE 'rpplot.inc'
C.......................................................................
C     Auxiliary storage locations:
      INTEGER I1,I2,J1
      INTEGER IADDRP
C     I1 ... Controls main loop over rays.
C     I2 ... Controls the loop over parameters of ray I1.
C     J1 ... address of the last used record of array RAM.
C
C-----------------------------------------------------------------------
      J1=IRAM(IRAYMI-1-ORAYA)
      IADDRP=J1
C     Loop for the rays:
      DO 20, I1=IRAYMI,IRAY
        IF (IRAM(I1-ORAYE).GE.(IRAY-1)) THEN
C         This ray is not to be erased:
          DO 10, I2=IADDRP+1,IRAM(I1-ORAYA)
            J1=J1+1
            RAM(J1)=RAM(I2)
   10     CONTINUE
          IADDRP=IRAM(I1-ORAYA)
          IRAM(I1-ORAYA)=J1
        ELSE
C         This ray is to be erased:
          IRAM(I1-ORAYE)=0
          IADDRP=IRAM(I1-ORAYA)
          IRAM(I1-ORAYA)=J1
        ENDIF
   20 CONTINUE
      NRAMP=J1
      RETURN
      END
C=======================================================================
C
      SUBROUTINE RPTPL(IT1)
C
C----------------------------------------------------------------------
C Subroutine for plotting the triangle formed by the rays
C IRAM(IT1), IRAM(IT1+1), IRAM(IT1+2).
C
      INTEGER IT1
C Input:
C     IT1 ... The address of the index of the first ray of the triangle
C             to be plotted.
C No output.
C
C Coded by Petr Bulant
C
C ...........................
C Common block /RAM/ to store the information about triangles and rays:
      INCLUDE 'rpplot.inc'
C.......................................................................
C Auxiliary storage locations:
      INTEGER I1,ILIN,ISH
      REAL P1A,P2A,P1B,P2B,P1C,P2C,P1,P2,P1O,P2O
      REAL P1MI,P1MA,P2MI,P2MA,DP1,DP2
      LOGICAL LPL
C-----------------------------------------------------------------------
      IF ((IRAM(IRAM(IT1  )-ORAYA).EQ.0).OR.
     *    (IRAM(IRAM(IT1  )-ORAYA).EQ.0).OR.
     *    (IRAM(IRAM(IT1  )-ORAYA).EQ.0)) THEN
C       RPPLOT-01
        CALL ERROR('RPPLOT-01: Parameters of a ray not found in memory')
C       This error may be caused by K2P 
C       not equal to zero, then only two-point rays are stored in
C       output files of CRT.
C       It may occur also in case, when a file with points along rays is
C       specified instead of input file with initial points of rays.
C       This also happens when there is less than NWAVE elementary waves
C       in the file with initial points.
       ENDIF
C
      P1A=RAM(IRAM(IRAM(IT1  )-ORAYA-1)+1)
      P2A=RAM(IRAM(IRAM(IT1  )-ORAYA-1)+2)
      P1B=RAM(IRAM(IRAM(IT1+1)-ORAYA-1)+1)
      P2B=RAM(IRAM(IRAM(IT1+1)-ORAYA-1)+2)
      P1C=RAM(IRAM(IRAM(IT1+2)-ORAYA-1)+1)
      P2C=RAM(IRAM(IRAM(IT1+2)-ORAYA-1)+2)
      ISH=IRAM(IRAM(IRAM(IT1+2)-ORAYA-1)+4)
C
      IF ((ISHP.NE.0).AND.(ISHP.NE.ISH)) RETURN
C
      IF ((.NOT.LSANG).AND.(ISH.LE.0)) RETURN
C
      P1MI=AMIN1(P1A,P1B,P1C)
      P1MA=AMAX1(P1A,P1B,P1C)
      P2MI=AMIN1(P2A,P2B,P2C)
      P2MA=AMAX1(P2A,P2B,P2C)
C
      ILIN=50*INT(AMAX1(((P1MA-P1MI)/P1DIF),((P2MA-P2MI)/P2DIF),1.))
      IF ((P1MI.GE.PLIMIT(1)).AND.(P1MA.LE.PLIMIT(2)).AND.
     *    (P2MI.GE.PLIMIT(3)).AND.(P2MA.LE.PLIMIT(4))) ILIN=1
C
      P1=P1A
      P2=P2A
      IF ((P1.GE.PLIMIT(1)).AND.(P1.LE.PLIMIT(2)).AND.
     *    (P2.GE.PLIMIT(3)).AND.(P2.LE.PLIMIT(4))) THEN
        IF ((ISUC.EQ.1).AND.(ISH.LT.0)) THEN
          CALL NEWPEN(1)
        ELSE
          CALL NEWPEN(ICOL(MIN0(MCOL,IABS(ISH))))
        ENDIF
        CALL PLOT(HOR*((P1-PLIMIT(1))/P1DIF)     +DO,
     *            VER*((P2-PLIMIT(3))/P2DIF)     +DO,3)
        P1O=P1
        P2O=P2
        LPL=.TRUE.
      ELSE
        LPL=.FALSE.
      ENDIF
      DP1=(P1B-P1A)/ILIN
      DP2=(P2B-P2A)/ILIN
      DO 30, I1=1,ILIN
        P1=P1+DP1
        P2=P2+DP2
        IF ((P1.GE.PLIMIT(1)).AND.(P1.LE.PLIMIT(2)).AND.
     *      (P2.GE.PLIMIT(3)).AND.(P2.LE.PLIMIT(4))) THEN
          IF (LPL) THEN
            CALL PLOT(HOR*((P1-PLIMIT(1))/P1DIF)     +DO,
     *                VER*((P2-PLIMIT(3))/P2DIF)     +DO,2)
            P1O=P1
            P2O=P2
          ELSE
            IF ((ISUC.EQ.1).AND.(ISH.LT.0)) THEN
              CALL NEWPEN(1)
            ELSE
              CALL NEWPEN(ICOL(MIN0(MCOL,IABS(ISH))))
            ENDIF
            CALL PLOT(HOR*((P1-PLIMIT(1))/P1DIF)     +DO,
     *                VER*((P2-PLIMIT(3))/P2DIF)     +DO,3)
            P1O=P1
            P2O=P2
            LPL=.TRUE.
          ENDIF
        ELSE
          IF (LPL) THEN
            CALL PLOT(HOR*((P1O-PLIMIT(1))/P1DIF)     +DO,
     *                VER*((P2O-PLIMIT(3))/P2DIF)     +DO,3)
            LPL=.FALSE.
          ENDIF
        ENDIF
  30  CONTINUE
      DP1=(P1C-P1B)/ILIN
      DP2=(P2C-P2B)/ILIN
      DO 32, I1=1,ILIN
        P1=P1+DP1
        P2=P2+DP2
        IF ((P1.GE.PLIMIT(1)).AND.(P1.LE.PLIMIT(2)).AND.
     *      (P2.GE.PLIMIT(3)).AND.(P2.LE.PLIMIT(4))) THEN
          IF (LPL) THEN
            CALL PLOT(HOR*((P1-PLIMIT(1))/P1DIF)     +DO,
     *                VER*((P2-PLIMIT(3))/P2DIF)     +DO,2)
            P1O=P1
            P2O=P2
          ELSE
            IF ((ISUC.EQ.1).AND.(ISH.LT.0)) THEN
              CALL NEWPEN(1)
            ELSE
              CALL NEWPEN(ICOL(MIN0(MCOL,IABS(ISH))))
            ENDIF
            CALL PLOT(HOR*((P1-PLIMIT(1))/P1DIF)     +DO,
     *                VER*((P2-PLIMIT(3))/P2DIF)     +DO,3)
            P1O=P1
            P2O=P2
            LPL=.TRUE.
          ENDIF
        ELSE
          IF (LPL) THEN
            CALL PLOT(HOR*((P1O-PLIMIT(1))/P1DIF)     +DO,
     *                VER*((P2O-PLIMIT(3))/P2DIF)     +DO,3)
            LPL=.FALSE.
          ENDIF
        ENDIF
  32  CONTINUE
      DP1=(P1A-P1C)/ILIN
      DP2=(P2A-P2C)/ILIN
      DO 34, I1=1,ILIN
        P1=P1+DP1
        P2=P2+DP2
        IF ((P1.GE.PLIMIT(1)).AND.(P1.LE.PLIMIT(2)).AND.
     *      (P2.GE.PLIMIT(3)).AND.(P2.LE.PLIMIT(4))) THEN
          IF (LPL) THEN
            CALL PLOT(HOR*((P1-PLIMIT(1))/P1DIF)     +DO,
     *                VER*((P2-PLIMIT(3))/P2DIF)     +DO,2)
            P1O=P1
            P2O=P2
          ELSE
            IF ((ISUC.EQ.1).AND.(ISH.LT.0)) THEN
              CALL NEWPEN(1)
            ELSE
              CALL NEWPEN(ICOL(MIN0(MCOL,IABS(ISH))))
            ENDIF
            CALL PLOT(HOR*((P1-PLIMIT(1))/P1DIF)     +DO,
     *                VER*((P2-PLIMIT(3))/P2DIF)     +DO,3)
            P1O=P1
            P2O=P2
            LPL=.TRUE.
          ENDIF
        ELSE
          IF (LPL) THEN
            CALL PLOT(HOR*((P1O-PLIMIT(1))/P1DIF)     +DO,
     *                VER*((P2O-PLIMIT(3))/P2DIF)     +DO,3)
            LPL=.FALSE.
          ENDIF
        ENDIF
  34  CONTINUE
      RETURN
      END
C=======================================================================
C
      SUBROUTINE RPRPL(IR1)
C
C----------------------------------------------------------------------
C Subroutine for plotting the ray IR1.
C
      INTEGER IR1
C Input:
C     IR1 ... Index of the ray to be plotted.
C No output.
C
C Coded by Petr Bulant
C
C ...........................
C Common block /RAM/ to store the information about triangles and rays:
      INCLUDE 'rpplot.inc'
C.......................................................................
C Auxiliary storage locations:
      INTEGER ISH,IREC
      REAL P1,P2
C-----------------------------------------------------------------------
      IF (IRAM(IR1-ORAYA).EQ.0) THEN
C       RPPLOT-02
        CALL ERROR('RPPLOT-02: Parameters of a ray not found in memory')
C       This error may be caused by K2P 
C       not equal to zero, then only two-point rays are stored in
C       output files of CRT.
C       It may occur also in case, when a file with points along rays is
C       specified instead of input file with initial points of rays.
C       This also happens when there is less than NWAVE elementary waves
C       in the file with initial points.
      ENDIF
C
      P1=RAM(IRAM(IR1-ORAYA-1)+1)
      P2=RAM(IRAM(IR1-ORAYA-1)+2)
      IREC=IRAM(IRAM(IR1-ORAYA-1)+3)
      ISH =IRAM(IRAM(IR1-ORAYA-1)+4)
C
      IF ((ISHP.NE.0).AND.(ISHP.NE.ISH)) RETURN
C
      IF ((.NOT.LSANG).AND.(ISH.LE.0)) RETURN
C
      IF ((P1.GE.PLIMIT(1)).AND.(P1.LE.PLIMIT(2)).AND.
     *    (P2.GE.PLIMIT(3)).AND.(P2.LE.PLIMIT(4))) THEN
        IF ((ISUC.EQ.1).AND.(ISH.LT.0)) THEN
          CALL NEWPEN(1)
        ELSE
          CALL NEWPEN(ICOL(MIN0(MCOL,IABS(ISH))))
        ENDIF
        IF (LRBAS.AND.(IREC.EQ.0)) THEN
C         Basic ray:
          CALL RPSYMB(ISYMB(MIN0(MSYMB,IABS(ISH))),
     *              HOR*((P1-PLIMIT(1))/P1DIF)+DO,
     *              VER*((P2-PLIMIT(3))/P2DIF)+DO,HRBAS)
          RETURN
        ENDIF
        IF (LRTWO.AND.(IREC.GT.0)) THEN
C         Two-point ray:
          CALL RPSYMB(ISYMB(MIN0(MSYMB,IABS(ISH))),
     *              HOR*((P1-PLIMIT(1))/P1DIF)+DO,
     *              VER*((P2-PLIMIT(3))/P2DIF)+DO,HRTWO)
          RETURN
        ENDIF
        IF (LRAUX.AND.(IREC.EQ.-1)) THEN
C         Auxiliary ray:
          CALL RPSYMB(ISYMB(MIN0(MSYMB,IABS(ISH))),
     *              HOR*((P1-PLIMIT(1))/P1DIF)+DO,
     *              VER*((P2-PLIMIT(3))/P2DIF)+DO,HRAUX)
          RETURN
        ENDIF
      ENDIF
      END
C=======================================================================
C
      SUBROUTINE RPPLER
C
C----------------------------------------------------------------------
C     RPPLOT-03
      CALL ERROR('RPPLOT-03: Array (I)RAM too small')
C     This error may be caused by too small dimension of array
C     RAM. Try to enlarge the parameter MRAM in common block
C     RAM.
      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
      INCLUDE 'calcops.for'
C     calcops.for
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'rpsymb.for'
C     rpsymb.for
C
C=======================================================================
C