GEANT321/ggeom/gnopco.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:52  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.29  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GNOPCO (X, PAR, IACT, SNEXT, SNXT, SAFE)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *       COMPUTE DISTANCE UP TO INTERSECTION WITH 'PCON' VOLUME,  *
  17 C.    *       FROM OUTSIDE POINT X(1-3) ALONG DIRECTION X(4-6)         *
  18 C.    *                                                                *
  19 C.    *       PAR   (input)  : volume parameters                       *
  20 C.    *       IACT  (input)  : action flag                             *
  21 C.    *         = 0  Compute SAFE only                                 *
  22 C.    *         = 1  Compute SAFE, and SNXT only if SNEXT .GT.new SAFE *
  23 C.    *         = 2  Compute both SAFE and SNXT                        *
  24 C.    *         = 3  Compute SNXT only                                 *
  25 C.    *       SNEXT (input)  : see IACT = 1                            *
  26 C.    *       SNXT  (output) : distance to volume boundary             *
  27 C.    *       SAFE  (output) : shortest distance to any boundary       *
  28 C.    *                                                                *
  29 C.    *    ==>Called by : GNEXT, GTNEXT                                *
  30 C.    *         Author  A.McPherson,  P.Weidhaas  *********            *
  31 C.    *                                                                *
  32 C.    ******************************************************************
  33 C.
  34 #include "geant321/gconsp.inc"
  35       DIMENSION X(6), PAR(9), PT(7), XT(6)
  36
  37       EQUIVALENCE (PT(1),PT1), (PT(2),PT2), (PT(3),PT3), (PT(4),PT4)
  38       EQUIVALENCE (PT(5),PT5), (PT(6),PT6), (PT(7),PT7)
  39 C.
  40 C.   --------------------------------------------------------
  41 C.
  42
  43       SNXT = BIG
  44       SAFE = BIG
  45       R2 = X(1)*X(1) + X(2)*X(2)
  46       R  = SQRT (R2)
  47       NZ  = PAR(3)
  48       NSEC= NZ - 1
  49       PT6=PAR(1)
  50       PT7=PAR(1)+PAR(2)
  51
  52       IFLAG = 2
  53       IF(PAR(2).EQ.360.0) IFLAG = 1
  54
  55       IF (IACT .LT. 3) THEN
  56
  57 C       -------------------------------------------------
  58 C       |  Compute safety-distance 'SAFE' (P.Weidhaas)  |
  59 C       -------------------------------------------------
  60
  61
  62         SAFEZ  = 0.0
  63 C
  64 C......  Obtain axial distance "SAFEZ".
  65 C
  66         ZMIN   = PAR(4)
  67         ZMAX   = PAR(3*NZ+1)
  68
  69         IF (X(3) .LT. ZMIN) THEN
  70           SAFEZ  = ZMIN - X(3)
  71         ELSEIF (X(3) .GT. ZMAX) THEN
  72           SAFEZ  = X(3) - ZMAX
  73         ENDIF
  74
  75
  76 C----------------------------------------------------
  77 C......  Prepare parameters for PHI-segmented cone:
  78 C----------------------------------------------------
  79
  80       IF (IFLAG .EQ. 2) THEN
  81
  82         PHI1 = MOD (PT6+360.0 , 360.0)
  83         PHI2 = MOD (PT7+360.0 , 360.0)
  84         IF ( X(1).NE.0.0 .OR. X(2).NE.0.0 ) THEN
  85            PHI = ATAN2( X(2), X(1) ) * RADDEG
  86         ELSE
  87            PHI = 0.0
  88         ENDIF
  89         PHI  = MOD (PHI+360.0 , 360.0)
  90
  91         SINPH1 = SIN(PHI1*DEGRAD)
  92         COSPH1 = COS(PHI1*DEGRAD)
  93         SINPH2 = SIN(PHI2*DEGRAD)
  94         COSPH2 = COS(PHI2*DEGRAD)
  95
  96
  97 C......  Set flag "IOPENP" if point (X(1),X(2),X(3)) lies in the
  98 C......  PHI-opening.
  99
 100         IOPENP = 0
 101         IF (PHI2 .GT. PHI1) THEN
 102           IF (PHI.GT.PHI2 .OR. PHI.LT.PHI1) IOPENP = 1
 103         ELSE
 104           IF (PHI.GT.PHI2 .AND. PHI.LT.PHI1) IOPENP = 1
 105         ENDIF
 106       ENDIF
 107
 108
 109
 110 C------------------   Start of loop over Z-sections  --------------
 111
 112         IPZ2=4
 113
 114         DO 150 IS=1,NSEC
 115
 116         IPZ1=IPZ2
 117         IPZ2=IPZ1+3
 118
 119         SAFSEG = 0.0
 120         SAFER  = 0.0
 121
 122         XT3=X(3) - 0.5 * (PAR(IPZ1)+PAR(IPZ2))
 123         DZ   = 0.5 * (PAR(IPZ2)-PAR(IPZ1))
 124
 125         PT2  = PAR(IPZ1+1)
 126         PT3  = PAR(IPZ1+2)
 127         PT4  = PAR(IPZ2+1)
 128         PT5  = PAR(IPZ2+2)
 129 C**** check DZ=0 segments
 130         IF (DZ.LE.0.) THEN
 131            IF ((R-PT2)*(R-PT4).LE.0. .OR. (R-PT3)*(R-PT5).LE.0.) THEN
 132               SAFER = ABS(XT3)
 133            ELSE
 134               GO TO 150
 135            ENDIF
 136            GO TO 100
 137         ENDIF
 138
 139         IF (PT2 .NE. PT4) GO TO 50
 140         IF (PT3 .NE. PT5) GO TO 50
 141
 142 C**********************************************************
 143 C
 144 C......  The segment is a tube;  invoke the algorithm
 145 C......  from "GNOTUB" inline to get "SAFER" and "SAFSEG".
 146 C
 147 C**********************************************************
 148
 149         IF (R .LT. PT2) SAFER = PT2 - R
 150         IF (R .GT. PT3) SAFER = R - PT3
 151
 152         IF (IFLAG .EQ. 2)  THEN
 153
 154 C********************************************************************
 155 C......  Handle the case in which we have a PHI-segment of a tube.
 156 C......  In addition to the radial distance (SAFER) and the
 157 C......  axial distance (SAFEZ) we compute here the distance (SAFSEG)
 158 C......  to the PHI-segment boundary that is closest to the point:
 159 C
 160 C......  SAFSEG is only calculated if PHI lies outside the interval
 161 C......  [PHI1, PHI2]. Here PHI is the angle to the given point
 162 C......  (thus we only consider SAFSEG if the point is outside the
 163 C......  PHI-segment).
 164 C
 165 C......  Algorithm to find SAFSEG (same as in routine "GNTUBE"):
 166 C
 167 C......  For each PHI-boundary we find the distance from the given
 168 C......  point to the outer (at RMAX) point of the segment boundary
 169 C......  (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
 170 C......  SAFSEG to be the distance to segment PHI1; else we set
 171 C......  SAFSEG to be the distance to segment PHI2.
 172 C********************************************************************
 173
 174
 175 C......  Next eliminate those points whose angle PHI places them
 176 C......  inside the given PHI-segment (IOPENP = 0).
 177
 178           IF (IOPENP .EQ. 0) GO TO 100
 179
 180
 181 C......  Get coordinates of outer endpoints (at RMAX) of both PHI-segments.
 182
 183           XS1  = PT3 * COSPH1
 184           YS1  = PT3 * SINPH1
 185           XS2  = PT3 * COSPH2
 186           YS2  = PT3 * SINPH2
 187
 188 C......  Get distances (squared) from the given point to each endpoint.
 189
 190           DISTS1 = (X(1) - XS1)**2  +  (X(2) - YS1)**2
 191           DISTS2 = (X(1) - XS2)**2  +  (X(2) - YS2)**2
 192
 193 C......  Get distance to that PHI-segment whose endpoint
 194 C......  is closest to the given point.
 195
 196           IF (DISTS1 .LE. DISTS2) THEN
 197             SAFSEG = ABS(X(1) * SINPH1 - X(2) * COSPH1)
 198           ELSE
 199             SAFSEG = ABS(X(1) * SINPH2 - X(2) * COSPH2)
 200           ENDIF
 201
 202         ENDIF
 203
 204         GO TO 100
 205
 206
 207    50   CONTINUE
 208
 209 C*********************************************************
 210 C
 211 C......  The segment is a cone; invoke the algorithm
 212 C......  from "GNOCON" inline to get "SAFER" and "SAFSEG".
 213 C
 214 C*********************************************************
 215
 216         ZLENI = 0.5 / DZ
 217
 218         FACT1  = (PT4 - PT2) * ZLENI
 219         FACT2  = (PT5 - PT3) * ZLENI
 220         RIN   = PT2 + FACT1 * (DZ + XT3)
 221         ROUT  = PT3 + FACT2 * (DZ + XT3)
 222
 223         IF (R .LT. RIN) THEN
 224           SAFER = (RIN - R) / SQRT(1.0 + FACT1 * FACT1)
 225         ELSE
 226           IF (R .GT. ROUT) THEN
 227             SAFER = (R - ROUT) / SQRT(1.0 + FACT2 * FACT2)
 228           ENDIF
 229         ENDIF
 230
 231
 232         IF (IFLAG .EQ. 2)  THEN
 233
 234 C********************************************************************
 235 C......  Handle the case in which we have a PHI-segment of a cone.
 236 C......  In addition to the radial distance (SAFER) and the
 237 C......  axial distance (SAFEZ) we compute here the distance (SAFSEG)
 238 C......  to the PHI-segment boundary that is closest to the point:
 239 C
 240 C......  SAFSEG is only calculated if PHI lies outside the interval
 241 C......  [PHI1, PHI2]. Here PHI is the angle to the given point
 242 C......  (thus we only consider SAFSEG if the point is outside the
 243 C......  PHI-segment).
 244 C
 245 C......  Algorithm to find SAFSEG (same as in routine "GNTUBE"):
 246 C
 247 C......  For each PHI-boundary we find the distance from the given
 248 C......  point to the outer (at ROUT) point of the segment boundary
 249 C......  (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
 250 C......  SAFSEG to be the distance to segment PHI1; else we set
 251 C......  SAFSEG to be the distance to segment PHI2.
 252 C********************************************************************
 253
 254
 255 C......  Next eliminate those points whose angle PHI places them
 256 C......  inside the given PHI-segment (IOPENP = 0).
 257
 258            IF (IOPENP .EQ. 0) GO TO 100
 259
 260 C......  Get coordinates of outer endpoints (at ROUT) of both PHI-segments.
 261
 262           IF (XT3 .LT. -DZ) THEN
 263             ROUT = PT3
 264           ELSEIF (XT3 .GT. DZ) THEN
 265             ROUT = PT5
 266           ENDIF
 267
 268           XS1  = ROUT * COSPH1
 269           YS1  = ROUT * SINPH1
 270           XS2  = ROUT * COSPH2
 271           YS2  = ROUT * SINPH2
 272
 273 C......  Get distances (squared) from the given point to each endpoint.
 274
 275           DISTS1 = (X(1) - XS1)**2  +  (X(2) - YS1)**2
 276           DISTS2 = (X(1) - XS2)**2  +  (X(2) - YS2)**2
 277
 278 C......  Obtain distance to that PHI-segment whose endpoint
 279 C......  is closest to the given point.
 280
 281           IF (DISTS1 .LE. DISTS2) THEN
 282             SAFSEG = ABS(X(1) * SINPH1 - X(2) * COSPH1)
 283           ELSE
 284             SAFSEG = ABS(X(1) * SINPH2 - X(2) * COSPH2)
 285           ENDIF
 286
 287         ENDIF
 288
 289
 290   100   CONTINUE
 291         TSAFE  = MAX (SAFEZ, SAFER, SAFSEG)
 292
 293         IF (TSAFE .GT. 0.0) THEN
 294           IF (TSAFE .LT. SAFE) SAFE  = TSAFE
 295         ENDIF
 296
 297         IF (TSAFE .EQ. 0.0) THEN
 298           IF (X(3) .LT. PAR(IPZ1)) GO TO 200
 299         ENDIF
 300   150   CONTINUE
 301
 302   200   CONTINUE
 303         IF (IACT .EQ. 0) GO TO 999
 304         IF (IACT .EQ. 1) THEN
 305           IF (SNEXT .LT. SAFE) GO TO 999
 306         ENDIF
 307       ENDIF
 308
 309
 310 C     ------------------------------------------------
 311 C     |  Compute vector-distance 'SNXT' (McPherson)  |
 312 C     ------------------------------------------------
 313
 314       TSNXT = BIG
 315
 316       DO 210 I=1, 6
 317         XT(I) = X(I)
 318   210 CONTINUE
 319
 320       IPZ2 = 4
 321
 322       DO 300  IS=1, NSEC
 323         IPZ1  = IPZ2
 324         IPZ2  = IPZ1 + 3
 325
 326         XT(3) = X(3) - 0.5 * (PAR(IPZ1) + PAR(IPZ2))
 327
 328         PT1 = 0.5 * (PAR(IPZ2) - PAR(IPZ1))
 329         IF (PT1 .LE. 0.0) GO TO 300
 330
 331         IF (PAR(IPZ1+1) .NE. PAR(IPZ2+1)) GO TO 250
 332         IF (PAR(IPZ1+2) .NE. PAR(IPZ2+2)) GO TO 250
 333
 334 C
 335 C......  This Z-section is a tube.
 336 C
 337         PT3 = PT1
 338         PT1 = PAR(IPZ1+1)
 339         PT2 = PAR(IPZ1+2)
 340         PT4 = PT6
 341         PT5 = PT7
 342
 343         CALL GNOTUB (XT, PT, 3, IFLAG, SNEXT, TSNXT, TSAFE)
 344         GO TO 280
 345
 346
 347   250   CONTINUE
 348 C
 349 C......   This Z-section is a cone.
 350 C
 351         PT2 = PAR(IPZ1+1)
 352         PT3 = PAR(IPZ1+2)
 353         PT4 = PAR(IPZ2+1)
 354         PT5 = PAR(IPZ2+2)
 355
 356         CALL GNOCON (XT, PT, 3, IFLAG, SNEXT, TSNXT, TSAFE)
 357
 358   280 CONTINUE
 359       IF (TSNXT .LT. SNXT)  SNXT = TSNXT
 360   300 CONTINUE
 361
 362   999 CONTINUE
 363       END