C
C Subroutine file 'code.for' - codes for elementary waves.
C
C Version: 6.60
C Date: 2011, June 8
C
C Coded by Ludek Klimes
C     Department of Geophysics, Charles University in Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/klimes.htm
C
C.......................................................................
C
C This file consists of:
C     Codes for elementary waves - general description
C             Codes for elementary waves
C     CODE1...Subroutine designed to read the input data for the codes
C             of elementary waves and to store them in the common block
C             /COD/.
C             CODE1
C     LCODE...Integer function returning the length of the code of the
C             current elementary wave.
C             LCODE
C     NELEM...Integer function returning the number of ray elements
C             corresponding to the given position in the code of the
C             elementary wave.
C             NELEM
C     CODE... Subroutine designed to transform the used numerical code
C             of the elementary wave under consideration into
C             instructions specifying the behaviour of the wave at the
C             initial points of rays and at all points of incidence of
C             the rays at interfaces (boundaries of complex blocks).
C             CODE
C
C.......................................................................
C
C                                                    
C Input data CODE for the codes of elementary waves:
C     The data are read in by the list directed input (free format).  In
C     the list of input data below, each numbered paragraph indicates
C     the beginning of a new input operation (new READ statement).  If
C     the first letter of the symbolic name of the input variable is
C     I-N, the corresponding value in input data must be of the type
C     INTEGER.  Otherwise (except TEXTC), the input parameter is of the
C     type REAL.
C (1) TEXTC
C     String describing the data.  Only the first 80 characters of the
C     string are significant.
C (2) NSKIP
C     The number of elementary waves to be skipped over.  The
C     (NSKIP+1)-th elementary wave in these input data will be the first
C     computed wave.  Default value: 0.
C For each elementary wave IWAVE the following data (3):
C (3) KODTYP,KODE
C     KODTYP..Determines the type of the code
C             1... Block-surface code             1
C             2... Basic block code               2
C             3... Compound-element block code    3
C             4... Generalized layer code         4
C             5... Reflection-transmission code   5
C     KODE... The code of the elementary wave in the form of a finite
C             sequence of integers terminated by a slash.
C     Attention:
C             Generally, different types of the code may be combined
C             within the same input data, but the user has then be very
C             careful.  If combining the different types of the code,
C             output parameters IWAVE0 and IKODE of subroutine CODE1 may
C             have no sense.  Then:
C             If the parameter MODCRT of the data set 'DCRT' (see source
C               code file 'ray.for') is zero, the output parameters
C               IWAVE0 and IKODE of subroutine CODE1 are hopefully
C               unused.
C             If the parameter MODCRT of the data set 'DCRT' is not
C               zero, mixing the types of the wave code should be
C               avoided unless the user finds a way how to do it
C               reasonably.
C (4) a slash.
C Example of data set CODE
C
C.......................................................................
C
C Storage in the memory:
C     The input data (2) for the computed elementary wave are stored in
C     the common block /COD/ defined in the include file 'code.inc'.
C     code.inc
C
C=======================================================================
C
C                                                   
C Codes for elementary waves - general description
C
C     We consider ordinary seismic body waves, such as refracted,
C     primarily or multiply reflected, possibly converted waves.  In
C     general, incidence of a wave at an interface (boundary of a
C     complex block) produces four waves, reflected P and S, and
C     transmitted P and S waves.  When performing complete ray tracing,
C     we must know a priori which of the four generated waves to follow.
C     This decision must be made at each interface.  The alphanumeric
C     string specifying the behaviour of a ray from its initial point to
C     its endpoint is a 'code'.  In the following, we shall consider the
C     code to be a sequence of nonzero integers.
C
C     The term 'elementary wave' does not have unique meaning in the
C     literature.  Here we apply the term elementary wave to that part
C     of the wavefield that is described by one specific code.  Since
C     there may be various types of codes, there is also a variety of
C     divisions of the wave field into elementary waves.
C
C     We introduce the term 'element of a ray', which has an important
C     meaning in the construction of codes.  By an element of a ray, we
C     denote that part of the ray that is situated in one complex block
C     between two successive points of reflection/transmission, or
C     between the initial point or endpoint of the ray and the closest
C     point of reflection/transmission, or between the initial point and
C     the endpoint of the ray, if the ray is entirely situated in one
C     complex block.
C
C     In the following, we present several possible types of numerical
C     codes of elementary waves.
C
C Examples of code types
C                                                   
C     (1) Block-surface code
C             Two integers are used for any element of the ray.  The
C             absolute value of the first integer specifies the index of
C             the complex block, in which the element of the ray is
C             situated, the sign specifies the type of the wave (plus
C             sign...  P wave, minus sign...  S wave).  The second
C             number is the index of the smooth surface, at which the
C             element terminates.  The code of the ray is a chain of the
C             above doublets describing the ray from its initial point
C             to its endpoint.  The index of a surface, at which the
C             last element of a ray terminates, need not be specified.
C             In this case the ray is allowed to terminate on any
C             surface bounding the complex block in which the last
C             element of the ray is situated.  Application of this code
C             is straightforward.
C
C             If any complex block is bounded by several smooth
C             surfaces, the block-surface code may divide artificially
C             the wave field into several elementary waves.  This would
C             result in repeated application of the complete ray tracing
C             to all of the elementary waves, what may be time
C             consuming.  In such situations, it is therefore desirable
C             to use other codes, for which the elementary wave has a
C             more general meaning.
C
C                                                   
C     (2) Basic block code
C             This code may be obtained from the block-surface code if
C             the indices of the surfaces, at which individual elements
C             of a ray terminate, are omitted.  Then, any element of the
C             ray is described by a single integer specifying the
C             complex block and wave type, and it may terminate at any
C             surface bounding the complex block.  The code of the ray
C             is again a chain of the numbers describing successively
C             its individual elements.
C
C             Elementary waves specified by the block-surface code,
C             passing through the same complex blocks, and reflected
C             from different surfaces bounding a complex block, are
C             united into one elementary wave in the basic block code.
C             Similarly, elementary waves specified by the block-surface
C             code and transmitted from a block into the neighbouring
C             block across different surfaces separating these blocks,
C             are united into one elementary wave by the basic block
C             code.
C
C                                                   
C     (3) Compound-element block code
C             Let us introduce the terms simple and compound element of
C             a ray.  For this purpose, we call the 'lower-index
C             boundary' of a complex block that part of its boundary,
C             which separates the complex block either from complex
C             blocks with lower indices or from a free space.  The
C             remaining part of the boundary is called the
C             'higher-index boundary' of the complex block.  The initial
C             point of a ray is treated in the same way as the points
C             situated at a lower-index boundary.
C
C             An element of a ray is called 'simple element' if its
C             initial point (e.g. the point where the wave enters the
C             corresponding complex block) and its endpoint (e.g. the
C             point at which the wave leaves the complex block) are
C             situated one on the lower-index and the other on the
C             higher-index boundary of the complex block.  The
C             'compound element' is such an element for which its
C             initial point and its endpoint are both situated either
C             on the lower-index or on the higher-index boundary of the
C             complex block.  The compound element is formally
C             considered as two simple elements.
C
C             Often it is convenient to work with waves, the rays of
C             which have either a compound element in a complex block or
C             two simple elements of the same wave type (unconverted
C             reflection) in the same block, as with one elementary
C             wave.  The compound-element block code makes a division of
C             the wave field into such elementary waves possible.
C
C             We introduce the compound-element block code as follows.
C             For any simple element of a ray one number is used.  The
C             absolute value of this number specifies the index of the
C             complex block, in which the simple element is situated.
C             The sign of this number specifies the type of the wave
C             (plus sign...  P wave, minus sign...  S wave).  For any
C             compound element of a ray, two identical numbers are used.
C             A compound element and a doublet of simple elements,
C             described by the same numbers, are not distinguished.  The
C             code of a ray is a chain of these numbers describing the
C             ray from its initial point to its endpoint.  The use of
C             the compound-element block code leads to more general
C             elementary waves.
C
C                                                   
C     (4) Generalized layer code
C             Another code may be obtained from the compound-element
C             block code if the modified interpretation of the model
C             structure is used.  The modified interpretation consists
C             in assuming an existence of fictitious parts of blocks of
C             a zero thickness situated between every neighbouring
C             complex blocks, the indices of which are not sequentially
C             ordered.  The number of fictitious blocks is considered to
C             be just the number, which is necessary to fill the gap
C             between indices of the two neighbouring complex blocks.
C             The ray can only pass through the fictitious blocks, no
C             reflection in fictitious blocks being allowed.  Similarly,
C             no conversion is allowed neither at interfaces between
C             fictitious blocks nor at the interface at which the ray
C             leaves a fictitious block (conversion may take place only
C             at the interface at which the ray enters the fictitious
C             block).  If any of these prohibited situations is
C             specified by the code, it leads to the termination of
C             computations of the corresponding ray.
C
C             This code is a generalization of the code used in the 2-D
C             program package 'SEIS83' and is described in Cerveny,
C             Molotkov and Psencik: Ray Method in Seismology, Charles
C             University, Prague 1977.  It enables any block structure
C             to be interpreted as locally layered structure with
C             fictitious layers.  Thus, the routines for automatic or
C             semiautomatic generation of ray codes for layered
C             structures may be used for this code even in general block
C             structures.  Its effective use is conditioned by proper
C             indexing of complex blocks, which minimizes the number of
C             fictitious layers.  For example, in a layered structure,
C             the blocks (layers) should be indexed sequentially from
C             the top to the bottom.
C
C                                                   
C     (5) Reflection-transmission code
C             For any element of a ray one integer is used.  The first
C             element of the ray, i.e.  the element containing the
C             initial point of the ray, is denoted by 1 for P and by -1
C             for S wave.  Any other element is denoted by 1 or -1 (P or
C             S wave) if the ray is transmitted at the initial point of
C             the element, and by 2 or -2 (P or S wave) if the ray is
C             reflected at the initial point of the element.  The code
C             of a ray is a chain of the above numbers corresponding to
C             individual elements of the ray from its initial point to
C             its endpoint.  In a layered model, this code is very
C             convenient for the description of refracted and primary
C             reflected waves.  The application of this code is
C             straightforward.
C
C Specification of an elementary wave
C     an elementary wave is specified by the following data:
C     (1) Integer KODTYP, which determines the type of the code
C             1... Block-surface code,
C             2... Basic block code,
C             3... Compound-element block code,
C             4... Generalized layer code,
C             5... Reflection-transmission code.
C     (2) Integer array KODE containing the code of the elementary
C             wave.  The code is a finite sequence of nonzero integers.
C             A zero indicates the end of a code in the input data.
C     The data (1) and (2) should be stored in common block
C     -------------------------------
C     COMMON /COD/ KODTYP,KODE(MKODE), see 'code.inc'.
C     code.inc
C     -------------------------------
C     The dimension MKODE of the array KODE may be adjusted by the user.
C
C Date: 1988, June 3
C Written by Vlastislav Cerveny, Ludek Klimes, Ivan Psencik
C
C=======================================================================
C
C     
C
      SUBROUTINE CODE1(LUN,IWAVE,IWAVE0,IKODE)
      INTEGER LUN,IWAVE,IWAVE0,IKODE
