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

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:20:53  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.30  by  S.Giani
*-- Author :
      SUBROUTINE GNPCON (X, PAR, IACT, SNEXT, SNXT, SAFE)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *       COMPUTE DISTANCE UP TO INTERSECTION WITH 'PCON' 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(9),  XT(6), PT(7)
 
      EQUIVALENCE (PT(1), DZ),  (PT(2), PT2), (PT(3), PT3)
      EQUIVALENCE (PT(4), PT4), (PT(5), PT5), (PT(6), PT6)
      EQUIVALENCE (PT(7), PT7)
 
      EQUIVALENCE (XT(3), XT3)
C.
C.   ---------------------------------------------------
C.
 
      SNXT = BIG
      R2  = X(1)*X(1) + X(2)*X(2)
      R   = SQRT (R2)
 
      NZ  = PAR(3)
      ZMIN  = PAR(4)
      ZMAX  = PAR(3*NZ+1)
      SAFZ1 = X(3) - ZMIN
      SAFZ2 = ZMAX - X(3)
      SAFEZ = MIN (SAFZ1, SAFZ2)
 
C
C......  First determine in which z-segment the particle is located.
C
      DO 10 JPH=7, 3*(NZ-1)+1,  3
        IF (X(3) .LE. PAR(JPH)) THEN
           IPH=JPH
           GO TO 20
        ENDIF
   10 CONTINUE
      IPH=3*NZ+1
 
   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 * (PAR(IPH) - PAR(IPL))
      PT2 = PAR(IPL+1)
      PT3 = PAR(IPL+2)
      PT4 = PAR(IPH+1)
      PT5 = PAR(IPH+2)
      PT6 = PAR(1)
      PT7 = PAR(1) + PAR(2)
      IF (PT7 .GT. 360.0)  PT7  = PT7 - 360.0
 
      XT3 = X(3) - 0.5 * (PAR(IPL) + PAR(IPH))
 
      XT(1)  = X(1)
      XT(2)  = X(2)
      XT(4)  = X(4)
      XT(5)  = X(5)
      XT(6)  = X(6)
 
      IND = 2
      IF (PAR(2) .EQ. 360.0)  IND = 1
 
      IF (IACT .LT. 3) THEN
 
C       -------------------------------------------------
C       |  Compute safety-distance 'SAFE' (P.Weidhaas)  |
C       -------------------------------------------------
        SAFZ1   = DZ - ABS(XT3)
        SAFEZ   = MIN (SAFEZ,SAFZ1)
        SAFSEG  = BIG
