C 
C Program PFATUBES to modify the file 'CRT-T' with indices of triangles
C in order to enable interpolation within ray tubes formed by rays
C written by program CRTPFA.
C
C Version: 7.00
C Date: 2013, August 14
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
C     SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Names of input files related to ray tracing (output of CRT):
C     CRTOUT='string'... File with the names of the output files of
C             program CRT.
C             For general description of file CRTOUT refer to file
C             writ.for.
C             Description specific to this program:
C               Just the first set of CRT output files is read from
C               file CRTOUT.  Only files
C               'CRT-R',
C               'CRT-I' and
C               'CRT-T'
C               are read.
C               File 'CRT-R' must contain all rays traced by CRT, not
C               only two-point rays.
C             Default: CRTOUT=' ' which means 'CRT-R'='r01.out',
C             'CRT-I'='r01i.out', 'CRT-T'='t01.out'
C Names of output files:
C     FTRGL1... Name of the output file with the triangles formed
C             by rays which belong each to other and thus enable
C             the interpolation within ray tubes.
C             No default, FTRGL1 must be specified
C     FTRGL2... Name of the output file with the triangles formed
C             by rays which do not belong each to other and do not
C             enable the interpolation within ray tubes.
C             No default, FTRGL2 must be specified
C Optional parameters specifying the form of the quantities
C written in the output formatted files:
C     MINDIG,MAXDIG=positive integers... See the description in file
C             forms.for.
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL LOWER,FORM2,LENGTH,ERROR,WARN,RSEP1,RSEP3T,
     *INDEXI,AP00,PTREAD,PTERAS,PTINDX,PTQW,PTPOLA
      INTEGER LENGTH
C     PTREAD,PTERAS,PTINDX,PTQW,PTPOLA... This file.
C     LOWER,LENGTH... File length.for.
C     ERROR,WARN... File error.for.
C     RSEP1,RSEP3T... File sep.for.
C     INDEXI... File indexi.for.
C     AP00... File ap.for.
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 /PTCB/ to store the information about triangles and rays:
      INCLUDE 'pfatubes.inc'
C pfatubes.inc.
C
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER IRAY1,IRAY2,IRAY3
      INTEGER I1,I2,I3,I4,IT1,IIT1,JWAVES
      INTEGER I
      INTEGER LU0,LU1,LU2,LU3,LU4,LU5
      PARAMETER (LU0=10,LU1=1,LU2=2,LU3=3,LU4=4,LU5=5)
      REAL AUX
      CHARACTER*80 FILSEP,FILCRT,FILE1,FILE2,FILE3,CH,FTRGL1,FTRGL2
      CHARACTER*52 TXTERR
      LOGICAL LENDT
C
C-----------------------------------------------------------------------
C
C     Reading a name of the file with the input data:
      FILSEP=' '
      WRITE(*,'(A)') '+PFATUBES: Enter input filename: '
      READ(*,*) FILSEP
      IF (FILSEP.EQ.' ') THEN
C       PFATUBES-01
        CALL ERROR('PFATUBES-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.
      ENDIF
      WRITE(*,'(A)') '+PFATUBES: Reading input data.   '
C
C     Reading all the data from the SEP file into the memory:
      CALL RSEP1(LU1,FILSEP)
C
C     Reading filenames of the files with computed rays
C     and triangles:
      CALL RSEP3T('CRTOUT',FILCRT,' ')
      FILE1='r01.out'
      FILE2='r01i.out'
      FILE3='t01.out'
      IF (FILCRT.NE.' ') THEN
        OPEN (LU2,FILE=FILCRT,FORM='FORMATTED',STATUS='OLD')
        READ (LU2,*) FILE1,CH,FILE2,FILE3
        CLOSE (LU2)
      ENDIF
C
C     Opening the CRT output files:
      OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU3,FILE=FILE3,FORM='FORMATTED',STATUS='OLD')
      IRAY=0
      IWAVE=0
C
C     Reading filenames of the output files:
      CALL RSEP3T('FTRGL1',FTRGL1,' ')
      CALL RSEP3T('FTRGL2',FTRGL2,' ')
      IF ((FTRGL1.EQ.' ').OR.(FTRGL2.EQ.' ')) THEN
C       PFATUBES-02
        CALL ERROR('PFATUBES-02: No output file specified')
C       Output files FTRGL1 and FTRGL2 must be specified,
C       there is no default filename.
      ENDIF
C     Opening the output files:
      OPEN(LU4,FILE=FTRGL1)
      OPEN(LU5,FILE=FTRGL2)