C
C This subroutine is called when starting the computation of a new
C elementary wave.  It stores the code of the elementary wave into the
C common block /COD/.
C
C Input:
C     LUN...  Logical unit number of the external input device
C             containing the input data.
C     IWAVE...Zero when starting the complete ray tracing program,
C             otherwise the index of the last elementary wave which has
C             been computed (i.e. the output from the previous
C             invocation of the subroutine CODE1).
C
C Output:
C     IWAVE...Zero if all required elementary waves are computed and the
C             complete ray tracing program will be terminated.
C             Otherwise, the index of the elementary wave which will be
C             computed (i.e. NSKIP+1 for the first computed elementary
C             wave, where NSKIP is the number of elementary waves having
C             been skipped over, otherwise the input value increased by
C             one).
C     IWAVE0..Index of the already computed elementary wave having the
C             most numerous common elements with the current elementary
C             wave.  In the case of several possibilities, the first
C             computed wave of them is taken.
C     IKODE...The length of the common part of the codes of the IWAVE-th
C             and IWAVE0-th elementary waves.
C
C Common block /COD/:
      INCLUDE 'code.inc'
C     code.inc
C All the storage locations of the common block are defined in this
C subroutine.
C
C Subroutines and external functions required:
      EXTERNAL LCODE
      INTEGER  LCODE
