[u/mrichter/AliRoot.git] / GEANT321 / ggeom / gnopg1.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:20:52  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.29  by  S.Giani
*-- Author :
      SUBROUTINE GNOPG1(X,P,SNXT)
*********************************************************************
***** GNOPG1 ********************************************************
*
*     GNOPG1 ... 20-DEC-1988
*     Version 1.1
*     Rolf Nierhaus
*
*********************************************************************
*
*     Copyright   CERN,   Geneva  1988  -  Copyright  and  any  other
*     appropriate legal protection of  these  computer  programs  and
*     associated  documentation  reserved  in  all  countries  of the
*     world.
*
*********************************************************************
*
*          Subroutine  GNOPG1 is called by GNOPGO for the computation
*     of SNXT, the distance from a point P  along  a  track  T  to  a
*     boundary  surface  of a Geant volume V of shape PGON. The point
*     P is outside the volume V. If the track T has  no  intersection
*     with the volume V, the vector distance SNXT is set to 1.E10.
*
*          V  is  generally  a composite volume consisting of several
*     sections. The sections have  boundary  surfaces  orthogonal  to
*     the   Z-axis.   Each  section  consists  generally  of  several
*     sectors. Each sector is an  "elementary"  convex  volume.  This
*     package  assumes it is either a hexahedron or a pentahedron. If
*     it is a pentahedron, it has 6 vertices, of  which  two  are  on
*     the  Z-axis.  All  sectors  of  the same section are congruent.
*     Each section has the same number of sectors.
*
*          GNOPG1 calls GNOPG2 for each section, GNOPG2  call  GNOPG3
*     for  each  sector.  GNOPG4  is  called  to  store  the  surface
*     parameters of  a  sector  in  the  common  block  GCQ1.  GNOPG6
*     computes  the  vector  distance  to  a  convex  volume.  GNOPG7
*     computes the  vector  distance  to  a  plane  surface.  Logical
*     function  GNOPG8  determines  if  a  point  is  inside a convex
*     volume, and logical function GNOPG9 determines if  a  point  is
*     inside a region delimited by a plane surface.
*
*          We  describe each surface by 6 parameters: the first three
*     are   the   coordinates   of   a   point   on    the    surface
*     XS(I),YS(I),ZS(I),  the  other  three are the components of the
*     normal vector of the surface XN(I),YN(I),ZN(I). I is the  index
*     of  the surface. We consider only one sector at a time, and the
*     number of boundary  surfaces  is  never  larger  then  6.  Each
*     surface  divides  the  space  into  two  regions:  the positive
*     region and the negative region. We choose the direction of  the
*     normal  vectors  of the boundary surfaces such that the bounded
*     volume is within the positive region of each surface, that  is,
*     the normal vector is pointing to the inside of the volume.
*
*          Logical   function   GNOPG9  returns  TRUE  if  the  point
*     (XP,YP,ZP) is within the positive region of  the  surface  with
*     index   I.   This   is  the  case  if  the  scalar  product  of
*     (XP-XS,YP-YS,ZP-ZS) and (XN,YN,ZN) is positive (or zero).
*
*          GNOPG8 is true if  the  point  (XP,YP,ZP)  is  within  the
*     positive region of all bounding surfaces.
*
*          To  find  the  distance  from  a  point (XP,YP,ZP) along a
*     track  with  direction  cosines   (XD,YD,ZD)   to   a   surface
*     (XS,YS,ZS)(XN,YN,ZN),  we  compute  first the scalar product of
*     the  vector  (XS-XP,YS-YP,ZS-ZP)   with   the   normal   vector
*     (XN,YN,ZN),  then  the scalar product of the vectors (XD,YD,ZD)
*     and (XN,YN,ZN).  The  first  scalar  product  is  the  shortest
*     distance  from  the  point  to  the  plane,  the  second scalar
*     product is the cosine of the angle between the  track  and  the
*     plane  normal.  The  quotient  is  the vector distance. If this
*     vector distance is  positive  (or  zero)  we  set  the  logical
*     variable  FLAG  TRUE.  GNOPG7  is  called with three parameters
*     I,FLAG and DIST. I is the index of the  surface,  and  DIST  is
*     the vector distance if FLAG is TRUE.
*
*          To  find the vector distance from a point to an elementary
*     volume, all bounding surfaces of the volume are considered.  If
*     the  point  is  in  the positive region of a surface, the track
*     could possibly exit through the surface, but  it  cannot  enter
*     through  it.  For  those surfaces which have the point in their
*     negative region, we  determine  if  the  track  intersects  the
*     surface,  and  what  is the distance to the intersection point.
*     Only the largest of these distances  is  a  candidate  for  the
*     distance  from  the point to the volume. It is however possible
*     that the track misses the elementary volume entirely.  This  we
*     find  out  by applying the function GNOPG8 (Inside volume test)
*     to the coordinates of the intersection point. GNOPG6 is  called
*     with  two parameters: a logical variable FLAG, which is TRUE if
*     the track  intersects  the  volume,  and  DIST,  which  is  the
*     distance if FLAG is TRUE.
*
*          The  coordinates of  the point P and the direction cosines
*     of the  track  T  are  stored  in  the  common  block  GCQ2  by
*     subroutine  GNOPG1. The parameter array P in the call to GNOPG1
*     contains in its first two members the lower  phi-limit  of  the
*     PGON  and  the range in phi. Both angles are in degrees. GNOPG1
*     convertes from degrees to radians and stores the phi-limits  of
*     the  first  sector  of  each  section  in the common block GCQ3
*     together with the number of sectors per section. The number  of
*     sectors  per  section  is  contained in the third member of the
*     parameter array P. The other members of P  have  the  following
*     meaning:  From  P(5)  onwards,  there are groups of three. Each
*     group describes a section boundary. The  first  member  is  the
*     Z-coordinate  of the boundary. The second and the third are the
*     distances from the Z-axis to the inner and outer edges at  that
*     boundary.  If  there is no inner edge, the sector is limited by
*     the Z-axis. In this case the second members of  the  group  are
*     zero.  Groups may be shared by adjacent sections. Otherwise the
*     Z-coordinates of two consecutive  groups  are  the  same.  P(4)
*     contains the number of groups.
*
*          GNOPG1  calls  GNOPG2  with  8 parameters. The first 6 are
*     input parameters, the first  3  refer  to  the  "left"  section
*     boundary,  the next 3 to the "right" section boundary. The last
*     two parameters are output. Logical variable  FLAG  is  TRUE  if
*     the  track  intersects  the  section.  In this case DIST is the
*     distance.
*
*          GNOPG2 calls GNOPG3 with four parameters. The first 2  are
*     input  parameters,  namely  the  phi-limits  of the sector. The
*     last two parameters again are  output:  FLAG  is  TRUE  if  the
*     track   intersects  the  sector.  In  this  case  DIST  is  the
*     distance. The other variables required by GNOPG3 are  the  same
*     for  all  sectors  of the same section and are stored by GNOPG2
*     in the common block GCQ4. GNOPG2  also  stores  in  the  common
*     block GCQ5 the number of boundary surfaces of each sector.
*
*          The  surfaces  orthogonal  to  the Z-axis are the same for
*     all sectors of a section. The surface parameters of  these  two
*     sections  are  stored  by  calls  to  GNOPG4  from  GNOPG2. The
*     surface parameters of the other  three  or  four  surfaces  are
*     stored from GNOPG3.
*
*          GNOPG3  sets  FLAG  TRUE,  if  the  track T intersects the
*     corresponding sector. GNOPG2 finds  the  shortest  distance  to
*     all  intersected  sectors,  and  set  FLAG TRUE, if the track T
*     intersects  the  corresponding  section.   GNOPG1   finds   the
*     shortest  distance  to  all intersected sections. If no section
*     intersects, the track FLAG is set FALSE, and 1.E10 is  returned
*     as SNXT.
*
***** Subroutine GNOPG1 *************************** 04-DEC-1988 *****
#if !defined(CERNLIB_SINGLE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (F=0.01745329251994330D0)
#endif
#if defined(CERNLIB_SINGLE)
      PARAMETER (F=0.01745329251994330)