C
C     Reading file with computed triangles,
C     sorting the rays in triangles:
      WRITE(*,'(A)') '+PFATUBES: Reading the file with triangles.'
      MRAMP=MRAM
      IFORM=99999
      FORMAT='(I6,I6,I6)'
      JWAVES=0
      LENDT=.FALSE.
      ISRCN =0
      IWAVEN=0
      READ(LU3,*,END=20) IRAY1,IRAY2,IRAY3
      WRITE(LU4,FORMAT) IRAY1,IRAY2,IRAY3
      WRITE(LU5,FORMAT) IRAY1,IRAY2,IRAY3
   8  CONTINUE
      NT=0
      NRAMP=0
      IRAYMI=999999
      IRAYMA=0
  10  CONTINUE
        IRAY1=0
        IRAY2=0
        IRAY3=0
        READ(LU3,*,END=20) IRAY1,IRAY2,IRAY3
        IF (IRAY1.EQ.0) THEN
C         New elementary wave:
          ISRCN =IRAY2
          IWAVEN=IRAY3
          GOTO 22
        ENDIF
        IF ((IRAY1.LE.0).OR.(IRAY2.LE.0).OR.(IRAY3.LT.0)) THEN
C         PFATUBES-03
          CALL ERROR('PFATUBES-03: Disorder in file ''CRT-T''')
C         The 'triangles' in the input file 'CRT-T' should be written
C         as three nonnegative integers (indices of rays), two of them
C         being positive (the third may equal zero in 2-D computations).
        ENDIF
        IF (IRAY1.GT.IRAYMA) IRAYMA=IRAY1
        IF (IRAY2.GT.IRAYMA) IRAYMA=IRAY2
        IF (IRAY3.GT.IRAYMA) IRAYMA=IRAY3
        IF (IRAY1.LT.IRAYMI)                    IRAYMI=IRAY1
        IF (IRAY2.LT.IRAYMI)                    IRAYMI=IRAY2
        IF ((IRAY3.NE.0).AND.(IRAY3.LT.IRAYMI)) IRAYMI=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
        IF (NRAMP+6.GT.MRAMP) THEN
C         PFATUBES-04
          CALL ERROR('PFATUBES-04: Small array RAM')
        ENDIF
        NRAMP=NRAMP+1
        IRAM(NRAMP)=IRAY1
        NRAMP=NRAMP+1
        IRAM(NRAMP)=IRAY2
        NRAMP=NRAMP+1
        IRAM(NRAMP)=IRAY3
        NRAMP=NRAMP+3
        NT=NT+1
      GOTO 10
  20  CONTINUE
      LENDT=.TRUE.
      CLOSE(LU3)
  22  CONTINUE
      JWAVES=JWAVES+1
      IRAYMI=2*IRAYMI-1
      IRAYMA=2*IRAYMA
      NR=IRAYMA-IRAYMI+1
C
C     Check for the 'triangles' in case of 2D computation:
      IF (IRAM(3).EQ.0) THEN
C        DO 25, I1=9,6*NT,6
C          IF (IRAM(I1).NE.0) THEN
CC           PFATUBES-XX
C            CALL ERROR('PFATUBES-XX: Disorder in triangles.')
CC           The first 'triangle' is formed by two rays with zero
CC           instead of the third ray. This identifies 2-D calculation
CC           and all the 'triangles' should be formed by two rays and
CC           zero instead of the third ray. This is not the case of the
CC           current triangle. Check the input file 'CRT-T'.
C          ENDIF
C  25    CONTINUE
C        L3D=.FALSE.
C        L2D=.TRUE.
C       PFATUBES-05
        CALL ERROR('PFATUBES-05: 2-D calculation')
C       The first 'triangle' is formed by two rays with zero instead
C       of the third ray. This identifies 2-D calculation, which
C       is not yet coded in this program.
      ELSE
        L3D=.TRUE.
        L2D=.FALSE.
      ENDIF
C
C     Sorting the triangles:
      WRITE(*,'(A)') '+PFATUBES: Sorting the triangles.          '
      IF (NRAMP+2*NT.GT.MRAMP) THEN
C         PFATUBES-06
          CALL ERROR('PFATUBES-06: Small array RAM')
      ENDIF
      DO 30, I1=NRAMP+1,NRAMP+NT
        IRAM(I1)=IRAM((I1-NRAMP-1)*6+1)
  30  CONTINUE
      CALL INDEXI(NT,IRAM(NRAMP+1),IRAM(NRAMP+NT+1))
      DO 40, I1=NRAMP+1,NRAMP+NT
        IRAM(I1)=IRAM(I1+NT)
  40  CONTINUE
      NRAMP=NRAMP+NT
C
C
C     Forming an auxiliary array with information about when can be rays
C     erased from the memory ("deleting array"):
      IF (NRAMP+NR.GT.MRAMP) THEN
C       PFATUBES-07
        CALL ERROR('PFATUBES-07: Small array RAM')
      ENDIF
C     Marking all the rays as "not to be used":
      DO 50, I1=NRAMP+1,NRAMP+NR
        IRAM(I1)=0
  50  CONTINUE
      NRAMP=NRAMP+NR
      ORAYE=IRAYMI-1-7*NT
C     Marking rays which will be used:
C     Loop over the sorted indices of triangles:
      DO 60, I2=6*NT+1,7*NT
C       Address of the first record of the current triangle:
        I1=(IRAM(I2)-1)*6+1