C     LCODE.. This file.
C
C Date: 2011, February 4
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Storage for previous codes:
      INTEGER MKODES
      PARAMETER (MKODES=1024)
      INTEGER KODES(0:MKODES)
      SAVE KODES
C
C     Auxiliary storage locations:
      CHARACTER*80 TEXTC
      INTEGER NSKIP,I,J,K,L
      SAVE NSKIP
C
C     TEXTC...The name of the data. String of 80 characters.
C     NSKIP...The number of elementary waves to be skipped over.
C     I,J...  Auxiliary loop variables.
C     K...    Auxiliary variable - index in array KODES.
C     L...    Auxiliary variable - length of the current code.
C
C.......................................................................
C
      IF(IWAVE.EQ.0) THEN
C       Reading the name of the input data
        READ(LUN,*) TEXTC
C
C       Reading the number of elementary waves to be skipped over
        NSKIP=0
        READ(LUN,*) NSKIP
C       No codes of elementary waves are read and stored now
        KODES(0)=0
      END IF
C
C     Reading the code of the elementary wave
   10 CONTINUE
        KODTYP=0
        DO 11 I=1,MKODE
          KODE(I)=0
   11   CONTINUE
        READ(LUN,*) KODTYP,KODE
        IF(KODTYP.EQ.0) THEN
          IWAVE=0
          RETURN
        ELSE
          K=KODES(IWAVE)+1
          L=LCODE()
          IWAVE=IWAVE+1
          IF(K+L.GT.MKODES) THEN
