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