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

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1999/05/18 15:55:17  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:20:54  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.30  by  S.Giani
*-- Author :
      SUBROUTINE GNPGON (X, PAR, IACT, SNEXT, SNXT, SAFE)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *       COMPUTE DISTANCE UP TO INTERSECTION WITH 'PGON' VOLUME,  *
C.    *       FROM INSIDE POINT X(1-3) ALONG DIRECTION X(4-6)          *
C.    *                                                                *
C.    *       PAR   (input)  : volume parameters                       *
C.    *       IACT  (input)  : action flag                             *
C.    *         = 0  Compute SAFE only                                 *
C.    *         = 1  Compute SAFE, and SNXT only if SNEXT .GT.new SAFE *
C.    *         = 2  Compute both SAFE and SNXT                        *
C.    *         = 3  Compute SNXT only                                 *
C.    *       SNEXT (input)  : see IACT = 1                            *
C.    *       SNXT  (output) : distance to volume boundary             *
C.    *       SAFE  (output) : shortest distance to any boundary       *
C.    *                                                                *
C.    *    ==>Called by : GNEXT, GTNEXT                                *
C.    *         Author  A.McPherson,  P.Weidhaas  *********            *
C.    *                                                                *
C.    ******************************************************************
C.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "geant321/gconsp.inc"
 
      REAL X(6), PAR(50), SNEXT, SNXT, SAFE
      DIMENSION TPAR(50)
 
C-----------------------------------------------------------------
 
      SNXT = BIG
      SAFSEG = BIG
      R2 = X(1)*X(1) + X(2)*X(2)
      R  = SQRT (R2)
 
      PDIV  = PAR(2) / PAR(3)
      DELPHI  = PDIV * DEGRAD
      DPHI2 = 0.5 * DELPHI
      CSDPH2 = COS (DPHI2)
 
      TPAR(1)  = PAR(1)
      TPAR(2)  = PAR(2)
      TPAR(3)  = PAR(3)
      TPAR(4)  = PAR(4)
      NZ  = PAR(4)
 
      DO 5 I=1, NZ
        I3  = 3*(I-1)
        TPAR(5+I3) = PAR(5+I3)
        TPAR(6+I3) = PAR(6+I3) / CSDPH2
        TPAR(7+I3) = PAR(7+I3)
    5 CONTINUE
 
C**************************************************************************
C
C......  Here we start the logic from "GNPCON" (which for reasons of
C......  efficiency and clarity has been implemented inline).
C
C**************************************************************************
 
      ZMIN  = TPAR(5)
      ZMAX  = TPAR(3*NZ+2)
      SAFZ1 = X(3) - ZMIN
      SAFZ2 = ZMAX - X(3)
      SAFEZ = MIN (SAFZ1, SAFZ2)
 
C
C......  First determine in which segment the particle is located.
C
      DO 10 JPH=8, 3*NZ-1,  3
        IF (X(3) .LT. TPAR(JPH)) THEN
           IPH=JPH
           GO TO 20
        ENDIF
   10 CONTINUE
      IPH = 3*NZ+2
 
   20 CONTINUE
C
C......  The particle is in the segment bounded by z-planes at
C......  Z1=PAR(IPL) and Z2=PAR(IPH), i.e.,  Z1 < X(3) < Z2.
C
C......  Set parameters for this segment and translate z-coordinate
C......  of point relative to center of this segment. this is done in
C......  preparation of invoking the algorithms used in "GNTUBE" and
C......  "GNCONE" (which for reasons of efficiency and clarity are
C......  implemented inline).
C
      IPL  = IPH - 3
      DZ    = 0.5 * (TPAR(IPH) - TPAR(IPL))
      PT2 = TPAR(IPL+1)
      PT3 = TPAR(IPL+2)
      PT4 = TPAR(IPH+1)
      PT5 = TPAR(IPH+2)
      PT6 = TPAR(1)
      PT7 = TPAR(1) + TPAR(2)
      IF (PT7 .GT. 360.0)  PT7  = PT7 - 360.0
 
      XT3 = X(3) - 0.5 * (TPAR(IPL) + TPAR(IPH))
 
      SAFZ2  = DZ + XT3
      ZLENI  = 0.5 / DZ
 
      PHI1  = PT6 * DEGRAD
      PHI2  = PT7 * DEGRAD
      IF (PHI2.LE.PHI1) PHI2=PHI2+TWOPI
 
      IF (IACT .LT. 3) THEN
 