#endif
      REAL X(6),P(49),SNXT
      LOGICAL FLAG,FLAG1
      DIMENSION XS(6), YS(6), ZS(6), XN(6), YN(6), ZN(6)
      LOGICAL FLAG1X, GNOP1X, GNOP2X, GNOP4X
      PARAMETER (ONE=1,HALF=ONE/2)
      XP=X(1)
      YP=X(2)
      ZP=X(3)
      XD=X(4)
      YD=X(5)
      ZD=X(6)
      NT=P(3)+.5
      P1=F*P(1)
      P2=F*P(2)/NT
      INDEX=5
      MINDEX=3.*P(4)+1.5
      FLAG=.FALSE.
      DIST=1.E10
   10 Z1=P(INDEX)
      D1N=P(INDEX+1)
      D1X=P(INDEX+2)
      Z2=P(INDEX+3)
      D2N=P(INDEX+4)
      D2X=P(INDEX+5)
C*****  Code Expanded From Routine:  GNOPG2
C*****  Code Expanded From Routine:  GNOPG4
      XS(1) = 0.
      YS(1) = 0.
      ZS(1) = Z1
      XN(1) = 0.
      YN(1) = 0.
      ZN(1) = 1.
C*****  End of Code Expanded From Routine:  GNOPG4
C*****  Code Expanded From Routine:  GNOPG4
      XS(2) = 0.
      YS(2) = 0.
      ZS(2) = Z2
      XN(2) = 0.
      YN(2) = 0.
      ZN(2) = -1.
