GEANT321/ggeom/ggperp.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1999/05/18 15:55:17  fca
   6 * AliRoot sources
   7 *
   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)
  19 C.
  20 C.    ****************************************************************
  21 C.    *                                                              *
  22 C.    *  This routine solves the general problem of calculating the  *
  23 C.    *  unit vector normal to the surface of the current volume at  *
  24 C.    *  the point X. The result is returned in the array U.  X is   *
  25 C.    *  assumed to be on or near a boundary of the current volume.  *
  26 C.    *  The current volume is indicated by the common /GCVOLU/.     *
  27 C.    *  U points from inside to outside in that neighbourhood.      *
  28 C.    *  If X is equidistant to more than one boundary (in a corner) *
  29 C.    *  an arbitrary choice is made based upon the order of         *
  30 C.    *  precedence implied by the IF statements below.  If the      *
  31 C.    *  routine fails to find the unit normal, it returns with      *
  32 C.    *  IERR=1, otherwise IERR=0.                                   *
  33 C.    *                                                              *
  34 C.    *   Called by : GSURFP, GDSTEP                                 *
  35 C.    *   Authors   : F.Carminati, R.Jones, F.Ohlsson-Malek          *
  36 C.    *                                                              *
  37 C.    ****************************************************************
  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"
  44       DIMENSION X(3),U(3),XL(3),UL(3),DXL(3),PAR(100),SPAR(100),ATT(20)
  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.)
  64 C.
  65 C.    ------------------------------------------------------------------
  66 C.
  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
 709 10000 FORMAT(' GGPERP - geometry check error for shape #',I2,'!')
 710 10100 FORMAT(' GGPERP - non implemented for shape #',I2)
 711   999 END