C       -------------------------------------------------
C       |  Compute safety-distance 'SAFE' (P.Weidhaas)  |
C       -------------------------------------------------
 
        IFLAG = 2
        IF (TPAR(2) .EQ. 360.0)  IFLAG = 1
        SAFZ1 = DZ - ABS(XT3)
        SAFEZ = MIN (SAFEZ,SAFZ1)
C
C......  Next determine whether the segment is a tube or a cone.
C
        IF (PT2 .NE. PT4) GO TO 50
        IF (PT3 .NE. PT5) GO TO 50
 
C*********************************************************
C
C......  The segment is a tube; invoke the algorithm
C......  from routine "GNTUBE" inline to get "SAFER".
C
C*********************************************************
 
        SAFR1  = BIG
        IF(PT2.GT.0.) SAFR1  = R - PT2
        SAFR2  = PT3 - R
        SAFER  = MIN (SAFR1, SAFR2)
 
        IF (IFLAG .EQ. 2) GO TO 70
 
        GO TO 100
 
 
   50   CONTINUE
 
C*********************************************************
C
C......  The segment is a cone; invoke the algorithm
C......  from routine "GNCONE" inline to get "SAFER".
C
C*********************************************************
 
 
C......  Compute radial distance to inner wall.
 
        IF(PT2+PT4.GT.0.) THEN
        FACT  = (PT4 - PT2) * ZLENI
        RAD1  = PT2 + FACT * SAFZ2
        SAFR1 = (R - RAD1) / SQRT(1.0 + FACT*FACT)
        ELSE
        SAFR1 = BIG
        ENDIF
 
C......  Compute radial distance to outer wall.
 
        FACT  = (PT5 - PT3) * ZLENI
        RAD2  = PT3 + FACT * SAFZ2
        SAFR2 = (RAD2 - R) / SQRT(1.0 + FACT*FACT)
 
        SAFER  = MIN (SAFR1, SAFR2)
 
        IF (IFLAG .EQ. 1) GO TO 100
 
   70   CONTINUE
 