C*****  End of Code Expanded From Routine:  GNOPG4
      P3 = P1
      P4 = P1 + P2
      Z3 = HALF*(Z1 + Z2)
      D3X = HALF*(D1X + D2X)
      TH1 = ATAN((D2X - D1X)/(Z2 - Z1))
      COSTH1 = COS(TH1)
      SINTH1 = SIN(TH1)
      D3N = HALF*(D1N + D2N)
      ISMAX = 5
      IF (D3N .NE. 0.) THEN
         ISMAX = 6
         TH2 = ATAN((D2N - D1N)/(Z2 - Z1))
         COSTH2 = COS(TH2)
         SINTH2 = SIN(TH2)
      ENDIF
      FLAG1 = .FALSE.
      DIST1 = 1.E10
      DO 60  I = 1, NT
C*****  Code Expanded From Routine:  GNOPG3
C*****  Code Expanded From Routine:  GNOPG4
         XS(3) = 0.
         YS(3) = 0.
         ZS(3) = Z3
         XN(3) = -SIN(P3)
         YN(3) = COS(P3)
         ZN(3) = 0.
C*****  End of Code Expanded From Routine:  GNOPG4
C*****  Code Expanded From Routine:  GNOPG4
         XS(4) = 0.
         YS(4) = 0.
         ZS(4) = Z3
         XN(4) = SIN(P4)
         YN(4) = -COS(P4)
         ZN(4) = 0.
C*****  End of Code Expanded From Routine:  GNOPG4
         P1X = HALF*(P3 + P4)
         COSP = COS(P1X)
         SINP = SIN(P1X)
C*****  Code Expanded From Routine:  GNOPG4
         XS(5) = D3X*COSP
         YS(5) = D3X*SINP
         ZS(5) = Z3
         XN(5) = -COSP*COSTH1
         YN(5) = -SINP*COSTH1
         ZN(5) = SINTH1
C*****  End of Code Expanded From Routine:  GNOPG4
         IF (D3N .NE. 0.) THEN
C*****  Code Expanded From Routine:  GNOPG4
            XS(6) = D3N*COSP
            YS(6) = D3N*SINP
            ZS(6) = Z3
            XN(6) = COSP*COSTH2
            YN(6) = SINP*COSTH2
            ZN(6) = -SINTH2
C*****  End of Code Expanded From Routine:  GNOPG4
         ENDIF
C*****  Code Expanded From Routine:  GNOPG6
         FLAG1X = .FALSE.
         DIST2X = 0.
         DO 20  IS = 1, ISMAX
C*****  Code Expanded From Routine:  GNOPG9
*     TRUE if (XP,YP,ZP) in positive region of surface I
            GNOP1X = 0. .LE. (XP-XS(IS))*XN(IS)+(YP-YS(IS))*YN(IS)+(ZP-
     +         ZS(IS))*ZN(IS)
C*****  End of Code Expanded From Routine:  GNOPG9
C*****  Code Expanded From Routine:  GNOPG7
            IF (.NOT.GNOP1X) THEN
               SPPMSN = (XP - XS(IS))*XN(IS) + (YP - YS(IS))*YN(IS) + (
     +            ZP - ZS(IS))*ZN(IS)
               SPDN = XD*XN(IS) + YD*YN(IS) + ZD*ZN(IS)
               IF (SPDN .EQ. 0.) THEN
                  DIST1X = -.000001
               ELSE
                  DIST1X = -SPPMSN/SPDN
               ENDIF
