C Subroutine file 'rpar.for' to control the take-off parameters of rays.
C
C By Vlastislav Cerveny, Ludek Klimes, Ivan Psencik
C
C This file consists of the following subroutines:
C RPARB...Block data subroutine defining common block /RPARD/ to
C store the input data for the take-off parameters of rays,
C common block /RPARC/ to collect the quantities needed for
C boundary-value ray tracing, common block /RPARH/
C describing the histories of rays. The histories of rays
C are used by function IHIST to determine the ray-history
C index. Different values of ray-history index are attached
C to rays with different histories.
C RPAR1...Subroutine designed to read the input data for the
C take-off parameters, to prepare the parameters for the
C subroutine RPAR2 and to store them in the common block
C /RPARC/. It is called when starting the computation of a
C new elementary wave.
C RPAR2...Subroutine designed to specify the take-off parameters of
C individual rays.
C RPAR31..Subroutine called with constant step STORE of the
C independent variable along the ray, and at the points of
C intersection with interfaces either before and after the
C transformation, i.e. called simultaneously with the
C subroutine WRIT31 (see also C.R.T.5.5.1). It may be used
C to keep the quantities useful for the determination of the
C take-off parameters in the memory.
C RPAR32..Subroutine called at the points of intersection of the ray
C with specified interfaces, i.e. called simultaneously with
C the subroutine WRIT32 (see also C.R.T.5.5.2). It may be
C used to keep the quantities useful for the determination
C of the take-off parameters in the memory.
C RPAR33..Subroutine called at the endpoints of the elements of
C rays, situated at structural interfaces, i.e. called
C simultaneously with the subroutine WRIT33 (see also
C C.R.T.5.5.3). It may be used to keep the quantities
C useful for the determination of the take-off parameters in
C the memory.
C RPAR4...Subroutine called after finishing the computation of each
C ray. It defines four storage locations of the common
C block /INITC/ introduced in the file 'init.for'.
C XFUN... Subroutine returning the value and first derivatives of
C the X-functions describing the profiles. It is called by
C the subroutine RPAR32, and may be user-defined. In this
C version, it calls the subroutine SRFC2 from the subroutine
C file 'srfc.for', or directly returns one of the model
C coordinates.
C IHIST...Auxiliary integer function to the subroutine RPAR4 to
C determine the ray-history index. The ray-history index is
C used to distinguish the rays with different histories for
C the purposes of the subroutine RPAR2.
C SRPARH..Auxiliary subroutine to the subroutine RPAR31 to shift the
C storage locations of the common block /RPARH/.
C
C Note:
C The lines denoted by '*' in the first column contain an undebugged
C code intended for future extensions.
C
C Input data for the take-off parameters of rays:
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 TEXTP), the input parameter is of the
C type REAL.
C (1) TEXTP
C String describing the data. Only the first 80 characters of the
C string are significant.
C (2) ISRFR,ISRFX1,ISRFX2,NREC,XERR,AERR,PRM0(I),I=1,MPRM0,/
C ISRFR...Index of the reference surface at which the computed
C quantities are stored and the given profiles or receivers
C for two-point ray tracing are situated. The profiles are
C the crossections of the surface ISRFR and the surfaces
C defined by the X1-function. The receivers are the
C crossections of the surface ISRFR and the surfaces defined
C by the X1 and X2 functions. ISRFR must be one of the
C indices listed in the input data RAY-(7) (the indices for
C storing computed quantities), including the sign, see the
C subroutine file 'ray.for'. For ISRFR=0 there are no
C profiles for boundary-value ray tracing, ISRFR not listed
C in the input data RAY-(7) has the same effect.
C ISRFR=101,102,103,104,105, or 106 specify the boundary of
C the computational volume.
C Example: if the profile for the boundary-value ray tracing
C is situated at the model (earth's) surface, ISRFR is the
C index of the model surface. Its sign must correspond to
C the side of the surface ISRFR facing the material block.
C ISRFX1..Index of the function that becomes to be the X1-function.
C The X1-function is defined here by the subroutine XFUN.
C For ISRFX1=-1, -2, or -3, the X1-function coincides with
C the -ISRFX1-th model coordinate.
C ISRFX1=0: No boundary-value (two-point) ray tracing.
C The X1-function is zero and has no meaning.
C This option may be used only if the boundary-value rays
C are not searched for, i.e. if ISRFR.EQ.0, or NREC.EQ.0.,
C or IPOINT=0 at all input data lines (5.1).
C For ISRFX1 positive, the subroutine XFUN is an interface
C subroutine to the subroutine SRFC2, called with value
C ISRFX1 assigned to the dummy argument ISRFC, see the
C subroutine file 'srfc.for'. ISRFX1 must not exceed the
C number of the defined surfaces more than by 1.
C Example: The most often profile will be constructed as the
C intersection of the surface ISRFR with an auxiliary
C vertical plane. Using the basic version of the file
C 'srfc.for', the auxiliary vertical plane passing through
C the points (X1,X2,X3)=(X1A,X2A,X3A) and
C (X1,X2,X3)=(X1B,X2B,X3B) may be defined in the
C subroutine file 'srfc.for' by the following input data:
C For X1A.NE.X1B:
C 1 -2 0 0
C 2
C X1A X1B
C X2A X2B
C The interpolated function then returns the oriented
C distance from the vertical plane, measured in the
C direction of the X2-axis but with opposite sign.
C For X2A.NE.X2B:
C 2 -1 0 0
C 2
C X2A X2B
C X1A X1B
C The interpolated function then returns the oriented
C distance from the vertical plane, measured in the
C direction of the X1-axis but with opposite sign.
C ISRFX2..Index of the function that becomes to be the X2-function.
C ISRFX2=0: Initial-value ray traing or one-parametric
C shooting.
C Otherwise: Meaning of ISRFX2 is analogous to ISRFX1.
C NREC... Number of the defined receivers or profiles along the
C reference surface ISRFR for the two-point ray tracing.
C Zero if there are no profiles.
C Negative (e.g. -1) if the profile coordinates are to be
C read from a separate file, see (3.2). This is the
C recommended option.
C Positive values are left for the compatibility with the
C older versions, however, they are not recommended.
C XERR... Specifies the accuracy of the boundary-value ray tracing.
C If the distance of the ray from the given profile does not
C exceed XERR, the ray is considered to be a boundary-value
C ray. The distance of the ray from the given profile is
C measured by means of the difference of the X1-function.
C AERR... Specifies the accuracy of the determination of the
C boundary ray between the bundles of rays with different
C histories: the step between two consecutive rays with
C different histories is partitioned in order not to exceed
C the basic step in the first direction, multiplyed by AERR.
C No boundary rays are determined for AERR.GT.1 (e.g.
C AERR=999999.). Note that IPOINT=0 has the same effect as
C AERR.GT.1.
C PRM0(1)... Maximum alloved distance of the boundary ray from the
C shadow zone (measured on the reference surface).
C The distance is not measured if PRM0(1)=0 (default).
C Used only by two-parametric 3-D two-point ray tracing.
C PRM0(2)... Maximum alloved length of sides of homogeneous
C triangles (measured on the reference surface).
C The triangles are not measured if PRM0(2)=0 (default).
C Used only by two-parametric 3-D two-point ray tracing.
C PRM0(I),I=3,MPRM0... Reserved for future extension of the
C numerical parameters for two-point ray tracing.
C Default values:
C ISRFR=0, ISRFX1=0, ISRFX2=0, NREC=0, XERR=0.000, AERR=999999,
C PRM0(I)=0,I=1,MPRM0.
C The boundary rays are not determined if AERR.GT.1.
C Neither are they if IPOINT=0, see input data (5.1).
C The boundary-value rays are not determined if ISRFR.EQ.0
C or NREC.EQ.0.
C Neither are they if IPOINT=0, see input data (5.1).
C Thus, just the initial-value ray tracing is performed if
C (ISRFR.EQ.0 or NREC.LE.0) and AERR.GT.1.
C Just the initial-value ray tracing is also performed if
C IPOINT=0, see input data (5.1).
C During the initial-value ray tracing, only the basic net
C of rays with the given constant steps is generated.
C (3) Specification of the profiles for the one-parametric
C boundary-value ray tracing.
C This input is not performed if NREC=0.
C For NREC positive and ISRFX2.EQ.0 (3.1),
C For NREC positive and ISRFX2.NE.0 (3.2), for NREC negative (3.3):
C (3.1) XREC1(1),...,XREC1(NREC)
C XREC1(I)... The I-th profile is the intersection of the surface
C ISRFR with the isosurface of the X1-function,
C corresponding to the value XREC1(I). The values XREC1(I)
C may be specified in any order.
C (3.1) XREC1(1),XREC2(1),...,XREC1(NREC),XREC2(NREC)
C XREC1(I),XREC2(I)... The I-th receiver is the intersection of the
C surface ISRFR with the isosurface of the X1-function,
C corresponding to the value XREC1(I), and the isosurface of
C the X2-function, corresponding to the value XREC2(I).
C (3.3) 'FREC'
C 'FREC'..Name of the file with profile or receiver coordinates.
C Default: 'FREC'='rec.dat'.
C (4) The data specifying the functions No. ISRFX1 and ISRFX2. This
C auxiliary function is specified by the subroutine SRFC2 and the
C data are read by subroutine SRFC1. For the description of the
C data refer to the subroutine SRFC1 of the file 'srfc.for'. This
C input is performed only if ISRFX1 or ISRFX2 is positive and the
C function ISRFX1 or ISRFX2 is not defined in the data set
C 'model.dat' or 'dcrt.dat'. In such a case, the index of the first
C undefined function of ISRFX1 or ISRFX2 must equal the number NSRFC
C of the surfaces defined in the data set 'model.dat' + the number
C ISRFCA of the auxiliary surfaces defined in the data set
C 'dcrt.dat' + 1, and so on.
C (5) For each elementary wave IWAVE the following data (5.1) and (5.2):
C (5.1) Any times:
C PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM,IPOINT,
C PRM1(I),I=1,MPRM1,/
C PAR1L,PAR2L... Ray parameters corresponding to the origin (0,0) of
C the domain of the normalized ray parameters.
C PAR1A,PAR2A... Ray parameters corresponding to the right bottom
C corner (1,0) of the domain of the normalized ray
C parameters.
C PAR1B,PAR2B... Ray parameters corresponding to the left top corner
C (0,1) of the domain of the normalized ray parameters.
C ANUM,BNUM... Specify the 2*2 metric tensor G in the normalized ray
C parameter domain.
C ANUM,BNUM in RP2D (Initial-value ray tracing or one-parametric
C shooting):
C The components of the metric tensor are
C G11=ANUM*ANUM, G12=0., G22=BNUM*BNUM. The system of basic
C rays forms a regular rectangular grid in the normalized
C ray parameter domain. The grid cells are unit squares
C with respect to the above metric:
C For ANUM.GT.0 and BNUM.GT.0: the two-parametric basic net
C of INT(ANUM)*INT(ANUM) rays is given by the take-off ray
C parameters
C PAR1=PAR1L+A*(PAR1A-PAR1L)/ANUM+B*(PAR1B-PAR1L)/BNUM,
C PAR2=PAR2L+A*(PAR2A-PAR2L)/ANUM+B*(PAR2B-PAR2L)/BNUM,
C where
C A=0,1,2,...,INT(ANUM) and B=0,1,2,...,INT(BNUM).
C The basic net of rays corresponds to the initial-value
C ray tracing.
C Besides the basic net of rays, also additional rays are
C computed in order to find the boundary rays between the
C bundles of rays with different histories, and to find
C the boundary-value rays. All additional rays are given
C by
C B=0,1,2,...,INT(BNUM)
C and real a from the interval 0 to INT(ANUM), forming
C the cuts of the ray-parameter surface given by fixed B
C and parametrized by A. The boundary rays and the
C boundary-value rays are searched for along the mentioned
C cuts. Under the boundary-value rays we understand here
C the rays from the above mentioned cuts passing through
C the XERR-vicinity of the given profiles.
C For ANUM.GT.0 and BNUM.EQ.0: The one-parametric basic net
C of INT(ANUM) rays is given by the take-off ray
C parameters
C PAR1=PAR1L+A*(PAR1A-PAR1L)/ANUM,
C PAR2=PAR2L+A*(PAR2A-PAR2L)/ANUM,
C where
C A=0,1,2,...,INT(ANUM).
C Additional rays are computed as for ANUM.GT.0 and
C BNUM.GT.0.
C For ANUM.EQ.0 and BNUM.GT.0: Similarly for ANUM.GT.0 and
C BNUM.EQ.0, but no additional rays are computed.
C For ANUM.EQ.0 and BNUM.EQ.0: Just a single ray given by
C PAR1=PAR1L and PAR2=PAR2L is computed.
C ANUM,BNUM in RP3D (two-parametric shooting):
C The components of the metric tensor are
C G11=ANUM*ANUM*E11, G12=ANUM*BNUM*E12, G22=BNUM*BNUM*E22,
C where E is some basic metric tensor, e.g. the unit matrix
C or the Cartesian metric tensor on the unit sphere around
C the source. The system of basic rays forms basic
C triangles covering the normalized ray parameter domain.
C The basic triangles are constructed in a way to be roughly
C equilateral with unit sides with respect to metric G.
C IPOINT..IPOINT=0: No ray history is recorded, all rays have the
C same index 0 and no boundary-value ray tracing is
C performed.
C Otherwise: The IABS(IPOINT)-th point of intersection of
C a ray with the surface ISRFR is considered for the
C boundary-value ray tracing.
C IPOINT negative: Only those rays satisfying the whole
C code (IEND=10) are consideredred to be successful. All
C other rays are rejected (IMPORTANT OPTION).
C IPOINT=large: Ray history is recorded, boundary rays may
C be determined but no boundary-value ray tracing is
C performed because all rays are unsuccessful.
C In most cases, user will wish the value of IPOINT=1.
C For the structural interface (not auxiliary surface), the
C positive or negative side of the surface ISRFR is
C specified, then:
C No point of intersection is accounted if the reflected
C ray is situated at the other side of the surface ISRFR
C than the specified one, and
C two points of intersection are accounted if the
C reflected ray is situated at the specified side of the
C surface ISRFR.
C Just the initial-value ray tracing is performed if
C IPOINT=0.
C PRM1(I),I=1,MPRM1... Reserved for future extension of the
C numerical parameters for two-point ray tracing, specific
C to the calculated elementary wave.
C Default values:
C PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM,IPOINT:
C Zeroes when starting the complete ray tracing program,
C otherwise the previous values.
C PRM1(I)=0, I=1,MPRM1.
C (5.2) /(a slash)
C
C Input file 'FREC' containing profile (or receiver) coordinates:
C This file is used only if NREC.LT.0, see the input data (2) above.
C (1) Several strings terminated by / (a slash).
C The simplest way is to submit just the /.
C (2) Several times (2.1):
C (2.1) 'NAMER(IR)',XR(1,IR),XR(2,IR),XR(3,IR),XR(4,IR),XR(5,IR),/
C 'NAMER(IR)'... String reserved for the name of the receiver.
C No meaning here, anything in apostrophes may be submitted.
C XR(1,IR),XR(2,IR),XR(3,IR)... Model coordinates of the receiver.
C For ISRFX1.LT.0: XREC1(IR)=XR(-ISRFX1,IR) is the first
C coordinate of the point along the reference surface.
C For ISRFX2.LT.0: XREC2(IR)=XR(-ISRFX2,IR) is the second
C coordinate of the point along the reference surface.
C Values of other model coordinates have no meaning here.
C XR(4,IR),XR(5,IR)... Coordinates of the receiver along the
C reference surface other than the model coordinates.
C For ISRFX1.GT.0: XREC1(IR)=XR(4,IR) is the first
C coordinate of the point along the reference surface.
C For ISRFX2.GT.0: XREC2(IR)=XR(5,IR) is the second
C coordinate of the point along the reference surface.
C Other values have no meaning here and need not be
C specified.
C Default: X1R(IR)=0, X2R(IR)=0, X3R(IR)=0, X4R(IR)=0, X5R(IR)=0.
C (3) / (a slash).
C
C Storage in the memory:
C The input data (2) are stored in common block /RPARD/ defined
C together with common blocks /RPARC/ and /RPARH/ in the following
C subroutine:
C ------------------------------------------------------------------
BLOCK DATA RPARB
INCLUDE 'rpard.inc'
INCLUDE 'rparc.inc'
END
C ------------------------------------------------------------------
C
C Date: 1996, September 30
C Coded by Ludek Klimes
C
C=======================================================================
C
SUBROUTINE RPAR1(LUN,IWAVE)
INTEGER LUN,IWAVE
C
C This subroutine is called when starting the computation of a new
C elementary wave. It prepares the parameters for the subroutine RPAR2
C to specify the take-off parameters of the individual rays.
C
C Input:
C LUN... Logical unit number of the external input device
C containing the input data.
C IWAVE...Index of the elementary wave that is going to be computed,
C i.e. the output from the last invocation of the subroutine
C CODE1.
C
C No output.
C
C Common block /DCRT/:
INCLUDE 'dcrt.inc'
C None of the storage locations of the common block are altered.
C
C Common block /RPARD/:
INCLUDE 'rpard.inc'
C Nearly all the storage locations of the common block are defined in
C this subroutine.
C
C Common block /RPARC/:
INCLUDE 'rparc.inc'
C None of the storage locations of the common block, except LURPAR, are
C altered. LURPAR is set to LUN.
C
C Subroutines and external functions required:
INTEGER NSRFC
EXTERNAL RPARB,NSRFC
C RPARB.. Block data subroutine of this file.
C NSRFC... File 'model.for' of the package 'MODEL'.
C
C Error messages:
C 650... NREC negative:
C Negative input value of NREC, data line (2), is not
C allowed.
C 651... Too small array XREC in /RPARD/:
C The number of given profiles (specified by the values of
C X-functions, see input data 3) is exceeding the dimension
C mrec of the array XREC in the common block /RPARD/. Thus,
C mrec should be adjusted in all declarations of /RPARD/.
C 652... Wrong value of ISRFX1:
C The function no. ISRFX1 specified by the routines of the
C file 'srfc.for' is defined neither in the input data set
C 'model.dat' nor 'dcrt.dat' and cannot be defined in this
C data set. The index ISRFX1 is too large or is less than
C -3.
C
C Date: 1996, September 30
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Auxiliary storage locations:
INTEGER LUREC
PARAMETER (LUREC=9)
CHARACTER*80 TEXTP
INTEGER I,J
REAL AUX(5)
C
C LUREC...Logical unit number of the receiver file. The file is
C opened and closed here.
C TEXTP...The name of the data. String of 80 characters.
C I,J... Auxiliary loop variables, or temporary variables.
C AUX... Temporary storage locations to read the receiver
C coordinates.
C
C.......................................................................
C
IF(IWAVE.EQ.1) THEN
C (1) Reading the name of the input data:
READ(LUN,*) TEXTP
C (2):
ISRFR=0
ISRFX(1)=0
ISRFX(2)=0
NREC=0
XERR=0.000
AERR=999999.
DO 10 I=1,MPRM0
PRM0(I)=0.
10 CONTINUE
PRM0(2)=999999.
READ(LUN,*) ISRFR,ISRFX(1),ISRFX(2),NREC,XERR,AERR,PRM0
c IF(NREC.LT.0) THEN
c PAUSE 'Error 650 in RPAR1: NREC negative'
c STOP
c END IF
C (3):
IF(NREC.GT.0) THEN
IF(NREC.GT.MREC) THEN
PAUSE 'Error 651 in RPAR1: Too small array XREC in /RPARD/'
STOP
END IF
IF(ISRFX(2).EQ.0) THEN
C Initial-value or 1-paremetric shooting:
READ(LUN,*) (XREC(1,I),I=1,NREC)
ELSE
C 2-paremetric shooting:
READ(LUN,*) (XREC(1,I),XREC(2,I),I=1,NREC)
END IF
ELSE IF(NREC.LT.0) THEN
TEXTP='rec.dat'
READ(LUN,*) TEXTP
OPEN(LUREC,FILE=TEXTP,STATUS='OLD')
READ(LUREC,*) (TEXTP,I=1,20)
DO 11 I=1,MREC+1
TEXTP='$'
AUX(1)=0.
AUX(2)=0.
AUX(3)=0.
AUX(4)=0.
AUX(5)=0.
READ(LUREC,*) TEXTP,(AUX(J),J=1,5)
IF(TEXTP.EQ.'$') THEN
GO TO 12
END IF
IF(I.GT.MREC) THEN
PAUSE'Error 651 in RPAR1: Too small array XREC in /RPARD/'
STOP
END IF
IF(ISRFX(1).LT.0) THEN
XREC(1,I)=AUX(-ISRFX(1))
ELSE IF(ISRFX(1).GT.0) THEN
XREC(1,I)=AUX(4)
END IF
IF(ISRFX(2).LT.0) THEN
XREC(2,I)=AUX(-ISRFX(2))
ELSE IF(ISRFX(2).GT.0) THEN
XREC(2,I)=AUX(5)
END IF
11 CONTINUE
12 CONTINUE
NREC=I-1
CLOSE(LUREC)
END IF
C (4) Auxiliary functions no. ISRFX(1) and ISRFX(2):
I=NSRFC()+NSRFCA+1
IF(ISRFX(1).EQ.I) THEN
CALL SRFC1(LUN,1)
I=I+1
ELSE IF(ISRFX(1).LT.-3.OR.ISRFX(1).GT.I) THEN
PAUSE 'Error 652 in RPAR1: Wrong value of ISRFX1'
STOP
END IF
IF(ISRFX(2).EQ.I) THEN
CALL SRFC1(LUN,1)
ELSE IF(ISRFX(2).LT.-3.OR.ISRFX(2).GT.I) THEN
PAUSE 'Error 653 in RPAR1: Wrong value of ISRFX2'
STOP
END IF
C
C Default values
PAR1L=0.
PAR2L=0.
PAR1A=0.
PAR2A=0.
PAR1B=0.
PAR2B=0.
ANUM=0.
BNUM=0.
IPOINT=0
END IF
C
LURPAR=LUN
RETURN
END
C
C=======================================================================
C
SUBROUTINE RPAR2(IRAY,PAR1,PAR2)
REAL PAR1,PAR2
INTEGER IRAY
C
C This subroutine determines the take-off parameters of the ray.
C
C Input:
C IRAY... Number of the already computed rays. IRAY=0 when
C beginning the computation of a new elementary wave.
C Otherwise, the output from the previous invocation of
C RPAR2.
C
C Output:
C IRAY... IRAY=0 if all rays are computed and the computation of the
C elementary wave will be terminated.
C Otherwise, input value increased by 1.
C PAR1,PAR2... Take-off parameters of the ray.
C
C Common block /RPARD/:
INCLUDE 'rpard.inc'
C Storage locations PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM,
C IPOINT,PRM1 of the common block may be redefined in this subroutine.
C
C Common blocks /RPARC/ and /RPARH/:
INCLUDE 'rparc.inc'
C Storage locations IRAY0,JRAY,JTYPE,G1,G2 of common block /RPARC/ are
C altered.
C Storage locations JPOINT and KMAH of common block /RPARH/ are
C annulled.
C Storage locations NRPARH and IRPARH(0) of common block /RPARH/ may be
C altered.
C
C Subroutines and external functions required:
EXTERNAL RP2D
C RP2D... One-parametric boundary-value ray tracing. This file.
C
C Error messages:
C 654... Zero value of ISRFX1:
C ISRFX1 must not be zero if the boundary-value ray tracing
C is to be performed, i.e. if ISRFR.NE.0 and NREC.NE.0 and
C IPOINT.NE.0 in the input data.
C
C Date: 1996, June 15
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Auxiliary storage locations:
INTEGER I,IAUX
REAL AUX1,AUX2,AUX3,AUX4,AUX5,AUX6,AUX7,AUX8
C
C.......................................................................
C
JPOINT=0
KMAH=0
C
IF(IRAY.EQ.0) THEN
C New elementary wave
IRAY0=0
IF(ISRFX(2).EQ.0) THEN
G2=0.
END IF
GO TO 90
END IF
C
C New ray - determination of shooting parameters G1,G2:
10 CONTINUE
JRAY=IRAY-IRAY0
IF(ISRFX(2).EQ.0) THEN
C One-parametric shooting (or initial-value ray tracing):
CALL RP2D(IRAY0,JRAY,G1)
IF(JRAY.EQ.0) THEN
C End of shooting - new one-parametric slice or new data
G2=G2+1.
IF(G2.GT.BNUM+0.001) THEN
C New data or new elementary wave
GO TO 90
ELSE
C New one-parametric slice
IRAY0=IRAY
GO TO 10
END IF
END IF
ELSE
C Two-parametric shooting:
CALL RP3D(JRAY,JTYPE,G1,G2)
IF(JRAY.EQ.0) THEN
C End of shooting - new data or new elementary wave
GO TO 90
END IF
END IF
IRAY=IRAY0+JRAY
C
C New ray - determination of take-off parameters PAR1,PAR2:
NRPARH=IRPARH(IRPARH(0))
IF(ISRFX(2).EQ.0) THEN
C One-parametric shooting:
C Projection of the shooting domain (0,ANUM)*(0,BNUM)
PAR1=PAR1L
PAR2=PAR2L
IF(ANUM.NE.0.) THEN
PAR1=PAR1+(PAR1A-PAR1L)*G1/ANUM
PAR2=PAR2+(PAR2A-PAR2L)*G1/ANUM
END IF
IF(BNUM.NE.0.) THEN
PAR1=PAR1+(PAR1B-PAR1L)*G2/BNUM
PAR2=PAR2+(PAR2B-PAR2L)*G2/BNUM
END IF
ELSE
C Two-parametric shooting:
C Projection of the shooting domain (0,1)*(0,1)
PAR1=PAR1L+(PAR1A-PAR1L)*G1+(PAR1B-PAR1L)*G2
PAR2=PAR2L+(PAR2A-PAR2L)*G1+(PAR2B-PAR2L)*G2
END IF
RETURN
C
C Reading input data
90 CONTINUE
AUX1=PAR1L
AUX2=PAR2L
AUX3=PAR1A
AUX4=PAR2A
AUX5=PAR1B
AUX6=PAR2B
AUX7=ANUM
AUX8=BNUM
IAUX=IPOINT
DO 91 I=1,MPRM1
PRM1(I)=0.
91 CONTINUE
READ(LURPAR,*)
* PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM,IPOINT,PRM1
IF(IRAY.NE.0) THEN
IF(PAR1L.EQ.AUX1.AND.PAR2L.EQ.AUX2.AND.PAR1A.EQ.AUX3.AND.
* PAR2A.EQ.AUX4.AND.PAR1B.EQ.AUX5.AND.PAR2B.EQ.AUX6.AND.
* ANUM.EQ.AUX7.AND.BNUM.EQ.AUX8.AND.IAUX.EQ.IPOINT) THEN
C New elementary wave
IRAY=0
RETURN
END IF
END IF
IF(ISRFX(1).EQ.0.AND.ISRFR.NE.0.AND.NREC.NE.0.AND.IPOINT.NE.0)
* THEN
PAUSE 'Error 654 in RPAR2: Zero value of ISRFX1'
STOP
END IF
C Initializing ray histories
INDMIN=0
IRPARH(0)=0
C New initial ray
IRAY0=IRAY
GO TO 10
C
END
C
C=======================================================================
C
SUBROUTINE RPAR31(YL,Y,YY,IY)
REAL YL(6),Y(35),YY(5)
INTEGER IY(12)
C
C This subroutine may be used to store the quantities useful for the
C determination of the take-off parameters in the memory. It is called
C with constant step STORE (see the input data in the subroutine file
C 'ray.for') of the independent variable along the ray (if STORE.NE.0),
C and at the points of intersection with interfaces either before and
C after the transformation (in any case).
C
C Input:
C YL... Array containing local quantities at the point of the ray.
C Y... Array containing basic quantities computed along the ray.
C YY... Array containing real auxiliary quantities computed along
C the ray.
C IY... Array containing integer auxiliary quantities computed
C along the ray.
C None of the input parameters are altered.
C
C No output.
C
C Common block /RPARD/:
INCLUDE 'rpard.inc'
C None of the storage locations of the common block are altered.
C
C Common block /RPARH/:
INCLUDE 'rparc.inc'
C Storage locations of the common block may be altered.
C
C Subroutines and external functions required:
EXTERNAL NSRFC
INTEGER NSRFC
C NSRFC... File 'model.for' of the package 'MODEL'.
C
C Date: 1996, September 29
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C No auxiliary storage locations.
C
IF(IPOINT.NE.0) THEN
IF(NRPARH.EQ.IRPARH(IRPARH(0))) THEN
C Initial element of the ray
NRPARH=NRPARH+1
CALL SRPARH
IRPARH(NRPARH)=IY(5)
END IF
IF(JPOINT.LT.IABS(IPOINT)) THEN
C4.10 IF(IY(6).NE.0.AND.IABS(IY(6)).LE.NSRFC()) THEN
IF(IY(6).NE.0) THEN
C Structural interface
NRPARH=NRPARH+1
CALL SRPARH
IF(MOD(NRPARH-IRPARH(IRPARH(0)),2).EQ.0) THEN
C Incidence at a structural interface
IRPARH(NRPARH)=IY(6)
ELSE
C Leaving a structural interface
IRPARH(NRPARH)=IY(5)
END IF
END IF
END IF
END IF
RETURN
END
C
C=======================================================================
C
SUBROUTINE RPAR32(ISTOR,YL,Y,YY,IY)
INTEGER ISTOR,IY(12)
REAL YL(6),Y(27),YY(5)
C
C This subroutine may be used to store the quantities useful for the
C determination of the take-off parameters in the memory. It is called
C at the points of intersection of the ray with interfaces specified in
C the input data in the subroutine file 'ray.for' for storing the
C computed quantities (see also C.R.T.5.5.2).
C
C Input:
C ISTOR...The sequential number in the iput data of the specified
C surface.
C YL... Array containing local quantities at the point of the ray.
C Y... Array containing basic quantities computed along the ray.
C YY... Array containing real auxiliary quantities computed along
C the ray.
C IY... Array containing integer auxiliary quantities computed
C along the ray.
C None of the input parameters are altered.
C
C No output.
C
C Common block /DCRT/ (subroutine file 'ray.for'):
INCLUDE 'dcrt.inc'
C None of the storage locations of the common block are altered.
C
C Common block /INITC/ (see subroutine file 'init.for'):
INCLUDE 'initc.inc'
C None of the storage locations of the common block are altered.
C
C Common block /RPARD/:
INCLUDE 'rpard.inc'
C None of the storage locations of the common block are altered.
C
C Common block /RPARC/ and /RPARH/:
INCLUDE 'rparc.inc'
C Storage locations X1,X2,X1G1,X2G1,X1G2,X2G2 of common block /RPARC/
C may be altered.
C Storage locations JPOINT and KMAH of common block /RPARH/ may be
C altered.
C
C Subroutines and external functions required:
EXTERNAL METRIC,SMVPRD,XFUN
C METRIC... File 'metric.for' of the package 'MODEL'.
C SMVPRD... File 'means.for' of the package 'MODEL'.
C XFUN... This file.
C
C Date: 1996, Septenber 30
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Substitution in the common block /RPARC/:
REAL X(2),XG(4)
EQUIVALENCE (X(1),X1),(XG(1),X1G1)
C
C Auxiliary storage locations:
INTEGER NX,IX
REAL GSQRD,G(12),AUX(18),AUX1,AUX2,AUX3
EQUIVALENCE (AUX(11),AUX1),(AUX(12),AUX2),(AUX(13),AUX3)
REAL H42,H52,H62,RH1,RH2,RR,XH1,XH2,XR,SIDE1,SIDE2
C
C NX... Number of defined X-functions.
C IX... Auxiliary loop variable
C GSQRD,G... See subroutine METRIC of the file 'metric.for'.
C AUX,AUX1,AUX2,AUX3... Auxiliary storage locations:
C AUX(1:18)... Christoffel symbols,
C AUX(1:4)... Value and gradient of the function specifying
C the reference surface ISRFR,
C AUX(5:10)... Second derivatives of the function specifying
C the reference surface ISRFR,
C AUX(5:7)... Gradient of the X-function.
C AUX(11:13)... AUX1,AUX2,AUX3,
C AUX1,AUX2,AUX3... Contravariant derivatives the functions
C etc.
C H42,H52,H62... Contravariant components of the 2-nd basis vector
C of R.C.C.S.
C RH1,RH2... Derivatives of the function describing reference the
C surface ISRFR in R.C.C.S.
C RR... Norm of the gradient of the function describing reference
C the surface ISRFR in R.C.C.S.
C XH1,XH2... Derivatives of the X-functions in R.C.C.S.
C XR... Projection of the gradient of the X-functions onto the
C normal to the reference surface.
C SIDE1,SIDE2... Vectorial side of the shooting domain, in the ray
C parameters PAR1,PAR2.
C
C.......................................................................
C
IF(KSTOR(ISTOR).EQ.ISRFR) THEN
DO 1 IX=1,ISTOR-1
IF(KSTOR(IX).EQ.ISRFR) THEN
C The surface is already encountered
GO TO 9
END IF
1 CONTINUE
JPOINT=JPOINT+1
IF(JPOINT.EQ.IABS(IPOINT)) THEN
C Succesful ray:
C
C KMAH index:
KMAH=IY(12)
C
C Metric tenzor
CALL METRIC(Y(3),GSQRD,G,AUX)
C
C Contravariant components of the 2-nd basis vector of R.C.C.S.:
CALL SMVPRD(G(7),Y(6),Y(7),Y(8),AUX1,AUX2,AUX3)
AUX1=GSQRD*SQRT(AUX1*Y(6)+AUX2*Y(7)+AUX3*Y(8))
H42=(Y(7)*Y(11)-Y(8)*Y(10))/AUX1
H52=(Y(8)*Y( 9)-Y(6)*Y(11))/AUX1
H62=(Y(6)*Y(10)-Y(7)*Y( 9))/AUX1
C
C Gradient of the function describing reference surface ISRFR:
C Covariant components
IF(IABS(ISRFR).LE.100) THEN
CALL SRFC2(IABS(ISRFR),Y(3),AUX)
ELSE
C Reference surface limits the computational volume
AUX(2)=0.
AUX(3)=0.
AUX(4)=0.
AUX((IABS(ISRFR)-97)/2)=1.
END IF
C Contravariant components
CALL SMVPRD(G(7),AUX(2),AUX(3),AUX(4),AUX1,AUX2,AUX3)
C Gradient in R.C.C.S.
RH1=AUX1*Y(9)+AUX2*Y(10)+AUX3*Y(11)
RH2=AUX(2)*H42+AUX(3)*H52+AUX(4)*H62
C Projection onto the normal to the reference surface
RR =AUX1*AUX(2)+AUX2*AUX(3)+AUX3*AUX(4)
C
C Number of X-functions:
IF(ISRFX(2).EQ.0) THEN
IF(ISRFX(1).EQ.0) THEN
NX=0
ELSE
NX=1
END IF
ELSE
NX=2
END IF
C
C Loop for X-functions:
DO 5 IX=1,NX
C
C Value and gradient of the X-functions:
C Covariant components
CALL XFUN(IX,Y(3),X(IX),AUX(5))
C contravariant components
CALL SMVPRD(G(7),AUX(5),AUX(6),AUX(7),AUX1,AUX2,AUX3)
C Gradient in R.C.C.S.
XH1=AUX1*Y(9)+AUX2*Y(10)+AUX3*Y(11)
XH2=AUX(5)*H42+AUX(6)*H52+AUX(7)*H62
C Projection onto the normal to the reference surface
XR =AUX1*AUX(2)+AUX2*AUX(3)+AUX3*AUX(4)
C
C Derivatives of the X-function along the reference surface
C in R.C.C.S.:
AUX1=XH1*RH2-XH2*RH1
AUX2=1.-(RH1*RH1+RH2*RH2)/RR
XH1=(XH1-(XR*RH1+AUX1*RH2)/RR)/AUX2
XH2=(XH2-(XR*RH2-AUX1*RH1)/RR)/AUX2
C
C Derivative of the X-function with respect to the shooting
C parameters:
SIDE1=PAR1A-PAR1L
SIDE2=PAR2A-PAR2L
XG(IX)= ((XH1*Y(12)+XH2*Y(13))*(YI(12)*SIDE1+YI(16)*SIDE2)
* +(XH1*Y(16)+XH2*Y(17))*(YI(13)*SIDE1+YI(17)*SIDE2)
* +(XH1*Y(20)+XH2*Y(21))*(YI(14)*SIDE1+YI(18)*SIDE2)
* +(XH1*Y(24)+XH2*Y(25))*(YI(15)*SIDE1+YI(19)*SIDE2))
SIDE1=PAR1B-PAR1L
SIDE2=PAR2B-PAR2L
XG(IX+2)=((XH1*Y(12)+XH2*Y(13))*(YI(12)*SIDE1+YI(16)*SIDE2)
* +(XH1*Y(16)+XH2*Y(17))*(YI(13)*SIDE1+YI(17)*SIDE2)
* +(XH1*Y(20)+XH2*Y(21))*(YI(14)*SIDE1+YI(18)*SIDE2)
* +(XH1*Y(24)+XH2*Y(25))*(YI(15)*SIDE1+YI(19)*SIDE2))
C
5 CONTINUE
C
C Undefined functional values:
DO 6 IX=NX+1,2
X(IX)=0.
XG(IX)=0.
XG(IX+2)=0.
6 CONTINUE
C
END IF
END IF
9 CONTINUE
RETURN
END
C
C=======================================================================
C
SUBROUTINE RPAR33(YL,Y,YY,IY)
REAL YL(6),Y(35),YY(5)
INTEGER IY(12)
C
C This subroutine may be used to store the quantities useful for the
C determination of the take-off parameters in the memory. It is called
C at the endpoints of the elements of rays, situated at structural
C interfaces (see also C.R.T.5.5.3).
C
C Input:
C YL... Array containing local quantities at the point of the ray.
C Y... Array containing basic quantities computed along the ray.
C YY... Array containing real auxiliary quantities computed along
C the ray.
C IY... Array containing integer auxiliary quantities computed
C along the ray. For the definition of arrays YL, Y, YY and
C IY see the file 'ray.for'.
C None of the input parameters are altered.
C
C No output.
C
C No subroutines and external functions required.
C
C Date: 1990, January 22
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
RETURN
END
C
C=======================================================================
C
SUBROUTINE RPAR4(IRAY,PAR1,PAR2,YL,Y,YY,IY,IEND,ISHEET,IREC)
C
REAL PAR1,PAR2,YL(6),Y(35),YY(5)
INTEGER IRAY,IY(12),IEND,ISHEET,IREC
C
C This subroutine is called after finishing the computation of a ray.
C It defines the storage locations YI(22) to YI(25) of the common block
C /INITC/.
C
C Input:
C IRAY... Index of the ray, i.e. the number of the already computed
C rays, i.e. the output from the last invocation of RPAR2.
C PAR1,PAR2... Take-off parameters of the ray.
C YL... Array containing local quantities at the point of the ray.
C Y... Array containing basic quantities computed along the ray.
C YY... Array containing real auxiliary quantities computed along
C the ray.
C IY... Array containing integer auxiliary quantities computed
C along the ray. For the definition of arrays YL, Y, YY and
C IY see the file 'ray.for'.
C IEND... Reason of the termination of the computation of the ray
C (see C.R.T.5.4). For a detailed description see
C subroutine RAY (subroutine file 'ray.for').
C
C Output:
C ISHEET..Ray-history index. The different ray histories are
c consecutively indexed by positive integers 1,2,3,...
C according to their appearance during ray tracing.
C The ray histories are indexed independently within each
C elementary wave.
C The ray-history indices are complemented with sign:
C Positive - successful ray (crossing reference surface),
C negative - unsuccessful ray (terminating before crossing
C reference surface).
C IREC... Index of the receiver for a two-point ray,
C otherwise zero.
C
C Common block /INITC/ (see subroutine file 'init.for'):
INCLUDE 'initc.inc'
C ------------------------------------------------------------------
C YI(22)..Area of the element of the ray-parameter surface,
C corresponding to the ray,
C YI(23),YI(24),YI(25)... Components 11, 12, 22 of the symmetric
C matrix inverse to the specific moment of the element of
C the ray-parameter surface corresponding to the ray, see
C C.R.T.6.1.
C Locations YI(22) to YI(25) of the common block are defined in this
C subroutine, other storage locations remain unchanged.
C
C Common block /RPARD/:
INCLUDE 'rpard.inc'
C None of the storage locations of the common block are altered.
C
C Common block /RPARC/:
INCLUDE 'rparc.inc'
C Storage locations of the common block are not altered.
C
C Subroutines and external functions required:
EXTERNAL IHIST
INTEGER IHIST
C IHIST...This file.
C
C Date: 1996, June 17
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Auxiliary storage locations:
REAL X1A
REAL RG11,RG12,RG22,SG11,SG12,SG22,SIDE1,SIDE2,AUX1,AUX2
REAL G1X1,G2X1,G1X2,G2X2
C
C X1A... Derivative of the X1-function with respect to the
C shooting parameter A.
C RG11,RG12,RG22... Ray take-off parameter metric tensor.
C SG11,SG12,SG22... Shooting parameter metric tensor.
C SIDE1,SIDE2,AUX1,AUX2... Temporary variables.
C G1X1,G1X2,G2X1,G2X2... Derivatives of shooting parameters with
C respect to the X-functions.
C
C.......................................................................
C
C Ray-history index: positive - successful ray,
C otherwise - unsuccessful ray
ISHEET=IHIST(IEND)
C
C Calculating shooting-parameter metric tensor, storing the results:
IF(ISRFX(2).EQ.0) THEN
C Storing the quantities for one-parametric shooting:
IF(ANUM.GT.0.) THEN
X1A=X1G1/ANUM
ELSE
X1A=0.
END IF
CALL RP2DM(JRAY,ISHEET,X1,X1A,IREC)
C X1,X1G1... Determined in RPAR32 (if defined),
C ISHEET... Determined in this subroutine.
C IREC... Determined in RP2DM.
ELSE
C (a) metric tensor per ray take-off parameters:
C Check for non-existing rays:
IF(IEND.LT.70) THEN
C Cartesian metric tensor (zero for a point source)
RG11= YI(12)*YI(12)+YI(13)*YI(13)
RG12= YI(12)*YI(16)+YI(13)*YI(17)
RG22= YI(16)*YI(16)+YI(17)*YI(17)
IF(ABS(YI(12)*YI(17)-YI(13)*YI(14)).LT.0.000001) THEN
C Plus unit-sphere metric tensor for point or line source
AUX1=YI(6)**2+YI(7)**2+YI(8)**2
RG11=RG11+(YI(14)*YI(14)+YI(15)*YI(15))/AUX1
RG12=RG12+(YI(14)*YI(18)+YI(15)*YI(19))/AUX1
RG22=RG22+(YI(18)*YI(18)+YI(19)*YI(19))/AUX1
END IF
C The above two metric tensors (equal along an exploding unit
C sphere source) should be combined in a more sophisticated way.
C **************************************************************
C The metric tensor should be determined in 'init.for' !
C **************************************************************
ELSE
RG11=1.
RG12=0.
RG22=1.
END IF
C (b) Metric tensor per shooting parameters:
C Scale renormalizes the unit-sphere metric in radians to the
C metric in some scaled units, proportional to radians by the
C factors of scale in the directons of both the shooting
C parameters.
SIDE1=PAR1A-PAR1L
SIDE2=PAR2A-PAR2L
SCALE=ANUM/SQRT(SIDE1*SIDE1+SIDE2*SIDE2)
SIDE1=SIDE1*SCALE
SIDE2=SIDE2*SCALE
AUX1=RG11*SIDE1+RG12*SIDE2
AUX2=RG12*SIDE1+RG22*SIDE2
SG11=SIDE1*AUX1+SIDE2*AUX2
SIDE1=PAR1B-PAR1L
SIDE2=PAR2B-PAR2L
SCALE=BNUM/SQRT(SIDE1*SIDE1+SIDE2*SIDE2)
SIDE1=SIDE1*SCALE
SIDE2=SIDE2*SCALE
SG12=SIDE1*AUX1+SIDE2*AUX2
AUX1=RG11*SIDE1+RG12*SIDE2
AUX2=RG12*SIDE1+RG22*SIDE2
SG22=SIDE1*AUX1+SIDE2*AUX2
C Another extremely simple metric, conformal with the metric used
C at one-parametric shooting, and coinciding with the above scaled
C metric for two-parametric shooting in the case of a small
C rectangular shooting domain situated at the equator:
C SG11=ANUM*ANUM
C SG12=0.
C SG22=BNUM*BNUM
C Regularization near the poles:
SG11=SG11+ANUM*ANUM/100.
SG22=SG22+BNUM*BNUM/100.
C (c) Storing the quantities for two-parametric shooting:
IF(ISHEET.GT.0) THEN
AUX1=X1G1*X2G2-X2G1*X1G2
IF(AUX1.NE.0.) THEN
G1X1= X2G2/AUX1
G2X1=-X2G1/AUX1
G1X2=-X1G2/AUX1
G2X2= X1G1/AUX1
ELSE
C Caustic:
AUX1=X1G1*X1G1+X2G1*X2G1+X1G2*X1G2+X2G2*X2G2
IF(AUX1.NE.0.) THEN
C Line caustic - generalized inverse matrix:
G1X1=X1G1/AUX1
G2X1=X1G2/AUX1
G1X2=X2G1/AUX1
G2X2=X2G2/AUX1
ELSE
C Point caustic:
G1X1=0.
G2X1=0.
G1X2=0.
G2X2=0.
END IF
END IF
ELSE
X1=0.
X2=0.
G1X1=0.
G2X1=0.
G1X2=0.
G2X2=0.
END IF
C ****************************************************************
C Determination of two-point rays (to be moved to RPMEM):
IREC=0
IF(JTYPE.LT.-1000) THEN
C Two-point ray is being searched for
IF(ISHEET.GT.0) THEN
IREC=-JTYPE-1000
IF(IREC.LT.1.OR.NREC.LT.IREC) THEN
PAUSE 'Error *** in RPAR4: Wrong index of the receiver'
STOP
END IF
DIST2=(X1-XREC(1,IREC))**2+(X2-XREC(2,IREC))**2
IF(DIST2.GT.XERR**2) THEN
JTYPE=-2
IREC=0
END IF
ELSE
JTYPE=-2
END IF
END IF
C ****************************************************************
CALL RPMEM(JRAY,JTYPE,ISHEET,G1,G2,SG11,SG12,SG22,
* X1,X2,G1X1,G2X1,G1X2,G2X2)
C JRAY,JTYPE,G1,G2... Determined in RPAR2 (output of RP3D),
C X1,X2,X1G1,X2G1,X1G2,X2G2... Determined in RPAR32 (if defined),
C ISHEET,SG11,SG12,SG22,G1X1,G2X1,G1X2,G2X2... Determined in this
C subroutine.
END IF
C
C Determination of the area YI(22) and inverse specific moment
C YI(23:25) of element of the ray-parameter surface, corresponding
C to the ray:
IF(ISRFX(2).EQ.0) THEN
CALL RP2DG(YI(22),YI(23),YI(24),YI(25))
ELSE
YI(22)=0.
YI(23)=0.
YI(24)=0.
YI(25)=0.
C CALL RP3DG(YI(22),YI(23),YI(24),YI(25)) *** for future extension
END IF
C
RETURN
END
C
C=======================================================================
C
SUBROUTINE XFUN(IFUN,COOR,X,XDER)
INTEGER IFUN
REAL COOR(3),X,XDER(3)
C
C Interface subroutine designed to return the value and first
C derivatives of the X-functions describing the profiles. It is called
C by the subroutine RPAR32, and may be user-defined.
C
C Input:
C IFUN... 1 for the first, and 2 for the second receiver parameter
C (X-function).
C COOR... Array containing coordinates X1, X2, X3 of the given
C point.
C None of the input parameters are altered.
C
C Output:
C X... Value of the corresponding X-function at the given point.
C XDER... Array containing the derivatives of the X-function with
C respect to the general coordinates X1, X2, X3.
C
C Common block /RPARD/:
INCLUDE 'rpard.inc'
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
EXTERNAL SRFC2
C SRFC2... File 'srfc.for' of the package 'MODEL'.
C
C Date: 1993, December 20
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Auxiliary storage locations:
REAL FAUX(10)
C
IF(ISRFX(IFUN).GT.0) THEN
CALL SRFC2(ISRFX(IFUN),COOR,FAUX)
X=FAUX(1)
XDER(1)=FAUX(2)
XDER(2)=FAUX(3)
XDER(3)=FAUX(4)
ELSE
X=0.
XDER(1)=0.
XDER(2)=0.
XDER(3)=0.
IF(ISRFX(IFUN).LT.0) THEN
X=COOR(-ISRFX(IFUN))
XDER(-ISRFX(IFUN))=1.
END IF
END IF
RETURN
END
C
C=======================================================================
C
INTEGER FUNCTION IHIST(IEND)
INTEGER IEND
C
C Auxiliary integer function to the subroutine RPAR4 designed to
C determine the ray-history index. The ray-history index is used to
C distinguish between the rays with different histories for the purposes
C of the subroutine RPAR2.
C
C The ray history is recorded as a sequence of integers. In this
C version, the history consists of:
C (A) Indices of complex blocks along the ray (including the sign),
C (B) Indices of structural interfaces intersected by the ray (including
C the sign), and index of the end surface limiting the computational
C volume, if intersected.
C (C) Reason IEND (see subroutine RAY2 of 'ray.for' and INIT2 of
C 'init.for') if the ray is not successful (does not reach the
C reference surface).
C Note - differences in versions 4.10 and older:
C (B) Index of the end surface not considered in the history.
C (C) IEND not considered if IEND.GE.40.
C
C Input:
C IEND... Reason of the termination of the computation of the ray
C (see C.R.T.5.4). For a detailed description see
C subroutine RAY (subroutine file 'ray.for').
C None of the input parameters are altered.
C
C Output:
C IHIST...Positive (1,2,3,...) for successful ray,
C otherwise (0 or -1,-2,-3,...) for unsuccessful ray.
C Succesfull ray is understood here to be the ray having at
C least IABS(IPOINT) points of intersection with the
C specified surface ISRFR, see the input data (2).
C Different values of IHIST correspond to rays with
C different histories.
C
C Common block /RPARD/:
INCLUDE 'rpard.inc'
C None of the storage locations of the common block are altered.
C
C Common block /RPARH/:
INCLUDE 'rparc.inc'
C Storage locations of the common block may be altered.
C
C Subroutines and external functions required:
EXTERNAL SRPARH
C SRPARH... This file.
C
C Date: 1995, September 29
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Auxiliary storage locations:
INTEGER I,K
C
IF(IPOINT.EQ.0) THEN
IHIST=0
ELSE
IF(JPOINT.LT.IABS(IPOINT)) THEN
C Unsuccessful ray:
C KMAH index is not included in the history of an unsuccessful
C ray. It saves CPU time but may cause missing two-point rays,
C especially close to the boundary of the computational volume.
NRPARH=NRPARH+1
CALL SRPARH
IRPARH(NRPARH)=IEND
ELSE
C Successful ray:
NRPARH=NRPARH+1
CALL SRPARH
IRPARH(NRPARH)=KMAH
END IF
DO 13 K=1,IRPARH(0)
IF(IRPARH(K)-IRPARH(K-1).EQ.NRPARH-IRPARH(IRPARH(0))) THEN
DO 11 I=IRPARH(K-1)+1,IRPARH(K)
IF(IRPARH(I).NE.IRPARH(I-IRPARH(K)+NRPARH)) THEN
GO TO 12
END IF
11 CONTINUE
C the history of the last ray is equal to the K-th ray history
IHIST=INDMIN+K
GO TO 20
END IF
12 CONTINUE
13 CONTINUE
C New ray history:
NRPARH=NRPARH+1
CALL SRPARH
DO 18 I=NRPARH,IRPARH(0)+2,-1
IRPARH(I)=IRPARH(I-1)
18 CONTINUE
DO 19 I=1,IRPARH(0)
IRPARH(I)=IRPARH(I)+1
19 CONTINUE
IRPARH(0)=IRPARH(0)+1
IRPARH(IRPARH(0))=NRPARH
IHIST=INDMIN+IRPARH(0)
C
20 CONTINUE
IF(JPOINT.LT.IABS(IPOINT)) THEN
IHIST=-IHIST
END IF
IF(IPOINT.LT.0) THEN
IF(IEND.NE.10) THEN
IHIST=-IABS(IHIST)
END IF
END IF
END IF
RETURN
END
C
C=======================================================================
C
SUBROUTINE SRPARH()
C
C Auxiliary subroutine to the subroutine RPAR31 designed to shift the
C storage locations of the common block /RPARH/.
C
C No input.
C
C No output.
C
C Common block /RPARH/:
INCLUDE 'rparc.inc'
C Storage locations of the common block may be altered.
C
C No subroutines and external functions required.
C
C Error messages:
C 655... Too small array IRPARH in /RPARH/:
C Array IRPARH of the common block /RPARH/ designed to store
C different ray histories is too small to contain the last
C three different ray histories. The dimension MRPARH of
C the array IRPARH should be adjusted in all declarations of
C /RPARH/.
C
C Date: 1991, May 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Auxiliary storage locations:
INTEGER ISHIFT,I
C
C ISHIFT..Shift in the array IRPARH.
C I... Auxiliary loop variable.
C
C.......................................................................
C
IF(NRPARH.GT.MRPARH) THEN
IF(IRPARH(0).LE.2) THEN
PAUSE 'Error 655 in SRPARH: Too small array IRPARH in /RPARH/'
STOP
END IF
INDMIN=INDMIN+1
IRPARH(0)=IRPARH(0)-1
ISHIFT=IRPARH(1)-IRPARH(0)
NRPARH=NRPARH-ISHIFT
DO 1 I=1,IRPARH(0)
IRPARH(I)=IRPARH(I+1)-ISHIFT
1 CONTINUE
DO 2 I=IRPARH(0)+1,NRPARH-1
IRPARH(I)=IRPARH(I+ISHIFT)
2 CONTINUE
END IF
RETURN
END
C
C=======================================================================
C