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