C       Marking all the six rays corresponding to the triangle:
        I3=2*IRAM(I1)-1
        I4=2*IRAM(I1)
        IRAM(I3-ORAYE)=I4
        IRAM(I3+1-ORAYE)=I4
        I3=2*IRAM(I1+1)-1
        IRAM(I3-ORAYE)=I4
        IRAM(I3+1-ORAYE)=I4
        IF (IRAM(I1+2).NE.0) THEN
          I3=2*IRAM(I1+2)-1
          IRAM(I3-ORAYE)=I4
          IRAM(I3+1-ORAYE)=I4
        ENDIF
  60  CONTINUE
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.MRAMP) THEN
C       PFATUBES-08
        CALL ERROR('PFATUBES-08: Small array RAM')
      ENDIF
      IRAM(NRAMP)=NRAMP+NR
C     All other rays:
      IF (NRAMP+NR.GT.MRAMP) THEN
C       PFATUBES-09
        CALL ERROR('PFATUBES-09: Small array RAM')
      ENDIF
      DO 70, I1=NRAMP+1,NRAMP+NR
        IRAM(I1)=0
  70  CONTINUE
      NRAMP=NRAMP+NR
      ORAYA=IRAYMI-1-7*NT-NR-1
C
C
C     Loop for all the triangles:
      IIT1=1
      I1=-1
  100 CONTINUE
C       Screen output:
        IT1=(IRAM(6*NT+IIT1)-1)*6+1
        I2=NINT((100.*IIT1)/(NT))
        IF (I2.NE.I1) THEN
          WRITE(*,'(''+'',79('' ''))')
          WRITE(*,'(A,1I2,A,1I4,A)')
     *    '+PFATUBES: Interpolating wave ',JWAVES,'. ',
     *                                    I2,'% of ray tubes completed.'
          I1=I2
        ENDIF
C
C       If necessary, reading new rays:
        IF
     *  (((IRAM(2*IRAM(IT1)-ORAYA+1).EQ.0).AND.(2*IRAM(IT1).NE.IRAYMA))
     *   .OR.
     *   ((IRAM(2*IRAM(IT1)-ORAYA  ).EQ.0).AND.(2*IRAM(IT1).EQ.IRAYMA)))
     *  CALL PTREAD(LU1,LU2,IT1)
C
C       Empting the array RAM to enable writing of possible interpolated
C       quantities:
        IF ((MRAMP-NRAMP).LT.(NINT(MAXRAM/10.))) CALL PTERAS
C
C       Calculating the indices of sides of the triangle:
        CALL PTINDX(IT1)
C
      IIT1=IIT1+1
      IF (IIT1.LE.NT) GOTO 100
C     End of the loop for all the triangles of the elementary wave.
C
C     Writing the resulting triangles:
      CALL PTQW(LU4,LU5)
C
      IF (.NOT.LENDT) THEN
C       Continuing with the next elementary wave:
        GOTO 8
      ENDIF
C
      CLOSE(LU1)
      CLOSE(LU2)
C
C
      WRITE(*,'(A)')
     *    '+PFATUBES: Done.                                            '
      STOP
      END
C
C
C=======================================================================
C
      SUBROUTINE PTQW(LU4,LU5)
C
C----------------------------------------------------------------------
C Subroutine to write the files with triangles of one elementary wave.
C
      INTEGER LU4,LU5
C Input:
C     LU4...  Number of logical unit corresponding to the file with
C             triangles whose rays belong each to other.
C     LU5...  Number of logical unit corresponding to the file with
C             triangles whose rays do not belong each to other.
C ...........................
C Common block /PTCB/ to store the information about triangles and rays:
      INCLUDE 'pfatubes.inc'
C pfatubes.inc.
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER II,I1,I2,I3,J1,J2,J3,IA
C
C-----------------------------------------------------------------------
C     Setting output format:
   1  CONTINUE
      IF (IRAYMA.GT.IFORM) THEN
        IFORM=IFORM*10+9
        IF (IFORM.GT.99999999) THEN
C         PFATUBES-10
          CALL ERROR('PFATUBES-10: Insufficient format for triangles')
C         Maximum space for indices of triangles is 9 positions
C         including the spaces.
        ENDIF
        FORMAT(3:3)=CHAR(ICHAR(FORMAT(3:3))+1)
        FORMAT(6:6)=FORMAT(3:3)
        FORMAT(9:9)=FORMAT(3:3)
        GOTO 1
      ENDIF
C
C     Writing the triangles of the elementary wave:
      DO 10, II=0,6*(NT-1),6
C       Indices of two sets of rays of the triangle:
        I1=IRAM(II+1)*2-1
        I2=IRAM(II+2)*2-1
        I3=IRAM(II+3)*2-1
        J1=IRAM(II+1)*2
        J2=IRAM(II+2)*2
        J3=IRAM(II+3)*2
        IF (IRAM(II+4)*IRAM(II+5)*IRAM(II+6).NE.1) THEN
