C 
C Program RPPLOT to plot ray parameters
C
C Version: 5.40
C Date: 2000, May 15
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 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 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 Integers describing, what is to be plotted:
C     IRBAS=integer ... 1 if basic rays are to be plotted. Other value
C             means that the basic rays are not to be plotted.
C             Default: IRBAS=1
C     IRTWO=integer ... 1 if two-point rays are to be plotted.
C             Other value disables plotting of the two-point rays.
C             Default: IRTWO=1
C     IRAUX=integer ... 1 if auxiliary rays are to be plotted.
C             Other value disables plotting of the auxiliary rays.
C             Default: IRAUX=1
C     ITHOM=integer ... 1 if homogeneous triangles are to be plotted.
C             Other value disables plotting of the triangles.
C             Default: ITHOM=1
C     ISANG=integer ... 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             Default: ISANG=1
C     ISHP=integer ... 0 if objects of all histories are to be plotted,
C             otherwise only the objects of the history equal to ISHP
C             will be plotted.
C             Default: ISHP=0
C     ISUC=integer ... 1 if objects of positive histories are to be
C             plotted in colors according to COLORS, and objects of
C             negative histories are to be ploted in black. Note that
C             positive ray histories correspond to the successful rays,
C             which pass the reference surface.
C             Otherwise objects of all histories are to be plotted
C             in colors according to file COLORS.
C             Default: ISUC=0
C
C Limits of the plot, either in normalized ray parameters (if ISANG=1),
C or in reference surface coordinates.
C     PLIM1=real ... Minimum of first parameter (coordinate).
C             Default: PLIM1=0.
C     PLIM2=real ... Maximum of first parameter (coordinate).
C             Default: PLIM2=1.
C     PLIM3=real ... Minimum of second parameter (coordinate).
C             Default: PLIM3=0.
C     PLIM4=real ... Maximum of second parameter (coordinate).
C             Default: PLIM4=1.
C
C Heights of symbols on the resulting plot in centimeters:
C     HRBAS=real ... Height of symbols for basic rays.
C             Default: HRBAS=0.2
C     HRTWO=real ... Height of symbols for two-point rays.
C             Default: HRTWO=0.2
C     HRAUX=real ... Height of symbols for auxiliary rays.
C             Default: HRAUX=0.2
C     HTEXT=real ... Height of the text along axes of the figure.
C             Default: HTEXT=0.5
C
C Dimensiones of the plot in centimeters:
C     HSIZE=real ... Horizontal dimension of the plot.
C             Default: HSIZE=15.
C     VSIZE=real ... Vertical dimension of the plot.
C             Default: VSIZE=15.
C
C Names of output and input files:
C     RPPLOT='string' ... Name of the output PostScript file with the
C             figure. If RPPLOT=' ' (default), files 'plot00.ps',
C             'plot01.ps', ... are generated.
C             Default: RPPLOT=' '
C     SYMBOLS='string' ... Name of the file with numbers of symbols.
C             This file has on each line a single 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             Default: SYMBOLS=' '
C     COLORS='string' ... Name of the file with numbers of colors.
C             Description of file COLORS.
C             Default: COLORS=' '
C     CRTOUT='string'...File with the names of the output files of
C             program CRT. If blank, default names are considered.
C             Only the files 'CRT-I' and 'CRT-T' are read by RPPLOT.
C             For general description of file CRTOUT refer to
C             file writ.for.
C             Default: CRTOUT=' ' which means 'CRT-I'='s01i.out' and
C               'CRT-T'='t01.out'
C     NWAVE=positive integer ... Index of the wave to be plotted. File
C             'CRT-I' with initial points may contain several waves, but
C             only one wave may be plotted by a single invocation
C             of RPPLOT.
C             Default: NWAVE=1
C
C Filename and parameters used for plotting receivers. Used only
C when ISANG=0.
C     REC='string'... If non-blank, the name of the file with the names
C             and coordinates of the receiver points.
C             Description of file REC.
C             Default: REC=' '
C     ICREC=positive integer ... Color to be used for plotting the
C             receivers. The color is chosen according to COLORS.
C             Default: ICREC=1
C     ISREC=positive integer ... Index of symbol to be used for plotting
C             the receivers according to the file SYMBOLS.
C             Default: ISREC=3
C     HREC=real ...  Height (in centimeters) of the symbols used for
C             plotting the receivers.
C             Default: HREC=0.3
C
C Example of the input parameters
C
C
C                                                 
C Input formatted file SYMBOLS:
C     The file contains in I-th line a single integer telling the index
C     of symbol 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 SYMBOLS. 
C Another example of data SYMBOLS.
C
C                                                  
C Input formatted file COLORS:
C     The file contains in I-th line a single integer telling the index
C     of colour to be used to plot the rays or triangles of the history
C     I or -I according to the file
C     calcops.rgb.
c     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. 
C
C ......................................................................
C Subroutines and external functions required:
      EXTERNAL RPPLER,CIREAD,CIERAS,RPTPL,RPRPL,RPSYMB,SYMBOL,
     *ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,PLOTS,PLOT,NEWPEN,PLOTN,AP00
