This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / ggeom / ggperp.F.ori
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)
13 C.
14 C.    ****************************************************************
15 C.    *                                                              *
16 C.    *  This routine solves the general problem of calculating the  *
17 C.    *  unit vector normal to the surface of the current volume at  *
18 C.    *  the point X. The result is returned in the array U.  X is   *
19 C.    *  assumed to be on or near a boundary of the current volume.  *
20 C.    *  The current volume is indicated by the common /GCVOLU/.     *
21 C.    *  U points from inside to outside in that neighbourhood.      *
22 C.    *  If X is equidistant to more than one boundary (in a corner) *
23 C.    *  an arbitrary choice is made based upon the order of         *
24 C.    *  precedence implied by the IF statements below.  If the      *
25 C.    *  routine fails to find the unit normal, it returns with      *
26 C.    *  IERR=1, otherwise IERR=0.                                   *
27 C.    *                                                              *
28 C.    *   Called by : GSURFP, GDSTEP                                 *
29 C.    *   Authors   : F.Carminati, R.Jones, F.Ohlsson-Malek          *
30 C.    *                                                              *
31 C.    ****************************************************************
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.)
58 C.
59 C.    ------------------------------------------------------------------
60 C.
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  
702 10000 FORMAT(' GGPERP - geometry check error for shape #',I2,'!')
703 10100 FORMAT(' GGPERP - non implemented for shape #',I2)
704   999 END