C           401
            CALL ERROR('401 in CODE1: Insufficient memory for codes')
C           The dimension MKODES of the array KODES(0:MKODES) in
C           this subroutine should be increased.  The dimension MKODES
C           should at least equal to the number of elementary waves
C           plus the total number of all indices forming the codes of
C           the waves.
C           Note:
C           If the parameter MODCRT of the data set 'DCRT' (see source
C           code file 'ray.for') is zero, the output parameters IWAVE0
C           and IKODE of this subroutine are likely unused.  Then,
C           instead of increasing MKODES, the above statements
C           generating this error message may be disabled to save
C           the computer memory.  Then this subroutine will work,
C           but considering only the first several elementary waves
C           when calculating output values of IWAVE0 and IKODE.
          ELSE
C           Storing the codes
            DO 12 I=L,1,-1
              KODES(K+I)=KODE(I)
   12       CONTINUE
            DO 13 I=K-1,IWAVE,-1
              KODES(I+1)=KODES(I)
   13       CONTINUE
            KODES(IWAVE)=K+L
            DO 14 I=IWAVE-1,0,-1
              KODES(I)=KODES(I)+1
   14       CONTINUE
          END IF
        END IF
      IF(IWAVE.LE.NSKIP) GO TO 10
C
C     Determining elementary wave having the most numerous common
C     elements with the current elementary wave
      IKODE=0
      IWAVE0=0
      L=LCODE()
      DO 39 J=1,MIN0(KODES(0),IWAVE-1)
        K=KODES(J-1)
        DO 31 I=1,MIN0(KODES(J)-K,L)
          IF(KODE(I).NE.KODES(K+I)) THEN
            GO TO 32
          END IF
   31   CONTINUE
   32   CONTINUE
        I=I-1
        IF(I.GT.IKODE) THEN
          IKODE=I
          IWAVE0=J
        END IF
   39 CONTINUE
      RETURN
      END
C
C=======================================================================
C
C     
C
      INTEGER FUNCTION LCODE()
C
C Integer function LCODE is designed to return the length of the code of
C the current elementary wave.
C
C No input.
C
C Output:
C     LCODE...Length of the code of the current elementary wave, i.e.
C             the count of the numeric items in the array KODE which
C             describes the behaviour of a ray at interfaces.  The array
C             is called here the code of elementary wave.
C
C Common block /COD/:
      INCLUDE 'code.inc'
C     code.inc
C None of the storage locations of the common block are altered.
C
C No subroutines and external functions required.
C
C Date: 1996, June 12
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
C
      DO 1 I=1,MKODE
        IF(KODE(I).EQ.0) THEN
          LCODE=I-1
          GO TO 2
        END IF
    1 CONTINUE
      LCODE=MKODE
    2 CONTINUE
      RETURN
      END
C
C=======================================================================
C
C     
C
      INTEGER FUNCTION NELEM(KODIND)
      INTEGER KODIND