C     RPPLER,CIREAD,CIERAS,RPTPL,RPRPL ... This file.
C     RPSYMB ... File rpsymb.for.
C     ERROR ... File error.for.
C     RSEP1,RSEP3I,RSEP3T,RSEP3R ...
C                           File sep.for.
C     PLOTS,PLOT,NEWPEN,PLOTN,SYMBOL ...
C                   File calcops.for.
C     AP00 ... File ap.for.
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,LU1,LU2,LU3,LU4,LU5,LU6
      LOGICAL LEND
      PARAMETER (LU0=10,LU1=1,LU2=2,LU3=3,LU4=4,LU5=5,LU6=6)
      CHARACTER*80 FILSEP,FILINI,FILTRI,FILSYM,FILCOL,FILCRT,FILREC,
     *             CH,CH1,FILEPS
C-----------------------------------------------------------------------
C
      LEND=.FALSE.
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+RPPLOT: Enter input filename: '
      FILSEP=' '
      READ(*,*) FILSEP
C
C     Reading all data from the SEP file into the memory:
      IF (FILSEP.NE.' ') THEN
        CALL RSEP1(LU0,FILSEP)
      ELSE
C       RPPLOT-04
        CALL ERROR('RPPLOT-04: 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

      WRITE(*,'(A)') '+RPPLOT: Working...            '
C
C     Reading quantities telling what is to be plotted:
      CALL RSEP3I('IRBAS',IRBAS,1)
      CALL RSEP3I('IRTWO',IRTWO,1)
      CALL RSEP3I('IRAUX',IRAUX,1)
      CALL RSEP3I('ITHOM',ITHOM,1)
      CALL RSEP3I('ISANG',ISANG,1)
      CALL RSEP3I('ISHP ',ISHP ,0)
      CALL RSEP3I('ISUC ',ISUC ,0)
      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 limits of the plot:
      CALL RSEP3R('PLIM1',PLIMIT(1),0.)
      CALL RSEP3R('PLIM2',PLIMIT(2),1.)
      CALL RSEP3R('PLIM3',PLIMIT(3),0.)
      CALL RSEP3R('PLIM4',PLIMIT(4),1.)
C
C     Reading heights of symbols and names of files:
      CALL RSEP3R('HRBAS',HRBAS,0.2)
      CALL RSEP3R('HRTWO',HRTWO,0.2)
      CALL RSEP3R('HRAUX',HRAUX,0.2)
      CALL RSEP3R('HSIZE',HOR,15.)
      CALL RSEP3R('VSIZE',VER,15.)
      CALL RSEP3R('HTEXT',HTEXT,0.5)
      CALL RSEP3T('SYMBOLS',FILSYM,' ')
      CALL RSEP3T('COLORS',FILCOL,' ')
      CALL RSEP3T('RPPLOT',FILEPS,' ')
C
C     Reading filenames of the files with computed triangles and rays
C     and index of the wave:
      CALL RSEP3T('CRTOUT',FILCRT,' ')
      FILINI='s01i.out'
      FILTRI='t01.out'
      IF (FILCRT.NE.' ') THEN
        OPEN(LU1,FILE=FILCRT,FORM='FORMATTED',STATUS='OLD')
        READ(LU1,*) CH,CH1,FILINI,FILTRI
        CLOSE(LU1)
      ENDIF
      CALL RSEP3I('NWAVE',NWAVE,1)
C
C     Reading quantities and filename used for plotting receivers:
      CALL RSEP3T('REC',FILREC,' ')
      CALL RSEP3I('ISREC',ISREC,3)
      CALL RSEP3I('ICREC',ICREC,1)
      CALL RSEP3R('HREC',HREC,0.3)
C
C
      IF (LTHOM) THEN
C       Reading file with computed triangles,
C       sorting the rays in triangles:
        NT=0
        NRAMP=0
        IRAYMI=1
        IRAYMA=0
        OPEN(LU3,FILE=FILTRI,FORM='FORMATTED',STATUS='OLD')
  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=FILSYM,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=FILCOL,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=FILREC,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)
      OPEN(LU2,FILE=FILINI,FORM='UNFORMATTED',STATUS='OLD')
      IR0=0
      IF (LTHOM) THEN
C       Loop for all the triangles:
        IT1=1
  100   CONTINUE
C
C         If necessary, reading 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)
      WRITE(*,'(A)') '+RPPLOT: Done.                 '
      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