C         Triangle where the rays do not belong each to other:
          WRITE(LU5,FORMAT) I1,I2,I3
          WRITE(LU5,FORMAT) J1,J2,J3
        ELSE
C         Triangle where the rays belong each to other:
          IF (IRAM(II+4).EQ.-1) THEN
            IA=I2
            I2=J2
            J2=IA
            IA=I3
            I3=J3
            J3=IA
          ENDIF
          IF (IRAM(II+5).EQ.-1) THEN
            IA=I3
            I3=J3
            J3=IA
          ENDIF
          WRITE(LU4,FORMAT) I1,I2,I3
          WRITE(LU4,FORMAT) J1,J2,J3
        ENDIF
  10  CONTINUE
C
      IF (ISRCN.NE.0) THEN
C       Writing index of the source and of the wave for the next
C       elementary wave:
        IFORM=99999
        FORMAT='(I6,I6,I6)'
        WRITE(LU4,FORMAT) 0,ISRCN,IWAVEN
        ISRCN =0
        IWAVEN=0
      ENDIF
C
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE PTREAD(LU1,LU2,IT1)
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 LU1,LU2,IT1
C Input:
C     LU1...  Number of logical unit corresponding to the file with
C             the quantities stored along rays.
C     LU2...  Number of logical unit corresponding to the file with
C             the quantities at the initial points of rays,
C             corresponding to file LU1.
C     IT1...  Position of the first ray of the actually processed
C             triangle in the array IRAM.
C No output.
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 /PTCB/ to store the information about triangles and rays:
      INCLUDE 'pfatubes.inc'
C pfatubes.inc.
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER IRAY0,IWAVE0,I1
      CHARACTER*52 TXTERR
C     IRAY0.. Index of the last ray read in into the array RAM.
C
C-----------------------------------------------------------------------
C     Loop for the points of rays:
   10 CONTINUE
        IF ((NRAMP+2*NQ).GT.MRAMP) THEN
C         Empting the memory:
          CALL PTERAS
          IF ((NRAMP+2*NQ).GT.MRAMP) THEN
C           PFATUBES-11
            CALL ERROR('PFATUBES-11: Small array RAM')
          ENDIF
        ENDIF
C       Reading the results of the complete ray tracing:
        IRAY0=IRAY
        IWAVE0=IWAVE
        IF ((IRAY0.EQ.0).AND.(IWAVE0.EQ.0)) THEN
C         Reading the first point on a ray of the first wave:
          CALL AP00(0,LU1,LU2)
          IF (IWAVE.LT.1) GOTO 20
        ELSEIF (IWAVE.EQ.IWAVE0) THEN
C         Reading the next point on a ray of the actual wave:
          CALL AP00(0,LU1,LU2)
          IF (IWAVE.NE.IWAVE0) GOTO 20
        ENDIF
        IF (IRAY.LT.IRAYMI) GOTO 10
        IF (IRAY.GT.IRAYMA) GOTO 10
        IF ((IRAY-IRAY0).GE.2) THEN
C         Some rays skipped by AP00:
          DO 15, I1=IRAY0+1,IRAY-1
            IF (I1.GE.IRAYMI) THEN
              IRAM(I1-ORAYA)=IRAM(IRAY0-ORAYA)
            ENDIF
  15      CONTINUE
        ENDIF
        IF (IRAM(IRAY-ORAYE).NE.0) THEN
C         Writing the results of the complete ray tracing - recording
C         new point on a ray:
          IF (IPT.LE.1) THEN
            IF ((ISHEET.EQ.0).AND.(IRAY.EQ.IRAYMI)) THEN
C             PFATUBES-12
              CALL WARN
     *        ('PFATUBES-12: a ray with history 0 observed')
C             All the rays are probably of the same history 0. Only the
C             initial-value ray tracing was performed.  Width and shape
C             of ray tubes were not controlled.  The interpolation
C             in such a case makes sense only in smooth models.
C             Check the value of the parameter IPOINT in CRT input data
C             RPAR (4).
            ENDIF
C           New ray - recording the initial point:
C           Coordinates:
            RAM(NRAMP+1)=YI(3)
            RAM(NRAMP+2)=YI(4)
            RAM(NRAMP+3)=YI(5)
C           Index of surface:
            IRAM(NRAMP+4)=0
C           Sequential index of point:
            IRAM(NRAMP+5)=0
C           Complex-valued polarization vector:
            RAM(NRAMP+ 6)=0.
            RAM(NRAMP+ 7)=0.
            RAM(NRAMP+ 8)=0.
            RAM(NRAMP+ 9)=0.
            RAM(NRAMP+10)=0.
            RAM(NRAMP+11)=0.
            NRAMP=NRAMP+NQ
          ENDIF
C         Recording the new point:
C         Coordinates:
          RAM(NRAMP+1)=Y(3)
          RAM(NRAMP+2)=Y(4)
          RAM(NRAMP+3)=Y(5)
