]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/ggeom/ggperp.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
100 parameters now allowed for geant shapes
[u/mrichter/AliRoot.git] / GEANT321 / ggeom / ggperp.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
d43b40e2 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:50 cernlib
9* Geant
10*
11*
12#include "geant321/pilot.h"
13*FCA : 05/01/99 09:58:02 by Federico Carminati
14* Effectively print the message when a shape is
15* not implemented
16*CMZ : 3.21/02 29/03/94 15.41.28 by S.Giani
17*-- Author :
18 SUBROUTINE GGPERP (X,U,IERR)
19C.
20C. ****************************************************************
21C. * *
22C. * This routine solves the general problem of calculating the *
23C. * unit vector normal to the surface of the current volume at *
24C. * the point X. The result is returned in the array U. X is *
25C. * assumed to be on or near a boundary of the current volume. *
26C. * The current volume is indicated by the common /GCVOLU/. *
27C. * U points from inside to outside in that neighbourhood. *
28C. * If X is equidistant to more than one boundary (in a corner) *
29C. * an arbitrary choice is made based upon the order of *
30C. * precedence implied by the IF statements below. If the *
31C. * routine fails to find the unit normal, it returns with *
32C. * IERR=1, otherwise IERR=0. *
33C. * *
34C. * Called by : GSURFP, GDSTEP *
35C. * Authors : F.Carminati, R.Jones, F.Ohlsson-Malek *
36C. * *
37C. ****************************************************************
38#include "geant321/gcvolu.inc"
39#include "geant321/gconsp.inc"
40#include "geant321/gcbank.inc"
41#include "geant321/gcshno.inc"
42#include "geant321/gctmed.inc"
43#include "geant321/gcunit.inc"
d43b40e2 44 DIMENSION X(3),U(3),XL(3),UL(3),DXL(3),PAR(100),SPAR(100),ATT(20)
fe4da5cc 45 DIMENSION PERP(10)
46#if !defined(CERNLIB_SINGLE)
47 DOUBLE PRECISION PERP,PMIN0
48 DOUBLE PRECISION PAR,DXL,RHO,R,RINV,PHI,THE
49 DOUBLE PRECISION PHI1,PHI2,THE1,THE2,XWID
50 DOUBLE PRECISION GUARD,DPHI,PHI0,SPHI0,CPHI0
51 DOUBLE PRECISION FACT,CALPH,SALPH,TALPH
52 DOUBLE PRECISION RAT,RATL,RATH,H,BL,TL,DX,DY,DZ,DU
53 DOUBLE PRECISION UU0,VV0,UU,W1,W2,W3,W4
54 DOUBLE PRECISION SEW1,SEW2,SEW3,SEW4
55 DOUBLE PRECISION TAN1,TAN2,TAN3,TAN4
56 DOUBLE PRECISION SEC1,SEC2,SEC3,SEC4
57 DOUBLE PRECISION U0,V0,U1,U1L,U2,U2L
58 DOUBLE PRECISION ONE,TWO
59 DOUBLE PRECISION DSECT,ZERO,FULL,FULL10,DBY2
60#endif
61 LOGICAL LNTFOU
62 PARAMETER (ONE=1,TWO=2)
63 PARAMETER (ZERO=0.,DBY2=0.5,FULL=360.,FULL10=3600.)
64C.
65C. ------------------------------------------------------------------
66C.
67 LNTFOU = .FALSE.
68*
69* *** Transform current point into local reference system
70 CALL GMTOD(X,XL,1)
71 DXL(1) = XL(1)
72 DXL(2) = XL(2)
73 DXL(3) = XL(3)
74*
75* *** Fetch the parameters of the current volume
76 JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
77 IN = LINDEX(NLEVEL)
78 IF (NLEVEL.GT.1) THEN
79 JVOM = LQ(JVOLUM-LVOLUM(NLEVEL-1))
80 JIN = LQ(JVOM-IN)
81 ENDIF
82 ISH = Q(JVO+2)
83 NIN = Q(JVO+3)
84 IF (NLEVEL.LT.NLDEV(NLEVEL)) THEN
85 JPAR = 0
86 ELSE
87* (case with structure JVOLUM locally developed)
88 JPAR = LQ(LQ(JVOLUM-LVOLUM(NLDEV(NLEVEL))))
89 IF (NLEVEL.EQ.NLDEV(NLEVEL)) GO TO 20
90 DO 10 ILEV = NLDEV(NLEVEL), NLEVEL-1
91 IF (IQ(JPAR+1).EQ.0) THEN
92 JPAR = LQ(JPAR-LINDEX(ILEV+1))
93 IF (JPAR.EQ.0) GO TO 20
94 ELSE IF (IQ(JPAR-3).GT.1) THEN
95 JPAR = LQ(JPAR-LINDEX(ILEV+1))
96 ELSE
97 JPAR = LQ(JPAR-1)
98 ENDIF
99 IF (ILEV.EQ.NLEVEL-1) THEN
100 JPAR = JPAR + 5
101 NPAR = IQ(JPAR)
102 CALL UCOPY (Q(JPAR+1), SPAR, NPAR)
103 DO 100 I=1,NPAR
104 PAR(I)=SPAR(I)
105 100 CONTINUE
106 ENDIF
107 10 CONTINUE
108 GO TO 30
109 ENDIF
110* (normal case)
111 20 CONTINUE
112 CALL GFIPAR(JVO,JIN,IN,NPAR,NATT,SPAR,ATT)
113 DO 101 I=1,NPAR
114 PAR(I)=SPAR(I)
115 101 CONTINUE
116 30 CONTINUE
117*
118* *** Case of the BOX:
119 IF (ISH.EQ.NSBOX) THEN
120 PERP(1) = ABS(ABS(DXL(1))-PAR(1))
121 PERP(2) = ABS(ABS(DXL(2))-PAR(2))
122 PERP(3) = ABS(ABS(DXL(3))-PAR(3))
123 PMIN0 = MIN(PERP(1),PERP(2),PERP(3))
124 IF (PERP(1).EQ.PMIN0) THEN
125 UL(1) = SIGN(ONE,DXL(1))
126 UL(2) = 0.
127 UL(3) = 0.
128 ELSE IF (PERP(2).EQ.PMIN0) THEN
129 UL(1) = 0.
130 UL(2) = SIGN(ONE,DXL(2))
131 UL(3) = 0.
132 ELSE IF (PERP(3).EQ.PMIN0) THEN
133 UL(1) = 0.
134 UL(2) = 0.
135 UL(3) = SIGN(ONE,DXL(3))
136 ELSE
137 LNTFOU=.TRUE.
138 ENDIF
139*
140* *** Case of the TUBE, TUBeSection:
141 ELSE IF (ISH.EQ.NSTUBE.OR.ISH.EQ.NSTUBS) THEN
142 RHO = SQRT(DXL(1)**2 + DXL(2)**2)
143 PERP(1) = ABS(RHO-PAR(1))
144 PERP(2) = ABS(RHO-PAR(2))
145 PERP(3) = ABS(ABS(DXL(3))-PAR(3))
146 IF (ISH.EQ.NSTUBE) THEN
147 PMIN0 = MIN(PERP(1),PERP(2),PERP(3))
148 ELSE
149 PHI = ATAN2(DXL(2),DXL(1))
150 IF (PHI.LT.0.) PHI = PHI+TWOPI
151 PHI1 = MOD(PAR(4)+FULL10,FULL)*DEGRAD
152 PERP(4) = ABS(PHI-PHI1)
153 IF (PERP(4).GT.PI) PERP(4) = TWOPI-PERP(4)
154 PHI2 = MOD(PAR(5)+FULL10,FULL)*DEGRAD
155 PERP(5) = ABS(PHI-PHI2)
156 IF (PERP(5).GT.PI) PERP(5) = TWOPI-PERP(5)
157 PERP(4) = PERP(4)*RHO
158 PERP(5) = PERP(5)*RHO
159 PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),PERP(5))
160 ENDIF
161 IF (PERP(1).EQ.PMIN0) THEN
162 UL(1) = -DXL(1)/RHO
163 UL(2) = -DXL(2)/RHO
164 UL(3) = 0.
165 ELSE IF (PERP(2).EQ.PMIN0) THEN
166 UL(1) = DXL(1)/RHO
167 UL(2) = DXL(2)/RHO
168 UL(3) = 0.
169 ELSE IF (PERP(3).EQ.PMIN0) THEN
170 UL(1) = 0.
171 UL(2) = 0.
172 UL(3) = SIGN(ONE,DXL(3))
173 ELSE IF (PERP(4).EQ.PMIN0) THEN
174 UL(1) = SIN(PHI1)
175 UL(2) = -COS(PHI1)
176 UL(3) = 0.
177 ELSE IF (PERP(5).EQ.PMIN0) THEN
178 UL(1) = -SIN(PHI2)
179 UL(2) = COS(PHI2)
180 UL(3) = 0.
181 ELSE
182 LNTFOU=.TRUE.
183 ENDIF
184*
185* *** Case of the CONE, CONeSection:
186 ELSE IF (ISH.EQ.NSCONE.OR.ISH.EQ.NSCONS) THEN
187 RHO = SQRT(DXL(1)**2 + DXL(2)**2)
188 TAN1 = (PAR(4)-PAR(2))/(TWO*PAR(1))
189 SEC1 = SQRT(ONE+TAN1**2)
190 U1 = RHO-DXL(3)*TAN1
191 U1L = PAR(4)-PAR(1)*TAN1
192 TAN2 = (PAR(5)-PAR(3))/(TWO*PAR(1))
193 SEC2 = SQRT(ONE+TAN2**2)
194 U2 = RHO-DXL(3)*TAN2
195 U2L = PAR(5)-PAR(1)*TAN2
196 PERP(1) = ABS(ABS(DXL(3))-PAR(1))
197 PERP(2) = ABS(U1-U1L)/SEC1
198 PERP(3) = ABS(U2-U2L)/SEC2
199 IF (ISH.EQ.NSCONE) THEN
200 PMIN0 = MIN(PERP(1),PERP(2),PERP(3))
201 ELSE
202 PHI = ATAN2(DXL(2),DXL(1))
203 IF (PHI.LT.0.) PHI = PHI+TWOPI
204 PHI1 = MOD(PAR(6)+FULL10,FULL)*DEGRAD
205 PERP(4) = ABS(PHI-PHI1)
206 IF (PERP(4).GT.PI) PERP(4) = TWOPI-PERP(4)
207 PHI2 = MOD(PAR(7)+FULL10,FULL)*DEGRAD
208 PERP(5) = ABS(PHI-PHI2)
209 IF (PERP(5).GT.PI) PERP(5) = TWOPI-PERP(5)
210 PERP(4) = PERP(4)*RHO
211 PERP(5) = PERP(5)*RHO
212 PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),PERP(5))
213 ENDIF
214 IF (PERP(1).EQ.PMIN0) THEN
215 UL(1) = 0.
216 UL(2) = 0.
217 UL(3) = SIGN(ONE,DXL(3))
218 ELSE IF (PERP(2).EQ.PMIN0) THEN
219 RHO = RHO*SEC1
220 UL(1) = -DXL(1)/RHO
221 UL(2) = -DXL(2)/RHO
222 UL(3) = TAN1/SEC1
223 ELSE IF (PERP(3).EQ.PMIN0) THEN
224 RHO = RHO*SEC2
225 UL(1) = DXL(1)/RHO
226 UL(2) = DXL(2)/RHO
227 UL(3) = -TAN2/SEC2
228 ELSE IF (PERP(4).EQ.PMIN0) THEN
229 UL(1) = SIN(PHI1)
230 UL(2) = -COS(PHI1)
231 UL(3) = 0.
232 ELSE IF (PERP(5).EQ.PMIN0) THEN
233 UL(1) = -SIN(PHI2)
234 UL(2) = COS(PHI2)
235 UL(3) = 0.
236 ELSE
237 LNTFOU=.TRUE.
238 ENDIF
239*
240* *** Case of the PolyCONe:
241 ELSE IF (ISH.EQ.NSPCON) THEN
242 PERP(1) = ABS(DXL(3)-PAR(4))
243 DO 400 I=7,NPAR,3
244 PERP(2) = ABS(DXL(3)-PAR(I))
245 IF (PERP(2).GT.PERP(1)) GOTO 401
246 PERP(1) = PERP(2)
247 400 CONTINUE
248 401 I = I-3
249 IF (I.GT.4) THEN
250 PERP(1) = 100.
251 RHO = SQRT(DXL(1)**2 + DXL(2)**2)
252 DZ = PAR(I)-PAR(I-3)+1.e-10
253 TAN1 = (PAR(I+1)-PAR(I-2))/DZ
254 SEC1 = SQRT(ONE+TAN1**2)
255 U1 = RHO-DXL(3)*TAN1
256 U1L = PAR(I+1)-PAR(I)*TAN1
257 TAN2 = (PAR(I+2)-PAR(I-1))/DZ
258 SEC2 = SQRT(ONE+TAN2**2)
259 U2 = RHO-DXL(3)*TAN2
260 U2L = PAR(I+2)-PAR(I)*TAN2
261 GUARD = MAX(DXL(3)-PAR(I),ZERO)
262 PERP(3) = ABS(U1-U1L)/SEC1 + GUARD*SEC1
263 PERP(4) = ABS(U2-U2L)/SEC2 + GUARD*SEC2
264 ELSE
265 PERP(3) = 100.
266 PERP(4) = 100.
267 ENDIF
268 IF (I.LT.NPAR-2) THEN
269 PERP(2) = 100.
270 RHO = SQRT(DXL(1)**2 + DXL(2)**2)
271 DZ = PAR(I+3)-PAR(I)+1.e-10
272 TAN3 = (PAR(I+4)-PAR(I+1))/DZ
273 SEC3 = SQRT(ONE+TAN3**2)
274 U1 = RHO-DXL(3)*TAN3
275 U1L = PAR(I+1)-PAR(I)*TAN3
276 TAN4 = (PAR(I+5)-PAR(I+2))/DZ
277 SEC4 = SQRT(ONE+TAN4**2)
278 U2 = RHO-DXL(3)*TAN4
279 U2L = PAR(I+2)-PAR(I)*TAN4
280 GUARD = MAX(PAR(I)-DXL(3),ZERO)
281 PERP(5) = ABS(U1-U1L)/SEC3 + GUARD*SEC3
282 PERP(6) = ABS(U2-U2L)/SEC4 + GUARD*SEC4
283 ELSE
284 PERP(5) = 100.
285 PERP(6) = 100.
286 ENDIF
287 PHI = ATAN2(DXL(2),DXL(1))
288 IF (PHI.LT.0.) PHI = PHI+TWOPI
289 PHI1 = MOD(PAR(1)+FULL10,FULL)*DEGRAD
290 PERP(7) = ABS(PHI-PHI1)
291 IF (PERP(7).GT.PI) PERP(7) = TWOPI-PERP(7)
292 PHI2 = MOD(PAR(1)+PAR(2)+FULL10,FULL)*DEGRAD
293 PERP(8) = ABS(PHI-PHI2)
294 IF (PERP(8).GT.PI) PERP(8) = TWOPI-PERP(8)
295 PERP(7) = PERP(7)*RHO
296 PERP(8) = PERP(8)*RHO
297 PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),
298 + PERP(5),PERP(6),PERP(7),PERP(8))
299 IF (PERP(1).EQ.PMIN0) THEN
300 UL(1) = 0.
301 UL(2) = 0.
302 UL(3) = -1.
303 ELSE IF (PERP(2).EQ.PMIN0) THEN
304 UL(1) = 0.
305 UL(2) = 0.
306 UL(3) = 1.
307 ELSE IF (PERP(3).EQ.PMIN0) THEN
308 RHO = RHO*SEC1
309 UL(1) = -DXL(1)/RHO
310 UL(2) = -DXL(2)/RHO
311 UL(3) = TAN1/SEC1
312 ELSE IF (PERP(4).EQ.PMIN0) THEN
313 RHO = RHO*SEC2
314 UL(1) = DXL(1)/RHO
315 UL(2) = DXL(2)/RHO
316 UL(3) = -TAN2/SEC2
317 ELSE IF (PERP(5).EQ.PMIN0) THEN
318 RHO = RHO*SEC3
319 UL(1) = -DXL(1)/RHO
320 UL(2) = -DXL(2)/RHO
321 UL(3) = TAN3/SEC3
322 ELSE IF (PERP(6).EQ.PMIN0) THEN
323 RHO = RHO*SEC4
324 UL(1) = DXL(1)/RHO
325 UL(2) = DXL(2)/RHO
326 UL(3) = -TAN4/SEC4
327 ELSE IF (PERP(7).EQ.PMIN0) THEN
328 UL(1) = SIN(PHI1)
329 UL(2) = -COS(PHI1)
330 UL(3) = 0.
331 ELSE IF (PERP(8).EQ.PMIN0) THEN
332 UL(1) = -SIN(PHI2)
333 UL(2) = COS(PHI2)
334 UL(3) = 0.
335 ELSE
336 LNTFOU=.TRUE.
337 ENDIF
338*
339* *** Case of the PolyGON:
340 ELSE IF (ISH.EQ.NSPGON) THEN
341 RHO = SQRT(DXL(1)**2+DXL(2)**2)
342 PHI = ATAN2(DXL(2),DXL(1))
343 IF (PHI.LT.0.) PHI = PHI+TWOPI
344 DPHI = MOD(PHI*RADDEG-PAR(1)+FULL10,FULL)
345 PDIV = PAR(2)/PAR(3)
346 DSECT = INT(DPHI/PDIV + ONE)
347 IF (DSECT.GT.PAR(3)) THEN
348 IF (DPHI.GT.(180.+PAR(2)*DBY2)) THEN
349 DSECT = ONE
350 ELSE
351 DSECT = PAR(3)
352 ENDIF
353 ENDIF
354 PHI0 = MOD(PAR(1)+(DSECT-DBY2)*PDIV+FULL10,FULL)*DEGRAD
355 SPHI0 = SIN(PHI0)
356 CPHI0 = COS(PHI0)
357 U0 = DXL(1)*CPHI0 + DXL(2)*SPHI0
358 V0 = DXL(2)*CPHI0 - DXL(1)*SPHI0
359 PERP(1) = ABS(DXL(3)-PAR(5))
360 DO 500 I=8,NPAR,3
361 PERP(2) = ABS(DXL(3)-PAR(I))
362 IF (PERP(2).GT.PERP(1)) GOTO 501
363 PERP(1) = PERP(2)
364 500 CONTINUE
365 501 I = I-3
366 IF (I.GT.5) THEN
367 PERP(1) = 100.
368 DZ = PAR(I)-PAR(I-3)+1.e-10
369 TAN1 = (PAR(I+1)-PAR(I-2))/DZ
370 SEC1 = SQRT(ONE+TAN1**2)
371 U1 = U0-DXL(3)*TAN1
372 U1L = PAR(I+1)-PAR(I)*TAN1
373 TAN2 = (PAR(I+2)-PAR(I-1))/DZ
374 SEC2 = SQRT(ONE+TAN2**2)
375 U2 = U0-DXL(3)*TAN2
376 U2L = PAR(I+2)-PAR(I)*TAN2
377 GUARD = MAX(DXL(3)-PAR(I),ZERO)
378 PERP(3) = ABS(U1-U1L)/SEC1 + GUARD*SEC1
379 PERP(4) = ABS(U2-U2L)/SEC2 + GUARD*SEC2
380 ELSE
381 PERP(3) = 100.
382 PERP(4) = 100.
383 ENDIF
384 IF (I.LT.NPAR-2) THEN
385 PERP(2) = 100.
386 DZ = PAR(I+3)-PAR(I)+1.e-10
387 TAN3 = (PAR(I+4)-PAR(I+1))/DZ
388 SEC3 = SQRT(ONE+TAN3**2)
389 U1 = U0-DXL(3)*TAN3
390 U1L = PAR(I+1)-PAR(I)*TAN3
391 TAN4 = (PAR(I+5)-PAR(I+2))/DZ
392 SEC4 = SQRT(ONE+TAN4**2)
393 U2 = U0-DXL(3)*TAN4
394 U2L = PAR(I+2)-PAR(I)*TAN4
395 GUARD = MAX(PAR(I)-DXL(3),ZERO)
396 PERP(5) = ABS(U1-U1L)/SEC3 + GUARD*SEC3
397 PERP(6) = ABS(U2-U2L)/SEC4 + GUARD*SEC4
398 ELSE
399 PERP(5) = 100.
400 PERP(6) = 100.
401 ENDIF
402 PHI1 = MOD(PAR(1)+FULL10,FULL)*DEGRAD
403 PERP(7) = ABS(PHI-PHI1)
404 IF (PERP(7).GT.PI) PERP(7) = TWOPI-PERP(7)
405 PHI2 = MOD(PAR(1)+PAR(2)+FULL10,FULL)*DEGRAD
406 PERP(8) = ABS(PHI-PHI2)
407 IF (PERP(8).GT.PI) PERP(8) = TWOPI-PERP(8)
408 PERP(7) = PERP(7)*RHO
409 PERP(8) = PERP(8)*RHO
410 PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),
411 + PERP(5),PERP(6),PERP(7),PERP(8))
412 IF (PERP(1).EQ.PMIN0) THEN
413 UL(1) = 0.
414 UL(2) = 0.
415 UL(3) = -1.
416 ELSE IF (PERP(2).EQ.PMIN0) THEN
417 UL(1) = 0.
418 UL(2) = 0.
419 UL(3) = 1.
420 ELSE IF (PERP(3).EQ.PMIN0) THEN
421 FACT = ONE/SEC1
422 UL(1) = -CPHI0*FACT
423 UL(2) = -SPHI0*FACT
424 UL(3) = TAN1*FACT
425 ELSE IF (PERP(4).EQ.PMIN0) THEN
426 FACT = ONE/SEC2
427 UL(1) = CPHI0*FACT
428 UL(2) = SPHI0*FACT
429 UL(3) = -TAN2*FACT
430 ELSE IF (PERP(5).EQ.PMIN0) THEN
431 FACT = ONE/SEC3
432 UL(1) = -CPHI0*FACT
433 UL(2) = -SPHI0*FACT
434 UL(3) = TAN3*FACT
435 ELSE IF (PERP(6).EQ.PMIN0) THEN
436 FACT = ONE/SEC4
437 UL(1) = CPHI0*FACT
438 UL(2) = SPHI0*FACT
439 UL(3) = -TAN4*FACT
440 ELSE IF (PERP(7).EQ.PMIN0) THEN
441 UL(1) = SIN(PHI1)
442 UL(2) = -COS(PHI1)
443 UL(3) = 0.
444 ELSE IF (PERP(8).EQ.PMIN0) THEN
445 UL(1) = -SIN(PHI2)
446 UL(2) = COS(PHI2)
447 UL(3) = 0.
448 ELSE
449 LNTFOU=.TRUE.
450 ENDIF
451*
452* *** Case of the SPHEre:
453 ELSE IF (ISH.EQ.NSSPHE) THEN
454 R = SQRT(DXL(1)**2+DXL(2)**2+DXL(3)**2)
455 RHO = SQRT(DXL(1)**2+DXL(2)**2)
456 THE = ATAN2(RHO,DXL(3))
457 PHI = ATAN2(DXL(2),DXL(1))
458 IF (PHI.LT.0.) PHI = PHI+TWOPI
459 THE1 = MOD(PAR(3)+FULL10,FULL)*DEGRAD
460 THE2 = MOD(PAR(4)+FULL10,FULL)*DEGRAD
461 PHI1 = MOD(PAR(5)+FULL10,FULL)*DEGRAD
462 PHI2 = MOD(PAR(6)+FULL10,FULL)*DEGRAD
463 PERP(1) = ABS(R-PAR(1))
464 PERP(2) = ABS(R-PAR(2))
465 PERP(3) = ABS(THE-THE1)*R
466 PERP(4) = ABS(THE-THE2)*R
467 PERP(5) = ABS(PHI-PHI1)
468 IF (PERP(5).GT.PI) PERP(5) = TWOPI-PERP(5)
469 PERP(5) = PERP(5)*RHO
470 PERP(6) = ABS(PHI-PHI2)
471 IF (PERP(6).GT.PI) PERP(6) = TWOPI-PERP(6)
472 PERP(6) = PERP(6)*RHO
473 PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),PERP(5),PERP(6))
474 IF (PERP(1).EQ.PMIN0) THEN
475 RINV = ONE/R
476 UL(1) = -DXL(1)*RINV
477 UL(2) = -DXL(2)*RINV
478 UL(3) = -DXL(3)*RINV
479 ELSE IF (PERP(2).EQ.PMIN0) THEN
480 RINV = ONE/R
481 UL(1) = DXL(1)*RINV
482 UL(2) = DXL(2)*RINV
483 UL(3) = DXL(3)*RINV
484 ELSE IF (PERP(3).EQ.PMIN0) THEN
485 UL(1) = -COS(THE1)*COS(PHI)
486 UL(2) = -COS(THE1)*SIN(PHI)
487 UL(3) = +SIN(THE1)
488 ELSE IF (PERP(4).EQ.PMIN0) THEN
489 UL(1) = +COS(THE2)*COS(PHI)
490 UL(2) = +COS(THE2)*SIN(PHI)
491 UL(3) = -SIN(THE2)
492 ELSE IF (PERP(5).EQ.PMIN0) THEN
493 UL(1) = +SIN(PHI1)
494 UL(2) = -COS(PHI1)
495 UL(3) = 0
496 ELSE IF (PERP(6).EQ.PMIN0) THEN
497 UL(1) = -SIN(PHI2)
498 UL(2) = +COS(PHI2)
499 UL(3) = 0
500 ELSE
501 LNTFOU=.TRUE.
502 ENDIF
503*
504* *** Case of the PARAllelpiped:
505***************************************************************
506* Warning: the parameters for this shape are NOT stored in *
507* the data structure as the user supplies them. Rather, the *
508* user supplies PAR(4)=alph, PAR(5)=the, PAR(6)=phi, and the *
509* data structure contains PAR(4)=Tan(alph), PAR(5)=Tan(the)* *
510* Cos(phi), PAR(6)=Tan(the)*Sin(phi). *
511***************************************************************
512 ELSE IF (ISH.EQ.NSPARA) THEN
513 DX = PAR(5)
514 DY = PAR(6)
515 U0 = DXL(1)-DX*DXL(3)
516 V0 = DXL(2)-DY*DXL(3)
517 CALPH = ONE/SQRT(ONE+PAR(4)**2)
518 SALPH = -CALPH*PAR(4)
519 U1 = U0*CALPH+V0*SALPH
520 U1L = PAR(1)*CALPH
521 PERP(1) = ABS(ABS(U1)-U1L)
522 PERP(2) = ABS(ABS(V0)-PAR(2))
523 PERP(3) = ABS(ABS(DXL(3))-PAR(3))
524 PMIN0 = MIN(PERP(1),PERP(2),PERP(3))
525 IF (PERP(1).EQ.PMIN0) THEN
526 DU = DX*CALPH+DY*SALPH
527 FACT = SIGN(ONE/SQRT(ONE+DU**2),U1)
528 UL(1) = CALPH*FACT
529 UL(2) = SALPH*FACT
530 UL(3) = -DU*FACT
531 ELSE IF (PERP(2).EQ.PMIN0) THEN
532 FACT = SIGN(ONE/SQRT(ONE+DY**2),V0)
533 UL(1) = 0.
534 UL(2) = FACT
535 UL(3) = -DY*FACT
536 ELSE IF (PERP(3).EQ.PMIN0) THEN
537 UL(1) = 0.
538 UL(2) = 0.
539 UL(3) = SIGN(ONE,DXL(3))
540 ELSE
541 LNTFOU=.TRUE.
542 ENDIF
543*
544* *** Case of the trapezoid TRD1
545 ELSE IF (ISH.EQ.NSTRD1) THEN
546 DZ = TWO*PAR(4)+1.e-10
547 TAN1 = (PAR(2)-PAR(1))/DZ
548 SEC1 = SQRT(ONE+TAN1**2)
549 U1 = ABS(DXL(1))-DXL(3)*TAN1
550 U1L = PAR(2)-PAR(4)*TAN1
551 PERP(1) = ABS(U1-U1L)/SEC1
552 PERP(2) = ABS(ABS(DXL(2))-PAR(3))
553 PERP(3) = ABS(ABS(DXL(3))-PAR(4))
554 PMIN0 = MIN(PERP(1),PERP(2),PERP(3))
555 IF (PERP(1).EQ.PMIN0) THEN
556 FACT = ONE/SEC1
557 UL(1) = SIGN(FACT,DXL(1))
558 UL(2) = 0.
559 UL(3) = -TAN1*FACT
560 ELSE IF (PERP(2).EQ.PMIN0) THEN
561 UL(1) = 0.
562 UL(2) = SIGN(ONE,DXL(2))
563 UL(3) = 0.
564 ELSE IF (PERP(3).EQ.PMIN0) THEN
565 UL(1) = 0.
566 UL(2) = 0.
567 UL(3) = SIGN(ONE,DXL(3))
568 ELSE
569 LNTFOU=.TRUE.
570 ENDIF
571*
572* *** Case of the trapezoid TRD2
573 ELSE IF (ISH.EQ.NSTRD2) THEN
574 DZ = TWO*PAR(5)+1.e-10
575 TAN1 = (PAR(2)-PAR(1))/DZ
576 SEC1 = SQRT(ONE+TAN1**2)
577 U1 = ABS(DXL(1))-DXL(3)*TAN1
578 U1L = PAR(2)-PAR(5)*TAN1
579 TAN2 = (PAR(4)-PAR(3))/DZ
580 SEC2 = SQRT(ONE+TAN2**2)
581 U2 = ABS(DXL(2))-DXL(3)*TAN2
582 U2L = PAR(4)-PAR(5)*TAN2
583 PERP(1) = ABS(U1-U1L)/SEC1
584 PERP(2) = ABS(U2-U2L)/SEC2
585 PERP(3) = ABS(ABS(DXL(3))-PAR(5))
586 PMIN0 = MIN(PERP(1),PERP(2),PERP(3))
587 IF (PERP(1).EQ.PMIN0) THEN
588 FACT = ONE/SEC1
589 UL(1) = SIGN(FACT,DXL(1))
590 UL(2) = 0.
591 UL(3) = -TAN1*FACT
592 ELSE IF (PERP(2).EQ.PMIN0) THEN
593 FACT = ONE/SEC2
594 UL(1) = 0.
595 UL(2) = SIGN(FACT,DXL(2))
596 UL(3) = -TAN2*FACT
597 ELSE IF (PERP(3).EQ.PMIN0) THEN
598 UL(1) = 0.
599 UL(2) = 0.
600 UL(3) = SIGN(ONE,DXL(3))
601 ELSE
602 LNTFOU=.TRUE.
603 ENDIF
604*
605* *** Case of the TRAPezoid
606***************************************************************
607* Warning: the parameters for this shape are NOT stored in *
608* the data structure as the user supplies them. Rather, the *
609* user supplies PAR(2)=thet, PAR(3)=phi, PAR(7)=alp1, and *
610* PAR(11)=alp2, while the data structure contains PAR(2)= *
611* Tan(thet)*Cos(phi), PAR(3)=Tan(thet)*Sin(phi), PAR(7)= *
612* Tan(alp1), and PAR(11)=Tan(alp2). *
613***************************************************************
614 ELSE IF (ISH.EQ.NSTRAP) THEN
615 PERP(1) = ABS(ABS(DXL(3))-PAR(1))
616 DX = PAR(2)
617 DY = PAR(3)
618 U0 = DX*DXL(3)
619 V0 = DY*DXL(3)
620 UU0 = DX*PAR(1)
621 VV0 = DY*PAR(1)
622 RAT = DXL(3)/PAR(1)
623 RATL = (ONE-RAT)/TWO
624 RATH = (ONE+RAT)/TWO
625 H = PAR(4)*RATL+PAR(8)*RATH
626 BL = PAR(5)*RATL+PAR(9)*RATH
627 TL = PAR(6)*RATL+PAR(10)*RATH
628 TALPH = PAR(7)*RATL+PAR(11)*RATH
629 XWID = (TL+BL)/TWO
630 TAN1 = TALPH+(TL-BL)/(TWO*H)
631 SEC1 = SQRT(ONE+TAN1**2)
632 U1 = DXL(1)-DXL(2)*TAN1
633 U1L = U0+XWID-V0*TAN1
634 TAN2 = TAN1-TWO*TALPH
635 SEC2 = SQRT(ONE+TAN2**2)
636 U2 = DXL(1)+DXL(2)*TAN2
637 U2L = U0-XWID+V0*TAN2
638 IF (DXL(3).LT.0) THEN
639 DZ = PAR(1)-DXL(3)+1.e-10
640 UU = UU0+(PAR(9)+PAR(10))/TWO
641 W1 = (UU-VV0*TAN1-U1L)/DZ
642 UU = TWO*UU0-UU
643 W2 = (UU+VV0*TAN2-U2L)/DZ
644 ELSE
645 DZ = -PAR(1)-DXL(3)+1.e-10
646 UU = -UU0+(PAR(5)+PAR(6))/TWO
647 W1 = (UU+VV0*TAN1-U1L)/DZ
648 UU = -TWO*UU0-UU
649 W2 = (UU-VV0*TAN2-U2L)/DZ
650 ENDIF
651 W3 = DY+(PAR(8)-PAR(4))/(TWO*PAR(1))
652 W4 = TWO*DY-W3
653 SEW1 = SQRT(ONE+W1**2)
654 SEW2 = SQRT(ONE+W2**2)
655 SEW3 = SQRT(ONE+W3**2)
656 SEW4 = SQRT(ONE+W4**2)
657 PERP(2) = ABS(U1-U1L)/(SEC1*SEW1)
658 PERP(3) = ABS(U2-U2L)/(SEC2*SEW2)
659 PERP(4) = ABS(DXL(2)-V0-H)/SEW3
660 PERP(5) = ABS(DXL(2)-V0+H)/SEW4
661 PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),PERP(5))
662 IF (PERP(1).EQ.PMIN0) THEN
663 UL(1) = 0.
664 UL(2) = 0.
665 UL(3) = SIGN(ONE,DXL(3))
666 ELSE IF (PERP(2).EQ.PMIN0) THEN
667 FACT = ONE/(SEC1*SEW1)
668 UL(1) = FACT
669 UL(2) = -TAN1*FACT
670 UL(3) = -W1/SEW1
671 ELSE IF (PERP(3).EQ.PMIN0) THEN
672 FACT = ONE/(SEC2*SEW2)
673 UL(1) = -FACT
674 UL(2) = -TAN2*FACT
675 UL(3) = W2/SEW2
676 ELSE IF (PERP(4).EQ.PMIN0) THEN
677 FACT = ONE/SEW3
678 UL(1) = 0.
679 UL(2) = FACT
680 UL(3) = -W3*FACT
681 ELSE IF (PERP(5).EQ.PMIN0) THEN
682 FACT = ONE/SEW4
683 UL(1) = 0.
684 UL(2) = -FACT
685 UL(3) = W4*FACT
686 ELSE
687 LNTFOU=.TRUE.
688 ENDIF
689*
690* *** everything else (currently NOT IMPLEMENTED)
691 ELSE
692 WRITE(CHMAIL,10100) ISH
693 CALL GMAIL(0,0)
694 IERR = 1
695 GOTO 999
696 ENDIF
697
698 IF(LNTFOU) THEN
699 WRITE(CHMAIL,10000) ISH
700 CALL GMAIL(0,0)
701 IERR = 1
702 ELSE
703*
704* *** Transform back into the MCS
705 CALL GDTOM(UL,U,2)
706 IERR = 0
707 ENDIF
708
70910000 FORMAT(' GGPERP - geometry check error for shape #',I2,'!')
71010100 FORMAT(' GGPERP - non implemented for shape #',I2)
711 999 END