C*****  End of Code Expanded From Routine:  GNOPG7
               IF ((-1.E-5) .LE. DIST1X) THEN
                  FLAG1X = .TRUE.
                  DIST2X = MAX(DIST1X,DIST2X)
               ENDIF
            ENDIF
   20    CONTINUE
         IF (.NOT.FLAG1X) GO TO 50
         T = DIST2X + .001
         XQ = XP + T*XD
         YQ = YP + T*YD
         ZQ = ZP + T*ZD
C*****  Code Expanded From Routine:  GNOPG8
*     TRUE if (XP,YP,ZP) in volume
         GNOP2X = .FALSE.
         DO 30  IS1X = 1, ISMAX
C*****  Code Expanded From Routine:  GNOPG9
*     TRUE if (XP,YP,ZP) in positive region of surface I
            GNOP4X = 0. .LE. (XQ-XS(IS1X))*XN(IS1X)+(YQ-YS(IS1X))*YN(
     +         IS1X)+(ZQ-ZS(IS1X))*ZN(IS1X)
            IF (.NOT.GNOP4X) GO TO 40
C*****  End of Code Expanded From Routine:  GNOPG9
   30    CONTINUE
         GNOP2X = .TRUE.
   40    CONTINUE
         FLAG1X = GNOP2X
C*****  End of Code Expanded From Routine:  GNOPG8
   50    CONTINUE
C*****  End of Code Expanded From Routine:  GNOPG3
         IF (FLAG1X) THEN
            FLAG1 = .TRUE.
            DIST1 = MIN(DIST2X,DIST1)
         ENDIF
         P3 = P4
         P4 = P4 + P2
   60 CONTINUE
C*****  End of Code Expanded From Routine:  GNOPG2
      IF (FLAG1) THEN
         FLAG=.TRUE.
         IF (DIST1.LT.DIST) DIST=DIST1
      END IF
      INDEX=INDEX+3
      IF (MINDEX.LT.INDEX) THEN
         IF(FLAG) THEN
            SNXT=DIST
         ELSE
            SNXT=1.E10
         ENDIF
      ELSE
         IF (P(INDEX+3).EQ.Z2) INDEX=INDEX+3
         GO TO 10
      ENDIF
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1995/10/24 10:20:52 cernlib
	6	* Geant
	7	*
	8	*
	9	#include "geant321/pilot.h"
	10	*CMZ : 3.21/02 29/03/94 15.41.29 by S.Giani
	11	*-- Author :
	12	SUBROUTINE GNOPG1(X,P,SNXT)
	13	*********************************************************************
	14	*** GNOPG1 ******************************************************
	15	*
	16	* GNOPG1 ... 20-DEC-1988
	17	* Version 1.1
	18	* Rolf Nierhaus
	19	*
	20	*********************************************************************
	21	*
	22	* Copyright CERN, Geneva 1988 - Copyright and any other
	23	* appropriate legal protection of these computer programs and
	24	* associated documentation reserved in all countries of the
	25	* world.
	26	*
	27	*********************************************************************
	28	*
	29	* Subroutine GNOPG1 is called by GNOPGO for the computation
	30	* of SNXT, the distance from a point P along a track T to a
	31	* boundary surface of a Geant volume V of shape PGON. The point
	32	* P is outside the volume V. If the track T has no intersection
	33	* with the volume V, the vector distance SNXT is set to 1.E10.
	34	*
	35	* V is generally a composite volume consisting of several
	36	* sections. The sections have boundary surfaces orthogonal to
	37	* the Z-axis. Each section consists generally of several
	38	* sectors. Each sector is an "elementary" convex volume. This
	39	* package assumes it is either a hexahedron or a pentahedron. If
	40	* it is a pentahedron, it has 6 vertices, of which two are on
	41	* the Z-axis. All sectors of the same section are congruent.
	42	* Each section has the same number of sectors.
	43	*
	44	* GNOPG1 calls GNOPG2 for each section, GNOPG2 call GNOPG3
	45	* for each sector. GNOPG4 is called to store the surface
	46	* parameters of a sector in the common block GCQ1. GNOPG6
	47	* computes the vector distance to a convex volume. GNOPG7
	48	* computes the vector distance to a plane surface. Logical
	49	* function GNOPG8 determines if a point is inside a convex
	50	* volume, and logical function GNOPG9 determines if a point is
	51	* inside a region delimited by a plane surface.
	52	*
	53	* We describe each surface by 6 parameters: the first three
	54	* are the coordinates of a point on the surface
	55	* XS(I),YS(I),ZS(I), the other three are the components of the
	56	* normal vector of the surface XN(I),YN(I),ZN(I). I is the index
	57	* of the surface. We consider only one sector at a time, and the
	58	* number of boundary surfaces is never larger then 6. Each
	59	* surface divides the space into two regions: the positive
	60	* region and the negative region. We choose the direction of the
	61	* normal vectors of the boundary surfaces such that the bounded
	62	* volume is within the positive region of each surface, that is,
	63	* the normal vector is pointing to the inside of the volume.
	64	*