C         Index of surface:
          IRAM(NRAMP+4)=ISRF
C         Sequential index of point:
          IRAM(NRAMP+5)=IPT
C         Complex-valued polarization vector:
          CALL PTPOLA
          NRAMP=NRAMP+NQ
        ENDIF
        IRAM(IRAY-ORAYA)=NRAMP
CCC     IF ((IPT.LE.1).AND.(IRAY.GT.IRAM(IT1)).AND.(IRAY.NE.IRAYMA))
        IF ((NRAMP+2*NQ).GT.MRAMP) THEN
          IF (IRAY.LE.2*IRAM(IT1)) THEN
C           PFATUBES-13
            CALL ERROR('PFATUBES-13: Array RAM too small')
C           Reading of the rays needs to be terminated because
C           the array RAM is full, but the rays of the triangle IT1
C           were not yet read. The dimension MARAM of the array RAM
C           needs to be enlarged.
          ENDIF
          RETURN
        ENDIF
      GOTO 10
C
  20  CONTINUE
C     End of rays:
      IF (IRAY0.LT.IRAYMA) THEN
C       Some rays at the end of the CRT output file missing:
        DO 22, I1=IRAY0+1,IRAYMA
          IF (I1.GE.IRAYMI) THEN
            IRAM(I1-ORAYA)=NRAMP
          ENDIF
  22    CONTINUE
      ENDIF
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE PTPOLA
C
C----------------------------------------------------------------------
C Subroutine to calculate the complex-valued polarization vectors
C in the point of the ray which is currently stored in the common
C block POINTC.
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 /PTCB/ to store the information about triangles and rays:
      INCLUDE 'pfatubes.inc'
C pfatubes.inc.
C.......................................................................
C
C     Auxiliary storage locations:
      REAL PR11,PI11,PR21,PI21,PR12,PI12,PR22,PI22,VAR11,VAR22,VAR12,
     *     VAI12,VAR21,VAI21,AUX,AR1,AI1,AR2,AI2,G11,G12,G13,G21,G22,
     *     G23,ER1,ER2,ER3,EI1,EI2,EI3
C
C-----------------------------------------------------------------------
      PR11=Y(28)
      PI11=Y(29)
      PR21=Y(30)
      PI21=Y(31)
      PR12=Y(32)
      PI12=Y(33)
      PR22=Y(34)
      PI22=Y(35)
      VAR11=PR11**2+PI11**2 + PR12**2+PI12**2
      VAR22=PR21**2+PI21**2 + PR22**2+PI22**2
      VAR12=PR11*PR21+PI11*PI21 + PR12*PR22+PI12*PI22
      VAI12=PI11*PR21-PR11*PI21 + PI12*PR22-PR12*PI22
      VAR21=PR21*PR11+PI21*PI11 + PR22*PR12+PI22*PI12
      VAI21=PI21*PR11-PR21*PI11 + PI22*PR12-PR22*PI12
      IF (VAR11.GE.VAR22) THEN
        AUX=1./SQRT(VAR11**2+VAR21**2+VAI21**2)
        AR1=AUX*VAR11
        AI1=0.
        AR2=AUX*VAR21
        AI2=AUX*VAI21
      ELSE
        AUX=1./SQRT(VAR22**2+VAR12**2+VAI12**2)
        AR1=AUX*VAR12
        AI1=AUX*VAI12
        AR2=AUX*VAR22
        AI2=0.
      ENDIF
      IF (NPV.LT.6) THEN
C       PFATUBES-14
        CALL ERROR('PFATUBES-14: Missing polarization vectors')
C       The S-vave reference polarization vectors need to be defined.
      ENDIF
      G11=PV(1)
      G12=PV(2)
      G13=PV(3)
      G21=PV(4)
      G22=PV(5)
      G23=PV(6)
      ER1=G11*AR1+G21*AR2
      ER2=G12*AR1+G22*AR2
      ER3=G13*AR1+G23*AR2
      EI1=G11*AI1+G21*AI2
      EI2=G12*AI1+G22*AI2
      EI3=G13*AI1+G23*AI2
      RAM(NRAMP+ 6)=ER1
      RAM(NRAMP+ 7)=EI1
      RAM(NRAMP+ 8)=ER2
      RAM(NRAMP+ 9)=EI2
      RAM(NRAMP+10)=ER3
      RAM(NRAMP+11)=EI3
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE PTERAS
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 ...........................
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 /PTCB/ to store the information about triangles and rays:
      INCLUDE 'pfatubes.inc'
C pfatubes.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-----------------------------------------------------------------------
C     PFATUBES-15
      CALL ERROR('PFATUBES-15: PTERAS called')
C     This subroutine is not yet debugged.
      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
            IRAM(J1)=IRAM(I2)
   10     CONTINUE
        ELSE
C         This ray is to be erased:
          IRAM(I1-ORAYE)=0
        ENDIF
        IADDRP=IRAM(I1-ORAYA)
        IRAM(I1-ORAYA)=J1
   20 CONTINUE
      NRAMP=J1
      RETURN
      END
