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