C
C This function is designed to return the maximum possible number of ray
C elements corresponding to the given position in the code of the
C elementary wave.
C
C Input:
C     KODIND..Position in the code (index in array KODE)
C None of the input parameters are altered.
C
C Output:
C     NELEM...Number of ray elements between the initial point of the
C             ray and the end of the element corresponding to the given
C             position in the code.  Note that, for a particular ray,
C             some of these possible ray elements may be left out or
C             coupled into one virtual element.
C
C Common block /COD/:
      INCLUDE 'code.inc'
C     code.inc
C None of the storage locations of the common block are altered.
C
C No subroutines and external functions required.
C
C Date: 1989, December 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     No auxiliary storage locations.
C
      IF(KODTYP.EQ.1) THEN
        NELEM=(KODIND+1)/2
      ELSE
        NELEM=KODIND
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE CODE(IY,KODIND,ICBNEW,IEND)
      INTEGER IY(12),KODIND,ICBNEW,IEND
C
C This subroutine transforms the used numerical code of the elementary
C wave under consideration into instructions specifying the behaviour
C of the wave at the initial points of rays and at all points of
C incidence of the rays at interfaces (boundaries of complex blocks).
C Thus, it is first called at the initial point of a ray and then
C successively at all points of incidence.
C
C Input:
C     IY...   Integer array. Its numerical storage units 2, 3, 5, 6, 8
C             must be defined as follows:
C             IY(2)... Position in the code (index in array KODE)
C               corresponding to the last computed element of a ray.
C               IY(2)=0 at the initial point of the ray,
C               IY(2)=KODIND from the last invocation of subroutine CODE
C                 at other points of the ray.
C             IY(3)=ICB0... Index of the neighbouring complex block,
C               from which the ray entered the complex block in which
C               the last computed element of the ray is situated.
C               IY(3)=0 before leaving the complex block, in which the
C               initial point of the ray is situated.
C               At the initial point of the ray (IY(2)=0), IY(3) is
C               ignored.
C             IY(5)=ICB1... Index of the complex block containing the
C               computed element of the ray, supplemented by the sign
C               '+' for P wave and sign '-' for S wave.
C               ICB1 is ignored at the initial point of the ray
C               (IY(2)=0).
C             IY(6)=ISRF... Index of the surface at which the endpoint
C               of the last computed element of the ray is situated.
C               The sign of IY(6) is ignored.
C               IY(6)=0 at the initial point of the ray.
C             IY(8)=ICB2... Index of the complex block touching the
C               complex block ICB1 from the other side of the surface
C               ISRF at the endpoint of the last computed element of the
C               ray. IY(8)=0 for a free space.
C               At the initial point of the ray, ICB2 is the index of
C               the complex block containing the initial point.
C The input parameter is not altered.
C
C Output:
C     KODIND..Position in the code (index in the array KODE)
C             corresponding to the next element of the ray.
C     ICBNEW..The index of the complex block in which the next
C             element of the ray is to be situated, supplemented
C             by the sign "+" for P wave or "-" for S wave.
C     IEND... Information on the process of the interpretation
C             of the code:
C             IEND.EQ.0... Computation of the ray continues.
C             IEND.NE.0... The computation of the ray terminates.
C               Different values of IEND specify the reason for the
C               termination:
C               10... Ray satisfies the whole code (regular termination
C                 of the ray).
C               21... The point of incidence is situated at a different
C                 surface than that required by the code.
C               22... The next element of the ray is required by the
C                 code to be situated in a complex block that does not
C                 touch the point of incidence.
C               23... Transmission is required by the code at a free
C                 surface, or the initial point of the ray is situated
C                 in free space.
C               30... Reflection or wave conversion at the fictitious
C                 part of the interface.  This reason may occur for the
C                 generalized layer code only (KODTYP=4).
C
C Common block /COD/:
C     ------------------------------------------------------------------
      INCLUDE 'code.inc'
C     code.inc
C     ------------------------------------------------------------------
C     KODTYP..Determines the type of the code:
C             1... Block-surface code,
C             2... Basic block code,
C             3... Compound-element block code,
C             4... Generalized layer code,
C             5... Reflection-transmission code.
C     KODE... Array containing the code of the elementary wave.
C             The code is a finite sequence of nonzero integers.
C             Zero indicates the end of a code.
C None of the storage locations of the common block are altered.
C
C No subroutines and external functions required.
C
C Date: 1990, November 18
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I,J
C
      KODIND=IY(2)
      GO TO (10,20,30,30,50),KODTYP