C
C......  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  = R - PT2
        SAFR2  = PT3 - R
        SAFER  = MIN (SAFR1, SAFR2)
 
        IF (IND .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*********************************************************
 
        SAFZ2  = DZ + XT3
        ZLENI  = 0.5 / DZ
 
C......  Compute radial distance to inner wall.
 
        FACT  = (PT4 - PT2) * ZLENI
        R1  = PT2 + FACT * SAFZ2
        SAFR1 = (R - R1) / SQRT(1.0 + FACT*FACT)
 
C......  Compute radial distance to outer wall.
 
        FACT  = (PT5 - PT3) * ZLENI
        R2  = PT3 + FACT * SAFZ2
        SAFR2 = (R2 - R) / SQRT(1.0 + FACT*FACT)
        SAFER = MIN (SAFR1, SAFR2)
 
        IF (IND .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*********************************************************************
 
          PHI1  = PT6 * DEGRAD
          PHI2  = PT7 * DEGRAD
          IF (PHI2 .LT. PHI1)  PHI2 = PHI2 + TWOPI
 
          COSPH1  = COS (PHI1)
          COSPH2  = COS (PHI2)
          SINPH1  = SIN (PHI1)
          SINPH2  = SIN (PHI2)
 
C......  Get coordinates of outer endpoints (at R2) of both PHI-segments.
 
          XS1  = R2 * COSPH1
          YS1  = R2 * SINPH1
          XS2  = R2 * COSPH2
          YS2  = R2 * 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(X(1) * SINPH1 - X(2) * COSPH1)
          ELSE
            SAFSEG = ABS(X(1) * SINPH2 - X(2) * COSPH2)
          ENDIF
 
 
  100   CONTINUE
 
        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' (McPherson)  |
C     ------------------------------------------------
*
*     Avoid rounding effects induced by translation ********************
*
      IF (ABS(XT(3)).GE.DZ) XT(3) = (1.-0.000001)*XT(3)
*
      IF (PT2 .NE. PT4) GO TO 200
      IF (PT3 .NE. PT5) GO TO 200
 
      DELZ  = DZ
      PT(1) = PT2
      PT(2) = PT3
      PT(3) = DELZ
      PT(4) = PT6
      PT(5) = PT7
 
      CALL GNTUBE (XT, PT, 3, IND, SNEXT, SNXT, TSAFE)
      GO TO 999
 
  200 CONTINUE
      CALL GNCONE (XT, PT, 3, IND, SNEXT, SNXT, TSAFE)
 
  999 CONTINUE
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1995/10/24 10:20:53 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 GNPCON (X, PAR, IACT, SNEXT, SNXT, SAFE)
	13	C.
	14	C. ******************************************************************
	15	C. * *
	16	C. * COMPUTE DISTANCE UP TO INTERSECTION WITH 'PCON' 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	DIMENSION X(6), PAR(9), XT(6), PT(7)
	36
	37	EQUIVALENCE (PT(1), DZ), (PT(2), PT2), (PT(3), PT3)
	38	EQUIVALENCE (PT(4), PT4), (PT(5), PT5), (PT(6), PT6)
	39	EQUIVALENCE (PT(7), PT7)
	40
	41	EQUIVALENCE (XT(3), XT3)
	42	C.
	43	C. ---------------------------------------------------
	44	C.
	45
	46	SNXT = BIG
	47	R2 = X(1)X(1) + X(2)X(2)
	48	R = SQRT (R2)
	49
	50	NZ = PAR(3)
	51	ZMIN = PAR(4)
	52	ZMAX = PAR(3*NZ+1)
	53	SAFZ1 = X(3) - ZMIN
	54	SAFZ2 = ZMAX - X(3)
	55	SAFEZ = MIN (SAFZ1, SAFZ2)
	56
	57	C
	58	C...... First determine in which z-segment the particle is located.
	59	C
	60	DO 10 JPH=7, 3*(NZ-1)+1, 3
	61	IF (X(3) .LE. PAR(JPH)) THEN
	62	IPH=JPH
	63	GO TO 20
	64	ENDIF
65	10 CONTINUE
66	IPH=3*NZ+1
67
68	20 CONTINUE
69	C
70	C...... The particle is in the segment bounded by z-planes at
71	C...... Z1=PAR(IPL) and Z2=PAR(IPH), i.e., Z1 < X(3) < Z2.
72	C
73	C...... Set parameters for this segment and translate z-coordinate
74	C...... of point relative to center of this segment.This is done in
75	C...... preparation of invoking the algorithms used in "GNTUBE" and
76	C...... "GNCONE" (which for reasons of efficiency and clarity are
77	C...... implemented inline).
78	C
79	IPL = IPH - 3
80	DZ = 0.5 * (PAR(IPH) - PAR(IPL))
81	PT2 = PAR(IPL+1)
82	PT3 = PAR(IPL+2)
83	PT4 = PAR(IPH+1)
84	PT5 = PAR(IPH+2)
85	PT6 = PAR(1)
86	PT7 = PAR(1) + PAR(2)
87	IF (PT7 .GT. 360.0) PT7 = PT7 - 360.0
88
89	XT3 = X(3) - 0.5 * (PAR(IPL) + PAR(IPH))
90
91	XT(1) = X(1)
92	XT(2) = X(2)
93	XT(4) = X(4)
94	XT(5) = X(5)
95	XT(6) = X(6)
96
97	IND = 2
98	IF (PAR(2) .EQ. 360.0) IND = 1
99
100	IF (IACT .LT. 3) THEN
101
102	C -------------------------------------------------
103	C \| Compute safety-distance 'SAFE' (P.Weidhaas) \|
104	C -------------------------------------------------
105	SAFZ1 = DZ - ABS(XT3)
106	SAFEZ = MIN (SAFEZ,SAFZ1)
107	SAFSEG = BIG
108	C
109	C...... Determine whether the segment is a tube or a cone.
110	C
111
112	IF (PT2 .NE. PT4) GO TO 50
113	IF (PT3 .NE. PT5) GO TO 50
114
115	C*********************************************************
116	C
117	C...... The segment is a tube: invoke the algorithm
118	C...... from routine "GNTUBE" inline to get "SAFER".
119	C
120	C*********************************************************
121
122	SAFR1 = R - PT2
123	SAFR2 = PT3 - R
124	SAFER = MIN (SAFR1, SAFR2)
125
126	IF (IND .EQ. 2) GO TO 70
127
128	GO TO 100
129
130
131	50 CONTINUE
132
133	C*********************************************************
134	C
135	C...... The segment is a cone: invoke the algorithm
136	C...... from routine "GNCONE" inline to get "SAFER".
137	C
138	C*********************************************************
139
140	SAFZ2 = DZ + XT3
141	ZLENI = 0.5 / DZ
142
143	C...... Compute radial distance to inner wall.
144
145	FACT = (PT4 - PT2) * ZLENI
146	R1 = PT2 + FACT * SAFZ2
147	SAFR1 = (R - R1) / SQRT(1.0 + FACT*FACT)
148
149	C...... Compute radial distance to outer wall.
150
151	FACT = (PT5 - PT3) * ZLENI
152	R2 = PT3 + FACT * SAFZ2
153	SAFR2 = (R2 - R) / SQRT(1.0 + FACT*FACT)
154	SAFER = MIN (SAFR1, SAFR2)
155
156	IF (IND .EQ. 1) GO TO 100
157
158	70 CONTINUE
159
160	C********************************************************************
161	C...... Here we handle the case of a PHI-segment of a tube or cone.
162	C...... In addition to the radial distances (SAFR1, SAFR2) and the
163	C...... axial distances (SAFZ1, SAFZ2) we compute here the distance
164	C...... to the PHI-segment boundary that is closest to the point.
165	C
166	C...... For each PHI-boundary we find the distance from the given
167	C...... point to the outer (at R2) point of the segment boundary
168	C...... (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
169	C...... "SAFSEG" to be the distance to segment PHI1, else we set
170	C...... "SAFSEG" to be the distance to segment PHI2.
171	C*********************************************************************
172
173	PHI1 = PT6 * DEGRAD
174	PHI2 = PT7 * DEGRAD
175	IF (PHI2 .LT. PHI1) PHI2 = PHI2 + TWOPI
176
177	COSPH1 = COS (PHI1)
178	COSPH2 = COS (PHI2)
179	SINPH1 = SIN (PHI1)
180	SINPH2 = SIN (PHI2)
181
182	C...... Get coordinates of outer endpoints (at R2) of both PHI-segments.
183
184	XS1 = R2 * COSPH1
185	YS1 = R2 * SINPH1
186	XS2 = R2 * COSPH2
187	YS2 = R2 * SINPH2
188
189	C...... Get distances (squared) from given point to each endpoint.
190
191	DISTS1 = (X(1) - XS1)2 + (X(2) - YS1)2
192	DISTS2 = (X(1) - XS2)2 + (X(2) - YS2)2
193
194	C...... Get distance to that PHI-segment whose endpoint
195	C...... is closest to the given point.
196
197	IF (DISTS1 .LE. DISTS2) THEN
198	SAFSEG = ABS(X(1) * SINPH1 - X(2) * COSPH1)
199	ELSE
200	SAFSEG = ABS(X(1) * SINPH2 - X(2) * COSPH2)
201	ENDIF
202
203
204	100 CONTINUE
205
206	SAFE = MIN (SAFEZ, SAFER, SAFSEG)
207
208	IF (IACT .EQ. 0) GO TO 999
209	IF (IACT .EQ. 1) THEN
210	IF (SNEXT .LT. SAFE) GO TO 999
211	ENDIF
212	ENDIF
213
214
215	C ------------------------------------------------
216	C \| Compute vector-distance 'SNXT' (McPherson) \|
217	C ------------------------------------------------
218	*
219	* Avoid rounding effects induced by translation ********************
220	*
221	IF (ABS(XT(3)).GE.DZ) XT(3) = (1.-0.000001)*XT(3)
222	*
223	IF (PT2 .NE. PT4) GO TO 200
224	IF (PT3 .NE. PT5) GO TO 200
225
226	DELZ = DZ
227	PT(1) = PT2
228	PT(2) = PT3
229	PT(3) = DELZ
230	PT(4) = PT6
231	PT(5) = PT7
232
233	CALL GNTUBE (XT, PT, 3, IND, SNEXT, SNXT, TSAFE)
234	GO TO 999
235
236	200 CONTINUE
237	CALL GNCONE (XT, PT, 3, IND, SNEXT, SNXT, TSAFE)
238
239	999 CONTINUE
240	END
241