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