C********************************************************************
C......  Here we handle the case of a phi-segment of a tube or cone.
C......  in addition to the radial distances (SAFR1, SAFR2) and the
C......  axial distances (SAFZ1, SAFZ2) we compute here the distance
C......  to the phi-segment boundary that is closest to the point:
C
C......  For each phi-boundary we find the distance from the given
C......  point to the outer (at R2) point of the segment boundary
C......  (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
C......  "SAFSEG" to be the distance to segment PHI1, else we set
C......  "SAFSEG" to be the distance to segment PHI2.
C*********************************************************************
 
 
          COSPH1  = COS (PHI1)
          SINPH1  = SIN (PHI1)
          COSPH2  = COS (PHI2)
          SINPH2  = SIN (PHI2)
 
C......  Get coordinates of outer endpoints (at R2) of both phi-segments.
 
          XS1  = R * COSPH1
          YS1  = R * SINPH1
          XS2  = R * COSPH2
          YS2  = R * SINPH2
 
C......  Get distances (squared) from given point to each endpoint.
 
          DISTS1 = (X(1) - XS1)**2  +  (X(2) - YS1)**2
          DISTS2 = (X(1) - XS2)**2  +  (X(2) - YS2)**2
 
C......  Get distance to that phi-segment whose endpoint
C......  is closest to the given point.
 
          IF (DISTS1 .LE. DISTS2) THEN
            SAFSEG = ABS(SINPH1 * X(1) - COSPH1 * X(2))
          ELSE
            SAFSEG = ABS(SINPH2 * X(1) - COSPH2 * X(2))
          ENDIF
 
  100   CONTINUE
 
 
        IF (SAFER .LE. 0.0)  THEN
 
C---------------------------------------------------------------------------
C
C......  Here we handle the case in which  SAFER < 0, i.e., the point is
C......  inside the polygon but outside the inscribed polycone. We must
C......  do an accurate calculation of "SAFER".
C
C---------------------------------------------------------------------------
 
          FACT = SAFZ2 * ZLENI
          RAD1 = PT2 + FACT * (PT4 - PT2)
          RAD2 = PT3 + FACT * (PT5 - PT3)
          RR1 = RAD1 * RAD1
          RR2 = RAD2 * RAD2
 
          IF(X(1).EQ.0.)THEN
             PHI   = ATAN2(X(2), 1.E-8)
          ELSE
             PHI   = ATAN2(X(2), X(1))
          ENDIF
          IF (PHI .LT. PHI1)  PHI = PHI + TWOPI
 
          DIST=0.
          IF (PHI .GE. PHI1  .AND.  PHI .LE. PHI2)  THEN
            PHIREL = PHI - PHI1
            NSECTR = INT(PHIREL / DELPHI) + 1
            PHICTR  = PHI1 + (2.0*NSECTR - 1.0) * DPHI2
            COSPHC = COS (PHICTR)
            SINPHC = SIN (PHICTR)
 
 
            IF (R2 .GE. RR2) THEN
              SR2 = RAD2
              DIST = ABS (COSPHC * X(1) + SINPHC * X(2) - SR2)
              FACTC  = (PT5 - PT3) * ZLENI
            ELSEIF (R2 .LE. RR1) THEN
              RIN  = RAD1 * CSDPH2
              SRIN = RIN
              DIST = ABS (COSPHC * X(1) + SINPHC * X(2) - SRIN)
              FACTC = (PT4 - PT2) * ZLENI
            ELSE
              PRINT *, 'GNPGON: FACTC not calculated!!!'
              FACTC=0
            ENDIF
 
          ENDIF
 
          SAFER  = DIST / SQRT(1.0 + FACTC*FACTC)
 
        ENDIF
 
        SAFE  = MIN (SAFEZ, SAFER, SAFSEG)
 
 
        IF (IACT .EQ. 0) GO TO 999
        IF (IACT .EQ. 1) THEN
          IF (SNEXT .LT. SAFE) GO TO 999
        ENDIF
      ENDIF
 
 
C     ------------------------------------------------
C     |  Compute vector-distance 'SNXT' (Nierhaus)  |
C     ------------------------------------------------
 
 
      CALL GNPGO1(X,PAR,SNXT)
C
  999 CONTINUE
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
d80f58a7	5	* Revision 1.1.1.1 1999/05/18 15:55:17 fca
	6	* AliRoot sources
	7	*
fe4da5cc	8	* Revision 1.1.1.1 1995/10/24 10:20:54 cernlib
	9	* Geant
	10	*
	11	*
	12	#include "geant321/pilot.h"
	13	*CMZ : 3.21/02 29/03/94 15.41.30 by S.Giani
	14	*-- Author :
	15	SUBROUTINE GNPGON (X, PAR, IACT, SNEXT, SNXT, SAFE)
	16	C.
	17	C. ******************************************************************
	18	C. * *
	19	C. * COMPUTE DISTANCE UP TO INTERSECTION WITH 'PGON' VOLUME, *
	20	C. * FROM INSIDE POINT X(1-3) ALONG DIRECTION X(4-6) *
	21	C. * *
	22	C. * PAR (input) : volume parameters *
	23	C. * IACT (input) : action flag *
	24	C. * = 0 Compute SAFE only *
	25	C. * = 1 Compute SAFE, and SNXT only if SNEXT .GT.new SAFE *
	26	C. * = 2 Compute both SAFE and SNXT *
	27	C. * = 3 Compute SNXT only *
	28	C. * SNEXT (input) : see IACT = 1 *
	29	C. * SNXT (output) : distance to volume boundary *
	30	C. * SAFE (output) : shortest distance to any boundary *
	31	C. * *
	32	C. * ==>Called by : GNEXT, GTNEXT *
	33	C. * Author A.McPherson, P.Weidhaas ********* *
	34	C. * *
	35	C. ******************************************************************
	36	C.
d80f58a7	37	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
fe4da5cc	38	#include "geant321/gconsp.inc"
fe4da5cc	39
d80f58a7	40	REAL X(6), PAR(50), SNEXT, SNXT, SAFE
d80f58a7	41	DIMENSION TPAR(50)
fe4da5cc	42
	43	C-----------------------------------------------------------------
	44
	45	SNXT = BIG
	46	SAFSEG = BIG
	47	R2 = X(1)X(1) + X(2)X(2)
	48	R = SQRT (R2)
	49
	50	PDIV = PAR(2) / PAR(3)
	51	DELPHI = PDIV * DEGRAD
	52	DPHI2 = 0.5 * DELPHI
	53	CSDPH2 = COS (DPHI2)
	54
	55	TPAR(1) = PAR(1)
	56	TPAR(2) = PAR(2)
	57	TPAR(3) = PAR(3)
	58	TPAR(4) = PAR(4)
	59	NZ = PAR(4)
	60
	61	DO 5 I=1, NZ
	62	I3 = 3*(I-1)
	63	TPAR(5+I3) = PAR(5+I3)
	64	TPAR(6+I3) = PAR(6+I3) / CSDPH2
	65	TPAR(7+I3) = PAR(7+I3)
	66	5 CONTINUE
	67
	68	C**************************************************************************
	69	C
	70	C...... Here we start the logic from "GNPCON" (which for reasons of
	71	C...... efficiency and clarity has been implemented inline).
	72	C
	73	C**************************************************************************
	74
	75	ZMIN = TPAR(5)
	76	ZMAX = TPAR(3*NZ+2)
	77	SAFZ1 = X(3) - ZMIN
	78	SAFZ2 = ZMAX - X(3)
	79	SAFEZ = MIN (SAFZ1, SAFZ2)
	80
	81	C
	82	C...... First determine in which segment the particle is located.
	83	C
	84	DO 10 JPH=8, 3*NZ-1, 3
	85	IF (X(3) .LT. TPAR(JPH)) THEN
	86	IPH=JPH
	87	GO TO 20
	88	ENDIF
	89	10 CONTINUE
	90	IPH = 3*NZ+2
	91
	92	20 CONTINUE
	93	C
	94	C...... The particle is in the segment bounded by z-planes at
	95	C...... Z1=PAR(IPL) and Z2=PAR(IPH), i.e., Z1 < X(3) < Z2.
	96	C
	97	C...... Set parameters for this segment and translate z-coordinate
	98	C...... of point relative to center of this segment. this is done in
	99	C...... preparation of invoking the algorithms used in "GNTUBE" and
	100	C...... "GNCONE" (which for reasons of efficiency and clarity are
	101	C...... implemented inline).
	102	C
	103	IPL = IPH - 3
	104	DZ = 0.5 * (TPAR(IPH) - TPAR(IPL))
	105	PT2 = TPAR(IPL+1)
106	PT3 = TPAR(IPL+2)
107	PT4 = TPAR(IPH+1)
108	PT5 = TPAR(IPH+2)
109	PT6 = TPAR(1)
110	PT7 = TPAR(1) + TPAR(2)
111	IF (PT7 .GT. 360.0) PT7 = PT7 - 360.0
112
113	XT3 = X(3) - 0.5 * (TPAR(IPL) + TPAR(IPH))
114
115	SAFZ2 = DZ + XT3
116	ZLENI = 0.5 / DZ
117
118	PHI1 = PT6 * DEGRAD
119	PHI2 = PT7 * DEGRAD
120	IF (PHI2.LE.PHI1) PHI2=PHI2+TWOPI
121
122	IF (IACT .LT. 3) THEN
123
124	C -------------------------------------------------
125	C \| Compute safety-distance 'SAFE' (P.Weidhaas) \|
126	C -------------------------------------------------
127
128	IFLAG = 2
129	IF (TPAR(2) .EQ. 360.0) IFLAG = 1
130	SAFZ1 = DZ - ABS(XT3)
131	SAFEZ = MIN (SAFEZ,SAFZ1)
132	C
133	C...... Next determine whether the segment is a tube or a cone.
134	C
135	IF (PT2 .NE. PT4) GO TO 50
136	IF (PT3 .NE. PT5) GO TO 50
137
138	C*********************************************************
139	C
140	C...... The segment is a tube; invoke the algorithm
141	C...... from routine "GNTUBE" inline to get "SAFER".
142	C
143	C*********************************************************
144
145	SAFR1 = BIG
146	IF(PT2.GT.0.) SAFR1 = R - PT2
147	SAFR2 = PT3 - R
148	SAFER = MIN (SAFR1, SAFR2)
149
150	IF (IFLAG .EQ. 2) GO TO 70
151
152	GO TO 100
153
154
155	50 CONTINUE
156
157	C*********************************************************
158	C
159	C...... The segment is a cone; invoke the algorithm
160	C...... from routine "GNCONE" inline to get "SAFER".
161	C
162	C*********************************************************
163
164
165	C...... Compute radial distance to inner wall.
166
167	IF(PT2+PT4.GT.0.) THEN
168	FACT = (PT4 - PT2) * ZLENI
169	RAD1 = PT2 + FACT * SAFZ2
170	SAFR1 = (R - RAD1) / SQRT(1.0 + FACT*FACT)
171	ELSE
172	SAFR1 = BIG
173	ENDIF
174
175	C...... Compute radial distance to outer wall.
176
177	FACT = (PT5 - PT3) * ZLENI
178	RAD2 = PT3 + FACT * SAFZ2
179	SAFR2 = (RAD2 - R) / SQRT(1.0 + FACT*FACT)
180
181	SAFER = MIN (SAFR1, SAFR2)
182
183	IF (IFLAG .EQ. 1) GO TO 100
184
185	70 CONTINUE
186
187	C********************************************************************
188	C...... Here we handle the case of a phi-segment of a tube or cone.
189	C...... in addition to the radial distances (SAFR1, SAFR2) and the
190	C...... axial distances (SAFZ1, SAFZ2) we compute here the distance
191	C...... to the phi-segment boundary that is closest to the point:
192	C
193	C...... For each phi-boundary we find the distance from the given
194	C...... point to the outer (at R2) point of the segment boundary
195	C...... (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
196	C...... "SAFSEG" to be the distance to segment PHI1, else we set
197	C...... "SAFSEG" to be the distance to segment PHI2.
198	C*********************************************************************
199
200
201	COSPH1 = COS (PHI1)
202	SINPH1 = SIN (PHI1)
203	COSPH2 = COS (PHI2)
204	SINPH2 = SIN (PHI2)
205
206	C...... Get coordinates of outer endpoints (at R2) of both phi-segments.
207
208	XS1 = R * COSPH1
209	YS1 = R * SINPH1
210	XS2 = R * COSPH2
211	YS2 = R * SINPH2
212
213	C...... Get distances (squared) from given point to each endpoint.
214
215	DISTS1 = (X(1) - XS1)2 + (X(2) - YS1)2
216	DISTS2 = (X(1) - XS2)2 + (X(2) - YS2)2
217
218	C...... Get distance to that phi-segment whose endpoint
219	C...... is closest to the given point.
220
221	IF (DISTS1 .LE. DISTS2) THEN
222	SAFSEG = ABS(SINPH1 * X(1) - COSPH1 * X(2))
223	ELSE
224	SAFSEG = ABS(SINPH2 * X(1) - COSPH2 * X(2))
225	ENDIF
226
227	100 CONTINUE
228
229
230	IF (SAFER .LE. 0.0) THEN
231
232	C---------------------------------------------------------------------------
233	C
234	C...... Here we handle the case in which SAFER < 0, i.e., the point is
235	C...... inside the polygon but outside the inscribed polycone. We must
236	C...... do an accurate calculation of "SAFER".
237	C
238	C---------------------------------------------------------------------------
239
240	FACT = SAFZ2 * ZLENI
241	RAD1 = PT2 + FACT * (PT4 - PT2)
242	RAD2 = PT3 + FACT * (PT5 - PT3)
243	RR1 = RAD1 * RAD1
244	RR2 = RAD2 * RAD2
245
246	IF(X(1).EQ.0.)THEN
247	PHI = ATAN2(X(2), 1.E-8)
248	ELSE
249	PHI = ATAN2(X(2), X(1))
250	ENDIF
251	IF (PHI .LT. PHI1) PHI = PHI + TWOPI
252
253	DIST=0.
254	IF (PHI .GE. PHI1 .AND. PHI .LE. PHI2) THEN
255	PHIREL = PHI - PHI1
256	NSECTR = INT(PHIREL / DELPHI) + 1
257	PHICTR = PHI1 + (2.0NSECTR - 1.0) DPHI2
258	COSPHC = COS (PHICTR)
259	SINPHC = SIN (PHICTR)
260
261
262	IF (R2 .GE. RR2) THEN
263	SR2 = RAD2
264	DIST = ABS (COSPHC * X(1) + SINPHC * X(2) - SR2)
265	FACTC = (PT5 - PT3) * ZLENI
266	ELSEIF (R2 .LE. RR1) THEN
267	RIN = RAD1 * CSDPH2
268	SRIN = RIN
269	DIST = ABS (COSPHC * X(1) + SINPHC * X(2) - SRIN)
270	FACTC = (PT4 - PT2) * ZLENI
d80f58a7	271	ELSE
	272	PRINT *, 'GNPGON: FACTC not calculated!!!'
	273	FACTC=0
fe4da5cc	274	ENDIF
	275
	276	ENDIF
	277
	278	SAFER = DIST / SQRT(1.0 + FACTC*FACTC)
	279
	280	ENDIF
	281
	282	SAFE = MIN (SAFEZ, SAFER, SAFSEG)
	283
	284
	285	IF (IACT .EQ. 0) GO TO 999
	286	IF (IACT .EQ. 1) THEN
	287	IF (SNEXT .LT. SAFE) GO TO 999
	288	ENDIF
	289	ENDIF
	290
	291
	292	C ------------------------------------------------
	293	C \| Compute vector-distance 'SNXT' (Nierhaus) \|
	294	C ------------------------------------------------
	295
	296
	297	CALL GNPGO1(X,PAR,SNXT)
	298	C
	299	999 CONTINUE
	300	END
	301