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