C=======================================================================
C
      SUBROUTINE PTINDX(IT1)
C
C----------------------------------------------------------------------
C Subroutine for calculation of indices of sides of the triangle
C
      INTEGER IT1
C Input:
C     IT1...  The address of the index of the first ray of the triangle,
C             for which the indices of sides are to be computed.
C No output.
C
C ...........................
      EXTERNAL PTLCOS
      LOGICAL PTLCOS
C ...........................
C Common block /PTCB/ to store the information about triangles and rays:
      INCLUDE 'pfatubes.inc'
C pfatubes.inc.
C.......................................................................
C     Auxiliary storage locations:
      INTEGER I1B,I2B,I3B,I1C,I2C,I3C,J1B,J2B,J3B,J1C,J2C,J3C
      INTEGER I1MA,I2MA,I3MA,J1,J2,ISID1,ISID2,ISID3
      LOGICAL LFIRST,LI1I2,LJ1J2,LI1J2,LJ1I2
C     I1B,I2B,I3B,I1C,I2C,I3C... Integers controling main loop over
C             points on the rays (along ray tube). Numbers 1,2,3 denote
C             the first, second and third ray, character B denotes
C             bottom of the ray cell and C denotes top of the ray cell.
C             Each integer contains the address just before
C             the parameters of the corresponding point:
C             the first parameter of the first point: RAM(I1B+1)
C     J1...   Counts the number of identical points of the ray cell
C             (J1=0,1,2,3).
C     J2...   Displays maximum sequential index of a point allowed
C             for actual ray cell.
C
C-----------------------------------------------------------------------
      IF (         (IRAM(2*IRAM(IT1  )-1-ORAYA).EQ.0) .OR.
     *             (IRAM(2*IRAM(IT1  )  -ORAYA).EQ.0) .OR.
     *             (IRAM(2*IRAM(IT1+1)-1-ORAYA).EQ.0) .OR.
     *             (IRAM(2*IRAM(IT1+1)  -ORAYA).EQ.0) .OR.
     *    (L3D.AND.(IRAM(2*IRAM(IT1+2)-1-ORAYA).EQ.0)).OR.
     *    (L3D.AND.(IRAM(2*IRAM(IT1+2)  -ORAYA).EQ.0))) THEN
C       PFATUBES-16
        CALL ERROR
     *     ('PFATUBES-16: Parameters of a ray not found in memory')
C       This error may be caused by wrong file CRT-R with
C       missing ray(s) (e.g. only two-point rays stored).
      ENDIF
C
      I1MA=IRAM(2*IRAM(IT1  )-ORAYA)-NQ
      I2MA=IRAM(2*IRAM(IT1+1)-ORAYA)-NQ
      IF (L3D) I3MA=IRAM(2*IRAM(IT1+2)-ORAYA)-NQ
C     The first ray cell:
      I1B=IRAM(2*IRAM(IT1  )-ORAYA-1)
      I2B=IRAM(2*IRAM(IT1+1)-ORAYA-1)
      IF (L3D) I3B=IRAM(2*IRAM(IT1+2)-ORAYA-1)
      J1B=IRAM(2*IRAM(IT1  )-1-ORAYA-1)
      J2B=IRAM(2*IRAM(IT1+1)-1-ORAYA-1)
      IF (L3D) J3B=IRAM(2*IRAM(IT1+2)-1-ORAYA-1)
      IF ((I1B.GT.I1MA).OR.(I2B.GT.I2MA).OR.(L3D.AND.(I3B.GT.I3MA)))
     *  RETURN
      I1C=I1B
      I2C=I2B
      IF (L3D) I3C=I3B
      J1C=J1B
      J2C=J2B
      IF (L3D) J3C=J3B
      LFIRST=.TRUE.
      ISID1=1
      ISID2=1
      ISID3=1
C     Loop for points on the rays (loop for ray cells):
C     J1 counts the number of points, which were not shifted.
C     J2 displays the sequential number of points, which is the
C        maximum allowed number for current cell.
  10  CONTINUE
        IF (L3D) THEN
          J1=3
          J2=MAX0(IRAM(I1B+5),IRAM(I2B+5),IRAM(I3B+5))
          IF ((IRAM(I1B+5).EQ.IRAM(I2B+5)).AND.
     *        (IRAM(I1B+5).EQ.IRAM(I3B+5)))  J2=J2+1
        ELSE
          J1=2
          J2=MAX0(IRAM(I1B+5),IRAM(I2B+5))
          IF (IRAM(I1B+5).EQ.IRAM(I2B+5))  J2=J2+1
        ENDIF
