GEANT321/ggeom/gvdthe.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:57  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.31  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GVDTHE(ISH,IROT,DX,PARS,CL,CH,IERR)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *    ROUTINE TO COMPUTE THE THETA LIMITS FOR VOLUME OF SHAPE     *
  17 C.    *    ISH ROTATED BY MATRIX IROT AND DISPLACED BY VECTOR DX.      *
  18 C.    *    THE VOLUME HAS NPAR PARAMETERS IN THE ARRAY PARS. THE LOWER *
  19 C.    *    LIMIT IS RETURNED IN CL THE HIGHER IN CH. IF THE            *
  20 C.    *    CALCULATION CANNOT BE MADE IERR IS SET TO 1 OTHERWISE IT    *
  21 C.    *    IS SET TO 0.                                                *
  22 C.    *                                                                *
  23 C.    *    ==>Called by : GVDLIM                                       *
  24 C.    *         Author  S.Giani  *********                             *
  25 C.    *                                                                *
  26 C.    ******************************************************************
  27 C.
  28 #include "geant321/gcbank.inc"
  29 #include "geant321/gconsp.inc"
  30 #include "geant321/gcshno.inc"
  31       DIMENSION DX(3),PARS(50),X(3),XT(3),X1(3),X2(3),XT1(3),XT2(3)
  32 C.
  33 C.          ----------------------------------------------
  34 C.
  35       IERR=1
  36 C
  37       DXS=DX(1)*DX(1)+DX(2)*DX(2)
  38       IF(DXS.GT.0.0) DXS=SQRT(DXS)
  39 C
  40       IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 130
  41 C
  42 C            CUBES, TRAPEZOIDS AND PARALLELEPIPEDS.
  43 C
  44       IERR=0
  45       CL=0.0
  46       CH=180.0
  47 C
  48       IF(DXS.LE.0.0) GO TO 999
  49 C
  50       TH=90.
  51       IF(ABS(DX(3)).LT.1.0E-06)GO TO 5
  52       TH=ATAN(DXS/DX(3))*RADDEG
  53       IF(TH.LT.0.0) TH=TH+180.0
  54    5  TL=TH
  55 C
  56       DO 50 IP=1,8
  57 C
  58 C           THIS IS A LOOP OVER THE 8 CORNERS.
  59 C           FIRST FIND THE LOCAL COORDINATES.
  60 C
  61       IF(ISH.EQ.28) THEN
  62          goto 999
  63 C            General twisted trapezoid.
  64 C
  65 *         IL=(IP+1)/2
  66 *         I0=IL*4+11
  67 *         IS=(IP-IL*2)*2+1
  68 *         X(3)=PARS(1)*IS
  69 *         X(1)=PARS(I0)+PARS(I0+2)*X(3)
  70 *         X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
  71 *         GO TO 20
  72 C
  73       ENDIF
  74 C
  75       IP3=ISH+2
  76       IF(ISH.EQ.10) IP3=3
  77       IF(ISH.EQ.4) IP3=1
  78       X(3)=PARS(IP3)
  79       IF(IP.LE.4) X(3)=-X(3)
  80       IP2=3
  81       IF(ISH.GT.2.AND.X(3).GT.0.0) IP2=4
  82       IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
  83       IF(ISH.EQ.4) IP2=4
  84       IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
  85       X(2)=PARS(IP2)
  86       IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
  87       IP1=1
  88       IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
  89       IF(ISH.EQ.4) IP1=5
  90       IF(ISH.EQ.4.AND.X(3).GT.0.0) IP1=IP1+4
  91       IF(ISH.EQ.4.AND.X(2).GT.0.0) IP1=IP1+1
  92       X(1)=PARS(IP1)
  93       IF(MOD(IP,2).EQ.1) X(1)=-X(1)
  94 C
  95       IF(ISH.NE.10) GO TO 10
  96       X(1)=X(1)+X(2)*PARS(4)+X(3)*PARS(5)
  97       X(2)=X(2)+X(3)*PARS(6)
  98    10 CONTINUE
  99 C
 100       IF(ISH.NE.4) GO TO 20
 101       IP4=7
 102       IF(X(3).GT.0.0) IP4=11
 103       X(1)=X(1)+X(2)*PARS(IP4)+X(3)*PARS(2)
 104       X(2)=X(2)+X(3)*PARS(3)
 105    20 CONTINUE
 106 C
 107 C          ROTATE.
 108 C
 109       JROT=LQ(JROTM-IROT)
 110       XT(1)=X(1)
 111       XT(2)=X(2)
 112       XT(3)=X(3)
 113       IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
 114 C
 115       Z=DX(3)+XT(3)
 116       R=(DX(1)+XT(1))**2+(DX(2)+XT(2))**2
 117       RT=R+(DX(3)+XT(3))**2
 118       IF(RT.LT.1.0E-10) GO TO 999
 119 C
 120       IF(R.GT.0.0) GO TO 30
 121       IF(Z.GT.0.0) T=0.0
 122       IF(Z.LT.0.0) T=180.0
 123       GO TO 40
 124 C
 125    30 CONTINUE
 126       T=90.0
 127       IF(ABS(Z).LT.1.0E-6) GO TO 40
 128 C
 129       R=SQRT(R)
 130       T=ATAN(R/Z)*RADDEG
 131       IF(T.EQ.0.0.AND.Z.LT.0.0) T=180.0
 132       IF(T.LT.0.0) T=T+180.0
 133 C
 134    40 CONTINUE
 135       IF(T.GT.TH) TH=T
 136       IF(T.LT.TL) TL=T
 137 C
 138    50 CONTINUE
 139 C
 140 C           THETA LIMITS SET FROM THE POINTS NOW DO THE EDGES.
 141 C
 142       DO 120 IL=1,12
 143 C
 144 C           FIND THE END POINT NUMBERS FOR EACH EDGE.
 145 C
 146       IF(IL.GT.4) GO TO 60
 147       IPP1=1
 148       IF(IL.GT.2) IPP1=4
 149       IPP2=2
 150       IF(MOD(IL,2).EQ.1) IPP2=3
 151 C
 152       GO TO 80
 153    60 CONTINUE
 154       IF(IL.LT.9) GO TO 70
 155 C
 156       IPP1=5
 157       IF(IL.GT.10) IPP1=8
 158       IPP2=6
 159       IF(MOD(IL,2).EQ.1) IPP2=7
 160 C
 161       GO TO 80
 162    70 CONTINUE
 163 C
 164       IPP1=IL-4
 165       IPP2=IL
 166 C
 167    80 CONTINUE
 168 C
 169 C           NOW GET THE POINTS AND ROTATE THEM.
 170 C
 171       IF(ISH.EQ.28) THEN
 172 C
 173 C            General twisted trapezoid.
 174 C
 175          ILL=IPP1
 176          IF(IPP1.EQ.3) ILL=4
 177          IF(IPP1.EQ.4) ILL=3
 178          I0=ILL*4+11
 179          X1(3)=PARS(1)
 180          IF(IPP1.LT.5) X1(3)=-X1(3)
 181          X1(1)=PARS(I0)+PARS(I0+2)*X1(3)
 182          X1(2)=PARS(I0+1)+PARS(I0+3)*X1(3)
 183          ILL=IPP2
 184          IF(IPP2.EQ.3) ILL=4
 185          IF(IPP2.EQ.4) ILL=3
 186          I0=ILL*4+11
 187          X2(3)=PARS(1)
 188          IF(IPP2.LT.5) X2(3)=-X2(3)
 189          X2(1)=PARS(I0)+PARS(I0+2)*X2(3)
 190          X2(2)=PARS(I0+1)+PARS(I0+3)*X2(3)
 191 C
 192          GO TO 100
 193 C
 194       ENDIF
 195 C
 196       IP3=ISH+2
 197       IF(ISH.EQ.10) IP3=3
 198       IF(ISH.EQ.4) IP3=1
 199       X1(3)=PARS(IP3)
 200       IF(IPP1.LE.4) X1(3)=-X1(3)
 201       IP2=3
 202       IF(ISH.GT.2.AND.X1(3).GT.0.0) IP2=4
 203       IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
 204       IF(ISH.EQ.4) IP2=4
 205       IF(ISH.EQ.4.AND.X1(3).GT.0.0) IP2=8
 206       X1(2)=PARS(IP2)
 207       IF(MOD(IPP1+3,4).LT.2) X1(2)=-X1(2)
 208       IP1=1
 209       IF(ISH.NE.1.AND.ISH.NE.10.AND.X1(3).GT.0.0) IP1=2
 210       IF(ISH.EQ.4) IP1=5
 211       IF(ISH.EQ.4.AND.X1(3).GT.0.0) IP1=IP1+4
 212       IF(ISH.EQ.4.AND.X1(2).GT.0.0) IP1=IP1+1
 213       X1(1)=PARS(IP1)
 214       IF(MOD(IPP1,2).EQ.1) X1(1)=-X1(1)
 215 C
 216       IP3=ISH+2
 217       IF(ISH.EQ.10) IP3=3
 218       IF(ISH.EQ.4) IP3=1
 219       X2(3)=PARS(IP3)
 220       IF(IPP2.LE.4) X2(3)=-X2(3)
 221       IP2=3
 222       IF(ISH.GT.2.AND.X2(3).GT.0.0) IP2=4
 223       IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
 224       IF(ISH.EQ.4) IP2=4
 225       IF(ISH.EQ.4.AND.X2(3).GT.0.0) IP2=8
 226       X2(2)=PARS(IP2)
 227       IF(MOD(IPP2+3,4).LT.2) X2(2)=-X2(2)
 228       IP1=1
 229       IF(ISH.NE.1.AND.ISH.NE.10.AND.X2(3).GT.0.0) IP1=2
 230       IF(ISH.EQ.4) IP1=5
 231       IF(ISH.EQ.4.AND.X2(3).GT.0.0) IP1=IP1+4
 232       IF(ISH.EQ.4.AND.X2(2).GT.0.0) IP1=IP1+1
 233       X2(1)=PARS(IP1)
 234       IF(MOD(IPP2,2).EQ.1) X2(1)=-X2(1)
 235 C
 236       IF(ISH.NE.10) GO TO 90
 237       X1(1)=X1(1)+X1(2)*PARS(4)+X1(3)*PARS(5)
 238       X1(2)=X1(2)+X1(3)*PARS(6)
 239       X2(1)=X2(1)+X2(2)*PARS(4)+X2(3)*PARS(5)
 240       X2(2)=X2(2)+X2(3)*PARS(6)
 241    90 CONTINUE
 242 C
 243       IF(ISH.NE.4) GO TO 100
 244       IP4=7
 245       IF(X1(3).GT.0.0) IP4=11
 246       X1(1)=X1(1)+X1(2)*PARS(IP4)+X1(3)*PARS(2)
 247       X1(2)=X1(2)+X1(3)*PARS(3)
 248       IP4=7
 249       IF(X2(3).GT.0.0) IP4=11
 250       X2(1)=X2(1)+X2(2)*PARS(IP4)+X2(3)*PARS(2)
 251       X2(2)=X2(2)+X2(3)*PARS(3)
 252   100 CONTINUE
 253 C
 254 C          ROTATE.
 255 C
 256       JROT=LQ(JROTM-IROT)
 257       XT1(1)=X1(1)
 258       XT1(2)=X1(2)
 259       XT1(3)=X1(3)
 260       IF(IROT.NE.0) CALL GINROT(X1,Q(JROT+1),XT1)
 261       XT2(1)=X2(1)
 262       XT2(2)=X2(2)
 263       XT2(3)=X2(3)
 264       IF(IROT.NE.0) CALL GINROT(X2,Q(JROT+1),XT2)
 265 C
 266 C           NOW WE HAVE THE END POINTS IN THE MASTER SYSTEM.
 267 C           FIND THE TANGENT POINT TO THE SET OF CONES OF CONSTANT
 268 C           THETA.
 269 C
 270       DS=(XT2(1)-XT1(1))**2+(XT2(2)-XT1(2))**2+(XT2(3)-XT1(3))**2
 271       IF(DS.LE.0.0) GO TO 120
 272 C
 273       DS=SQRT(DS)
 274 C
 275       X0=(XT2(1)+XT1(1))/2.0+DX(1)
 276       Y0=(XT2(2)+XT1(2))/2.0+DX(2)
 277       Z0=(XT2(3)+XT1(3))/2.0+DX(3)
 278       ALX=(XT2(1)-XT1(1))/DS
 279       ALY=(XT2(2)-XT1(2))/DS
 280       ALZ=(XT2(3)-XT1(3))/DS
 281 C
 282       T=90.0
 283       IF(Z0.EQ.0.0.AND.ALZ.EQ.0.0) GO TO 110
 284 C
 285       IF(ALX.EQ.0.0.AND.ALY.EQ.0.0) GO TO 120
 286 C             THIS CHECKS WHETHER THE LINE IS PARALLEL TO THE
 287 C             Z AXIS IN WHICH CASE THERE IS NO TANGENT AND
 288 C             THE END POINTS DETERMINE THE THETA RANGE.
 289 C
 290       SNUM=(X0*Z0*ALX+Y0*Z0*ALY+X0*X0*ALZ+Y0*Y0*ALZ)
 291       SDEN=(Z0*ALX*ALX-X0*ALX*ALZ+Z0*ALY*ALY-Y0*ALY*ALZ)
 292 C
 293       IF(ABS(SNUM).GT.0.5*DS*ABS(SDEN)) GO TO 120
 294 C
 295 C           TANGENT EXIST BETWEEN THE TWO ENDS.
 296 C
 297       S = SNUM/SDEN
 298       X0=X0+S*ALX
 299       Y0=Y0+S*ALY
 300       Z0=Z0+S*ALZ
 301 C
 302       R=X0*X0+Y0*Y0
 303       RT=R+Z0*Z0
 304 C
 305       IF(RT.LT.1.0E-10) GO TO 999
 306 C
 307       IF(R.GT.0.0) R=SQRT(R)
 308       T=90.0
 309       IF(ABS(Z0).LT.1.0E-06) GO TO 110
 310 C
 311       T=ATAN(R/Z0)*RADDEG
 312       IF(T.EQ.0.0.AND.Z0.LT.0.0) T=180.0
 313       IF(T.LT.0.0) T=T+180.0
 314 C
 315   110 CONTINUE
 316 C
 317       IF(T.LT.TL) TL=T
 318       IF(T.GT.TH) TH=T
 319 C
 320 C        CHECK FOR THE POSSIBILITY OF STRADDLING T=0.0 &/OR 180.0.
 321 C
 322       C=X0*DX(1)+Y0*DX(2)
 323       IF(C.GT.0.0) GO TO 120
 324 C
 325 C          CHECK IF SAME SIGN OF Z.
 326 C
 327       IF(Z0*DX(3).LT.0.0) GO TO 999
 328 C
 329       T=0.0
 330       IF(Z0.LT.0.0) T=180.0
 331       IF(T.LT.TL) TL=T
 332       IF(T.GT.TH) TH=T
 333 C
 334   120 CONTINUE
 335 C
 336 C            DONE SET THE LIMITS.
 337 C
 338       CL=TL
 339       CH=TH
 340 C
 341       GO TO 999
 342 C
 343   130 CONTINUE
 344 C
 345 C          TUBES, CONES ETC.
 346 C
 347       IF(IROT.NE.0.OR.DXS.GT.1.0E-05) GO TO 180
 348 C
 349 C          UNROTATED AND CENTERED.
 350 C
 351       IF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13
 352      +           .AND.ISH.NE.14)GO TO 170
 353 C
 354 C               TUBES AND CONES.
 355 C
 356       IERR=0
 357       DZ=PARS(3)
 358       RMN=PARS(1)
 359       RMX=PARS(2)
 360       IF(ISH.EQ.13) THEN
 361 **
 362 **         approxime to a cylinder whit radius
 363 **         equal to the ellipse major axis
 364 **
 365         RMN=0.0
 366         IF(PARS(1).GT.RMX) RMX=PARS(1)
 367         GOTO 140
 368       ENDIF
 369       IF(ISH.EQ.14) THEN
 370 C   not really sure of the function of these... keep RMN=PARS(1)
 371         RMX = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
 372         GO TO 140
 373       ENDIF
 374       IF(ISH.EQ.NSCTUB) THEN
 375         S1 = (1.0-PARS(8))*(1.0+PARS(8))
 376         IF( S1 .GT. 0.0) S1 = SQRT(S1)
 377         S2 = (1.0-PARS(11))*(1.0+PARS(11))
 378         IF( S2 .GT. 0.0) S2 = SQRT(S2)
 379         IF( S2 .GT. S1 ) S1 = S2
 380         DZ = DZ+RMX*S1
 381       ENDIF
 382       IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 140
 383       DZ=PARS(1)
 384       RMN=PARS(2)
 385       RMX=PARS(3)
 386       IF(PARS(4).LT.RMN) RMN=PARS(4)
 387       IF(PARS(5).GT.RMX) RMX=PARS(5)
 388   140 CONTINUE
 389 C
 390       IF(DZ.GT.ABS(DX(3))) GO TO 160
 391 C
 392 C          ALL FORWARD OR ALL BACK.
 393 C
 394       DSH=DX(3)-DZ
 395       DLN=DX(3)+DZ
 396       IF(DX(3).GT.0.0) GO TO 150
 397       DSS=DSH
 398       DSH=DLN
 399       DLN=DSS
 400   150 CONTINUE
 401 C
 402       CL=90.0
 403       CH=90.0
 404       IF(DLN.NE.0.0) CL=ATAN(RMN/DLN)*RADDEG
 405       IF(DSH.NE.0.0) CH=ATAN(RMX/DSH)*RADDEG
 406       IF(DX(3).GT.0.0) GO TO 999
 407       CS=CL
 408       CL=CH
 409       CH=CS
 410       IF(CH.EQ.0.0) CH=180.0
 411       IF(CL.LT.0.0) CL=CL+180.0
 412       IF(CH.LT.0.0) CH=CH+180.0
 413 C
 414       GO TO 999
 415 C
 416   160 CONTINUE
 417 C
 418       CL=90.0
 419       CH=90.0
 420       IF(DZ+DX(3).NE.0.0) CL=ATAN(RMN/(DZ+DX(3)))*RADDEG
 421       IF(-DZ+DX(3).NE.0.0) CH=ATAN(RMN/(-DZ+DX(3)))*RADDEG
 422       IF(CH.LE.0.0) CH=CH+180.0
 423 C
 424       GO TO 999
 425 C
 426   170 CONTINUE
 427       IF(ISH.GT.9) GO TO 999
 428 C
 429 C           SPHERE.
 430 C
 431       IERR=0
 432       CL=PARS(3)
 433       CH=PARS(4)
 434 C
 435       GO TO 999
 436   180 CONTINUE
 437 C
 438       IF(ISH.EQ.11.OR.ISH.EQ.12) GOTO 999
 439 **
 440       RM=PARS(2)
 441       IF(ISH.EQ.9) GO TO 200
 442 C
 443       DZ=PARS(3)
 444       IF(ISH.EQ.NSCTUB) THEN
 445         S1 = (1.0-PARS(8))*(1.0+PARS(8))
 446         IF( S1 .GT. 0.0) S1 = SQRT(S1)
 447         S2 = (1.0-PARS(11))*(1.0+PARS(11))
 448         IF( S2 .GT. 0.0) S2 = SQRT(S2)
 449         IF( S2 .GT. S1 ) S1 = S2
 450         DZ = DZ+RM*S1
 451       ENDIF
 452 **
 453       IF(ISH.EQ.13) THEN
 454         IF(PARS(1).GT.RM) RM=PARS(1)
 455         GOTO 190
 456       ENDIF
 457 **
 458       IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 190
 459 C
 460       DZ=PARS(1)
 461       RM=PARS(3)
 462       IF(PARS(5).GT.RM) RM=PARS(5)
 463 C
 464   190 CONTINUE
 465 C
 466       RM=SQRT(RM**2+DZ**2)
 467 C
 468   200 CONTINUE
 469 C
 470       CL=0.0
 471       CH=180.0
 472       IERR=0
 473       RC=DXS**2+DX(3)**2
 474       IF(RC.GT.0.0) RC=SQRT(RC)
 475       IF(RM.GE.RC) GO TO 999
 476 C
 477       TC=90.0
 478       IF(ABS(DX(3)).GT.0.0) TC=ATAN(DXS/DX(3))*RADDEG
 479       IF(TC.LT.0.0) TC=TC+180.0
 480 C
 481       DT=ASIN(RM/RC)*RADDEG
 482       CL=TC-DT
 483       CH=TC+DT
 484       IF(CL.LT.0.0) CL=0.0
 485       IF(CH.GT.180.0) CH=180.0
 486 C
 487   999 CONTINUE
 488       END