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

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