65	* Logical function GNOPG9 returns TRUE if the point
66	* (XP,YP,ZP) is within the positive region of the surface with
67	* index I. This is the case if the scalar product of
68	* (XP-XS,YP-YS,ZP-ZS) and (XN,YN,ZN) is positive (or zero).
69	*
70	* GNOPG8 is true if the point (XP,YP,ZP) is within the
71	* positive region of all bounding surfaces.
72	*
73	* To find the distance from a point (XP,YP,ZP) along a
74	* track with direction cosines (XD,YD,ZD) to a surface
75	* (XS,YS,ZS)(XN,YN,ZN), we compute first the scalar product of
76	* the vector (XS-XP,YS-YP,ZS-ZP) with the normal vector
77	* (XN,YN,ZN), then the scalar product of the vectors (XD,YD,ZD)
78	* and (XN,YN,ZN). The first scalar product is the shortest
79	* distance from the point to the plane, the second scalar
80	* product is the cosine of the angle between the track and the
81	* plane normal. The quotient is the vector distance. If this
82	* vector distance is positive (or zero) we set the logical
83	* variable FLAG TRUE. GNOPG7 is called with three parameters
84	* I,FLAG and DIST. I is the index of the surface, and DIST is
85	* the vector distance if FLAG is TRUE.
86	*
87	* To find the vector distance from a point to an elementary
88	* volume, all bounding surfaces of the volume are considered. If
89	* the point is in the positive region of a surface, the track
90	* could possibly exit through the surface, but it cannot enter
91	* through it. For those surfaces which have the point in their
92	* negative region, we determine if the track intersects the
93	* surface, and what is the distance to the intersection point.
94	* Only the largest of these distances is a candidate for the
95	* distance from the point to the volume. It is however possible
96	* that the track misses the elementary volume entirely. This we
97	* find out by applying the function GNOPG8 (Inside volume test)
98	* to the coordinates of the intersection point. GNOPG6 is called
99	* with two parameters: a logical variable FLAG, which is TRUE if
100	* the track intersects the volume, and DIST, which is the
101	* distance if FLAG is TRUE.
102	*
103	* The coordinates of the point P and the direction cosines
104	* of the track T are stored in the common block GCQ2 by
105	* subroutine GNOPG1. The parameter array P in the call to GNOPG1
106	* contains in its first two members the lower phi-limit of the
107	* PGON and the range in phi. Both angles are in degrees. GNOPG1
108	* convertes from degrees to radians and stores the phi-limits of
109	* the first sector of each section in the common block GCQ3
110	* together with the number of sectors per section. The number of
111	* sectors per section is contained in the third member of the
112	* parameter array P. The other members of P have the following
113	* meaning: From P(5) onwards, there are groups of three. Each
114	* group describes a section boundary. The first member is the
115	* Z-coordinate of the boundary. The second and the third are the
116	* distances from the Z-axis to the inner and outer edges at that
117	* boundary. If there is no inner edge, the sector is limited by
118	* the Z-axis. In this case the second members of the group are
119	* zero. Groups may be shared by adjacent sections. Otherwise the
120	* Z-coordinates of two consecutive groups are the same. P(4)
121	* contains the number of groups.
122	*
123	* GNOPG1 calls GNOPG2 with 8 parameters. The first 6 are
124	* input parameters, the first 3 refer to the "left" section
125	* boundary, the next 3 to the "right" section boundary. The last
126	* two parameters are output. Logical variable FLAG is TRUE if
127	* the track intersects the section. In this case DIST is the
128	* distance.
129	*
130	* GNOPG2 calls GNOPG3 with four parameters. The first 2 are
131	* input parameters, namely the phi-limits of the sector. The
132	* last two parameters again are output: FLAG is TRUE if the
133	* track intersects the sector. In this case DIST is the
134	* distance. The other variables required by GNOPG3 are the same
135	* for all sectors of the same section and are stored by GNOPG2
136	* in the common block GCQ4. GNOPG2 also stores in the common
137	* block GCQ5 the number of boundary surfaces of each sector.
138	*
139	* The surfaces orthogonal to the Z-axis are the same for
140	* all sectors of a section. The surface parameters of these two
141	* sections are stored by calls to GNOPG4 from GNOPG2. The
142	* surface parameters of the other three or four surfaces are
143	* stored from GNOPG3.
144	*
145	* GNOPG3 sets FLAG TRUE, if the track T intersects the
146	* corresponding sector. GNOPG2 finds the shortest distance to
147	* all intersected sectors, and set FLAG TRUE, if the track T
148	* intersects the corresponding section. GNOPG1 finds the
149	* shortest distance to all intersected sections. If no section
150	* intersects, the track FLAG is set FALSE, and 1.E10 is returned
151	* as SNXT.
152	*
153	*** Subroutine GNOPG1 *********************** 04-DEC-1988 ***
154	#if !defined(CERNLIB_SINGLE)
155	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
156	PARAMETER (F=0.01745329251994330D0)
157	#endif
158	#if defined(CERNLIB_SINGLE)
159	PARAMETER (F=0.01745329251994330)
160	#endif
161	REAL X(6),P(49),SNXT
162	LOGICAL FLAG,FLAG1
163	DIMENSION XS(6), YS(6), ZS(6), XN(6), YN(6), ZN(6)
164	LOGICAL FLAG1X, GNOP1X, GNOP2X, GNOP4X
165	PARAMETER (ONE=1,HALF=ONE/2)
166	XP=X(1)
167	YP=X(2)
168	ZP=X(3)
169	XD=X(4)
170	YD=X(5)
171	ZD=X(6)
172	NT=P(3)+.5
173	P1=F*P(1)
174	P2=F*P(2)/NT
175	INDEX=5
176	MINDEX=3.*P(4)+1.5
177	FLAG=.FALSE.
178	DIST=1.E10
179	10 Z1=P(INDEX)
180	D1N=P(INDEX+1)
181	D1X=P(INDEX+2)
182	Z2=P(INDEX+3)
183	D2N=P(INDEX+4)
184	D2X=P(INDEX+5)
185	C***** Code Expanded From Routine: GNOPG2
186	C***** Code Expanded From Routine: GNOPG4
187	XS(1) = 0.
188	YS(1) = 0.
189	ZS(1) = Z1
190	XN(1) = 0.
191	YN(1) = 0.
192	ZN(1) = 1.
193	C***** End of Code Expanded From Routine: GNOPG4
194	C***** Code Expanded From Routine: GNOPG4
195	XS(2) = 0.
196	YS(2) = 0.
197	ZS(2) = Z2
198	XN(2) = 0.
199	YN(2) = 0.
200	ZN(2) = -1.
201	C***** End of Code Expanded From Routine: GNOPG4
202	P3 = P1
203	P4 = P1 + P2
204	Z3 = HALF*(Z1 + Z2)
205	D3X = HALF*(D1X + D2X)
206	TH1 = ATAN((D2X - D1X)/(Z2 - Z1))
207	COSTH1 = COS(TH1)
208	SINTH1 = SIN(TH1)
209	D3N = HALF*(D1N + D2N)
210	ISMAX = 5
211	IF (D3N .NE. 0.) THEN
212	ISMAX = 6
213	TH2 = ATAN((D2N - D1N)/(Z2 - Z1))
214	COSTH2 = COS(TH2)
215	SINTH2 = SIN(TH2)
216	ENDIF
217	FLAG1 = .FALSE.
218	DIST1 = 1.E10
219	DO 60 I = 1, NT
220	C***** Code Expanded From Routine: GNOPG3
221	C***** Code Expanded From Routine: GNOPG4
222	XS(3) = 0.
223	YS(3) = 0.
224	ZS(3) = Z3
225	XN(3) = -SIN(P3)
226	YN(3) = COS(P3)
227	ZN(3) = 0.
228	C***** End of Code Expanded From Routine: GNOPG4
229	C***** Code Expanded From Routine: GNOPG4
230	XS(4) = 0.
231	YS(4) = 0.
232	ZS(4) = Z3
233	XN(4) = SIN(P4)
234	YN(4) = -COS(P4)
235	ZN(4) = 0.
236	C***** End of Code Expanded From Routine: GNOPG4
237	P1X = HALF*(P3 + P4)
238	COSP = COS(P1X)
239	SINP = SIN(P1X)
240	C***** Code Expanded From Routine: GNOPG4
241	XS(5) = D3X*COSP
242	YS(5) = D3X*SINP
243	ZS(5) = Z3
244	XN(5) = -COSP*COSTH1
245	YN(5) = -SINP*COSTH1
246	ZN(5) = SINTH1
247	C***** End of Code Expanded From Routine: GNOPG4
248	IF (D3N .NE. 0.) THEN
249	C***** Code Expanded From Routine: GNOPG4
250	XS(6) = D3N*COSP
251	YS(6) = D3N*SINP
252	ZS(6) = Z3
253	XN(6) = COSP*COSTH2
254	YN(6) = SINP*COSTH2
255	ZN(6) = -SINTH2
256	C***** End of Code Expanded From Routine: GNOPG4
257	ENDIF
258	C***** Code Expanded From Routine: GNOPG6
259	FLAG1X = .FALSE.
260	DIST2X = 0.
261	DO 20 IS = 1, ISMAX
262	C***** Code Expanded From Routine: GNOPG9
263	* TRUE if (XP,YP,ZP) in positive region of surface I
264	GNOP1X = 0. .LE. (XP-XS(IS))XN(IS)+(YP-YS(IS))YN(IS)+(ZP-
265	+ ZS(IS))*ZN(IS)
266	C***** End of Code Expanded From Routine: GNOPG9
267	C***** Code Expanded From Routine: GNOPG7
268	IF (.NOT.GNOP1X) THEN
269	SPPMSN = (XP - XS(IS))XN(IS) + (YP - YS(IS))YN(IS) + (
270	+ ZP - ZS(IS))*ZN(IS)
271	SPDN = XDXN(IS) + YDYN(IS) + ZD*ZN(IS)
272	IF (SPDN .EQ. 0.) THEN
273	DIST1X = -.000001
274	ELSE
275	DIST1X = -SPPMSN/SPDN
276	ENDIF
277	C***** End of Code Expanded From Routine: GNOPG7
278	IF ((-1.E-5) .LE. DIST1X) THEN
279	FLAG1X = .TRUE.
280	DIST2X = MAX(DIST1X,DIST2X)
281	ENDIF
282	ENDIF
283	20 CONTINUE
284	IF (.NOT.FLAG1X) GO TO 50
285	T = DIST2X + .001
286	XQ = XP + T*XD
287	YQ = YP + T*YD
288	ZQ = ZP + T*ZD
289	C***** Code Expanded From Routine: GNOPG8
290	* TRUE if (XP,YP,ZP) in volume
291	GNOP2X = .FALSE.
292	DO 30 IS1X = 1, ISMAX
293	C***** Code Expanded From Routine: GNOPG9
294	* TRUE if (XP,YP,ZP) in positive region of surface I
295	GNOP4X = 0. .LE. (XQ-XS(IS1X))XN(IS1X)+(YQ-YS(IS1X))YN(
296	+ IS1X)+(ZQ-ZS(IS1X))*ZN(IS1X)
297	IF (.NOT.GNOP4X) GO TO 40
298	C***** End of Code Expanded From Routine: GNOPG9
299	30 CONTINUE
300	GNOP2X = .TRUE.
301	40 CONTINUE
302	FLAG1X = GNOP2X
303	C***** End of Code Expanded From Routine: GNOPG8
304	50 CONTINUE
305	C***** End of Code Expanded From Routine: GNOPG3
306	IF (FLAG1X) THEN
307	FLAG1 = .TRUE.
308	DIST1 = MIN(DIST2X,DIST1)
309	ENDIF
310	P3 = P4
311	P4 = P4 + P2
312	60 CONTINUE
313	C***** End of Code Expanded From Routine: GNOPG2
314	IF (FLAG1) THEN
315	FLAG=.TRUE.
316	IF (DIST1.LT.DIST) DIST=DIST1
317	END IF
318	INDEX=INDEX+3
319	IF (MINDEX.LT.INDEX) THEN
320	IF(FLAG) THEN
321	SNXT=DIST
322	ELSE
323	SNXT=1.E10
324	ENDIF
325	ELSE
326	IF (P(INDEX+3).EQ.Z2) INDEX=INDEX+3
327	GO TO 10
328	ENDIF
329	END