| fe4da5cc |
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 |
git://git.uio.no / u / mrichter / AliRoot.git / blame