C
C     Block-surface code:
   10 CONTINUE
      IF(KODIND.GT.0) THEN
C       Check on the surface:
        KODIND=KODIND+1
        IF(KODE(KODIND).EQ.0) THEN
C         End of code
          IEND=10
          RETURN
        ELSE IF(KODE(KODIND).NE.IABS(IY(6))) THEN
C         Wrong surface
          IEND=21
          RETURN
        END IF
      END IF
C     The rest of interpretation is the same as in the basic block code
C
C     Basic block code:
   20 CONTINUE
      KODIND=KODIND+1
      ICBNEW=KODE(KODIND)
      IF(ICBNEW.EQ.0) THEN
C       End of code
        IEND=10
      ELSE IF(IABS(ICBNEW).EQ.IY(8)) THEN
C       Transmission
        IEND=0
      ELSE IF(KODIND.GT.1.AND.IABS(ICBNEW).EQ.IABS(IY(5))) THEN
C       Reflection
        IEND=0
      ELSE
C       Required complex block is not attainable
        IEND=22
      END IF
      RETURN
C
C     Compound-element block code:
   30 CONTINUE
      IF(KODIND.GT.0) THEN
        I=IABS(IY(5))
        J=(I-IABS(IY(3)))*(I-IY(8))
        DO 31 I=KODIND-1,1,-1
          IF(KODE(I).NE.KODE(KODIND)) GO TO 32
          J=-J
   31   CONTINUE
   32   CONTINUE
        IF(J.GT.0) THEN
C         Compound element:
          KODIND=KODIND+1
          IF(KODE(KODIND).NE.IY(5)) THEN
C           Wrong surface
            IEND=21
            RETURN
          END IF
        END IF
      END IF
C     Rest of interpretation is the same as in Basic block code
      IF(KODTYP.EQ.3) GO TO 20
C
C     Generalized layer code:
C     Interpretation of compound elements is the same as in
C     Compound-element block code, then:
      KODIND=KODIND+1
      ICBNEW=KODE(KODIND)
      IF(ICBNEW.EQ.0) THEN
C       End of code
        IEND=10
      ELSE IF(KODIND.EQ.1) THEN
C       Initial point of a ray:
        IF(IABS(ICBNEW).EQ.IY(8)) THEN
C         Initial point of a ray is situated in the required c.b.
          IEND=0
        ELSE
C         Initial point of a ray is not situated in the required c.b.
          IEND=22
        END IF
      ELSE IF(IABS(ICBNEW).EQ.IABS(IY(5))) THEN
C       Reflection
        IEND=0
      ELSE
C       Possible transmission:
        J=ISIGN(1,IY(8)-IABS(IY(5)))
        IF(IABS(ICBNEW).EQ.IABS(IY(5))+J) THEN
C         Loop for fictitious parts of blocks:
          DO 41 I=IABS(IY(5))+J+J,IY(8),J
            KODIND=KODIND+1
            ICBNEW=KODE(KODIND)
            IF(ICBNEW.EQ.0) THEN
C             End of code in the fictitious part of block
              IEND=10
              RETURN
            ELSE IF(ICBNEW.EQ.ISIGN(I,KODE(KODIND-1))) THEN
C             Transmission from the fictitious part of block
              IEND=0
            ELSE
C             Termination of the ray computation in the fictitious part
C             of the block
              IEND=30
              RETURN
            END IF
   41     CONTINUE
        ELSE
C         Required complex block is not attainable
          IEND=22
        END IF
      END IF
      RETURN
C
C     Reflection-transmission code:
   50 CONTINUE
      KODIND=KODIND+1
      I=KODE(KODIND)
      IF(IABS(I).EQ.1) THEN
C       Transmission:
        IF(IY(8).GT.0) THEN
C         Transmission into material block
          ICBNEW=ISIGN(IY(8),I)
          IEND=0
        ELSE
C         Transmission into free space
          IEND=23
        END IF
      ELSE IF(IABS(I).EQ.2.AND.KODIND.GT.1) THEN
C       Reflection
        ICBNEW=ISIGN(IABS(IY(5)),I)
        IEND=0
      ELSE
C       End of code
        IEND=10
      END IF
      RETURN
      END
C
C=======================================================================
C