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

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