C       Forming standard ray cells:
        IF ((IRAM(I1B+4).EQ.0).AND.(IRAM(I1B+5).LT.J2)) THEN
          I1C=I1B+NQ
          J1C=J1B+NQ
          J1=J1-1
        ENDIF
        IF ((IRAM(I2B+4).EQ.0).AND.(IRAM(I2B+5).LT.J2)) THEN
          I2C=I2B+NQ
          J2C=J2B+NQ
          J1=J1-1
        ENDIF
        IF (L3D) THEN
          IF ((IRAM(I3B+4).EQ.0).AND.(IRAM(I3B+5).LT.J2)) THEN
            I3C=I3B+NQ
            J3C=J3B+NQ
            J1=J1-1
          ENDIF
        ENDIF
        IF ((L3D.AND.(J1.EQ.3)).OR.(L2D.AND.(J1.EQ.2))) THEN
          IF (L3D.AND.(IRAM(I1B+4).EQ.IRAM(I2B+4)).AND.
     *                (IRAM(I1B+4).EQ.IRAM(I3B+4))) THEN
C           Crossing the interface in 3D:
            I1B=I1B+NQ
            I2B=I2B+NQ
            I3B=I3B+NQ
            J1B=J1B+NQ
            J2B=J2B+NQ
            J3B=J3B+NQ
            J1=0
            IF ((I1B.GT.I1MA).OR.(I2B.GT.I2MA).OR.(I3B.GT.I3MA))
     *        RETURN
            IRAM(I1B+4)=0
            IRAM(I2B+4)=0
            IRAM(I3B+4)=0
            IRAM(J1B+4)=0
            IRAM(J2B+4)=0
            IRAM(J3B+4)=0
            I1C=I1B
            I2C=I2B
            I3C=I3B
            J1C=J1B
            J2C=J2B
            J3C=J3B
CCC         GOTO 10
          ELSEIF (L2D.AND.(IRAM(I1B+4).EQ.IRAM(I2B+4))) THEN
C           Crossing the interface in 2D:
            I1B=I1B+NQ
            I2B=I2B+NQ
            J1B=J1B+NQ
            J2B=J2B+NQ
            J1=0
            IF ((I1B.GT.I1MA).OR.(I2B.GT.I2MA))
     *        RETURN
            IRAM(I1B+4)=0
            IRAM(I2B+4)=0
            IRAM(J1B+4)=0
            IRAM(J2B+4)=0
            I1C=I1B
            I2C=I2B
            J1C=J1B
            J2C=J2B
CCC         GOTO 10
          ELSE
C           Moving the points in a complex block:
C           Forming standard ray cells:
            IF (IRAM(I1B+4).EQ.0) THEN
              I1C=I1B+NQ
              J1C=J1B+NQ
              J1=J1-1
            ENDIF
            IF (IRAM(I2B+4).EQ.0) THEN
              I2C=I2B+NQ
              J2C=J2B+NQ
              J1=J1-1
            ENDIF
            IF (L3D) THEN
              IF (IRAM(I3B+4).EQ.0) THEN
                I3C=I3B+NQ
                J3C=J3B+NQ
                J1=J1-1
              ENDIF
            ENDIF
            IF ((L3D.AND.(J1.EQ.3)).OR.(L2D.AND.(J1.EQ.2))) THEN
              IF ((I1B.GE.I1MA).AND.(I2B.GE.I2MA).AND.(I3B.GE.I3MA).AND.
     *            L3D) RETURN
              IF ((I1B.GE.I1MA).AND.(I2B.GE.I2MA).AND.L2D) RETURN
C             PFATUBES-17
              CALL ERROR
     *        ('PFATUBES-17: The points reached different interfaces')
C             This error should not appear.
C             Please contact the author.
            ENDIF
          ENDIF
        ENDIF
        IF ((I1C.GT.I1MA).OR.(I2C.GT.I2MA).OR.
     *      (L3D.AND.(I3C.GT.I3MA))) THEN
C         PFATUBES-18
          CALL ERROR('PFATUBES-18: point exceeded the ray.')
C         This error should not appear.
C         Please contact the author.
        ENDIF
C
C       The ray cell whose top is formed by points, whose parameters
C       can be found just behind the addresses I1C,I2C,I3C,J1C,J2C,J3C
C       is prepared for comparison of polarization vectors:
        IF (L2D) THEN
C         2-D case:
C         PFATUBES-19
          CALL ERROR('PFATUBES-19: 2-D calculation')
C         Not yet coded.
        ENDIF
        IF (L3D) THEN
C         Comparison of polarization vectors:
          IF (ISID1.NE.0) THEN
C           First side:
            LI1I2=PTLCOS(I1C,I2C)
            LJ1J2=PTLCOS(J1C,J2C)
            LI1J2=PTLCOS(I1C,J2C)
            LJ1I2=PTLCOS(J1C,I2C)
            IF     (LI1I2.AND.LJ1J2.AND.
     *              (.NOT.LI1J2).AND.(.NOT.LJ1I2)) THEN
              ISID1=1
            ELSEIF (LI1J2.AND.LJ1I2.AND.
     *              (.NOT.LI1I2).AND.(.NOT.LJ1J2)) THEN
              ISID1=-1
            ELSE
              ISID1=0
            ENDIF
          ENDIF
          IF (ISID2.NE.0) THEN
C           Second side:
            LI1I2=PTLCOS(I2C,I3C)
            LJ1J2=PTLCOS(J2C,J3C)
            LI1J2=PTLCOS(I2C,J3C)
            LJ1I2=PTLCOS(J2C,I3C)
            IF     (LI1I2.AND.LJ1J2.AND.
     *              (.NOT.LI1J2).AND.(.NOT.LJ1I2)) THEN
              ISID2=1
            ELSEIF (LI1J2.AND.LJ1I2.AND.
     *              (.NOT.LI1I2).AND.(.NOT.LJ1J2)) THEN
              ISID2=-1
            ELSE
              ISID2=0
            ENDIF
          ENDIF
          IF (ISID3.NE.0) THEN
C           Third side:
            LI1I2=PTLCOS(I3C,I1C)
            LJ1J2=PTLCOS(J3C,J1C)
            LI1J2=PTLCOS(I3C,J1C)
            LJ1I2=PTLCOS(J3C,I1C)
            IF     (LI1I2.AND.LJ1J2.AND.
     *              (.NOT.LI1J2).AND.(.NOT.LJ1I2)) THEN
              ISID3=1
            ELSEIF (LI1J2.AND.LJ1I2.AND.
     *              (.NOT.LI1I2).AND.(.NOT.LJ1J2)) THEN
              ISID3=-1
            ELSE
              ISID3=0
            ENDIF
          ENDIF
          IF (LFIRST) THEN
C           First points of the ray tube:
            IRAM(IT1+3)=ISID1
            IRAM(IT1+4)=ISID2
            IRAM(IT1+5)=ISID3
            LFIRST=.FALSE.
          ELSE
C           Other than the first points of the ray tube:
            IF (IRAM(IT1+3).NE.ISID1) THEN
              ISID1=0
              IRAM(IT1+3)=ISID1
            ENDIF
            IF (IRAM(IT1+4).NE.ISID2) THEN
              ISID2=0
              IRAM(IT1+4)=ISID2
            ENDIF
            IF (IRAM(IT1+5).NE.ISID3) THEN
              ISID3=0
              IRAM(IT1+5)=ISID3
            ENDIF
          ENDIF
          IF ((ISID1.EQ.0).AND.(ISID2.EQ.0).AND.(ISID3.EQ.0)) THEN
C           All sides formed by rays which do not belong each to other:
            RETURN
          ENDIF
        ENDIF
C       Continuing with the next ray cell:
        I1B=I1C
        I2B=I2C
        IF (L3D) I3B=I3C
        J1B=J1C
        J2B=J2C
        IF (L3D) J3B=J3C
      GOTO 10
C     End of the loop along the ray tube.
      END
C
C=======================================================================
C
      LOGICAL FUNCTION PTLCOS(IP1,IP2)
C
C----------------------------------------------------------------------
      INTEGER IP1,IP2
C     Subroutine to calculate cosinus of the angle of the
C     polarization vectors of two points on the rays, and to decide
C     whether the two rays belong together or not.
C
C Input:
C     IP1,IP2 ... Addresses of the two points on the rays in array RAM.
C Output:
C     PTLCOS ... .TRUE. if the points belong each to other
C                .FALSE. otherwise
C ...........................
C Common block /PTCB/ to store the information about triangles and rays:
      INCLUDE 'pfatubes.inc'
C pfatubes.inc.
C.......................................................................
C     Auxiliary storage locations:
      REAL E1R1,E1I1,E1R2,E1I2,E1R3,E1I3,E2R1,E2I1,E2R2,E2I2,E2R3,E2I3,
     *     AUX1,AUX2,COSFI2
C     E1R1... First polarization vector, real part, first component
C     E2I3... Second polariz. vector, imaginary part, third component
C-----------------------------------------------------------------------
      E1R1=RAM(IP1+ 6)
      E1I1=RAM(IP1+ 7)
      E1R2=RAM(IP1+ 8)
      E1I2=RAM(IP1+ 9)
      E1R3=RAM(IP1+10)
      E1I3=RAM(IP1+11)
      E2R1=RAM(IP2+ 6)
      E2I1=RAM(IP2+ 7)
      E2R2=RAM(IP2+ 8)
      E2I2=RAM(IP2+ 9)
      E2R3=RAM(IP2+10)
      E2I3=RAM(IP2+11)
      AUX1=E1R1*E2R1+E1R2*E2R2+E1R3*E2R3+
     *     E1I1*E2I1+E1I2*E2I2+E1I3*E2I3
      AUX2=E1R1*E2I1+E1R2*E2I2+E1R3*E2I3-
     *    (E1I1*E2R1+E1I2*E2R2+E1I3*E2R3)
      COSFI2=AUX1**2+AUX2**2
      IF (COSFI2.GT.0.5) THEN
        PTLCOS=.TRUE.
      ELSE
        PTLCOS=.FALSE.
      ENDIF
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'indexi.for'
C     indexi.for.
C
C=======================================================================
C