GEANT321/ggeom/gvdphi.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/04 02/03/95  11.58.00  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GVDPHI(ISH,IROT,DX,PARS,CL,CH,IERR)
  13 C.
  14 C.    **********************************************************
  15 C.    *                                                        *
  16 C.    *    ROUTINE TO FIND THE PHI LIMITS OF THE OBJECT SHAPE  *
  17 C.    *    ISH ROTATED BY MATRIX IROT AND DISPLACED BY VECTOR  *
  18 C.    *    DX. IT HAS NPAR PARAMTERS IN THE ARRAY PARS. THE    *
  19 C.    *    LOWER LIMIT IS RETURNED IN CL AND THE HIGHER IN CH. *
  20 C.    *    NOTE THE OBJECT IS CONTAINED IN THE RANGE OF        *
  21 C.    *    INCREASING PHI FROM CL TO CH THOUGH CL AND CH ARE   *
  22 C.    *    FORCED TO LIE IN THE RANGE 0.0 TO 360.0 SO THAT THE *
  23 C.    *    VALUE OF CL CAN BE HIGHER THAN THAT OF CH.          *
  24 C.    *                                                        *
  25 C.    *    ==>Called by : GVDLIM modified                      *
  26 C.    *         Author  S.Giani  *********                     *
  27 C.    *                                                        *
  28 C.    **********************************************************
  29 C.
  30 #include "geant321/gcbank.inc"
  31 #include "geant321/gconsp.inc"
  32 #include "geant321/gcshno.inc"
  33       DIMENSION DX(3),PARS(50),X(3),XT(3)
  34 C.
  35 C.          -------------------------------------------
  36 C.
  37       IERR=1
  38 C
  39       DXS=DX(1)*DX(1)+DX(2)*DX(2)
  40       IF(DXS.GT.0.0) DXS=SQRT(DXS)
  41 C
  42       IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 40
  43 C
  44 C
  45 C           CUBOIDS, TRAPEZOIDS, PARALLELEPIPEDS.
  46 C
  47       IERR=0
  48       CL=0.0
  49       CH=360.0
  50 C
  51 C           IF IN DOUBT SET IT TO FULL RANGE.
  52 C
  53       IF(DXS.LE.0.0) GO TO 999
  54 C
  55       PHC=90.
  56       IF(DX(1).NE.0.)PHC=ATAN2(DX(2),DX(1))*RADDEG
  57       IF(DX(1).EQ.0..AND.DX(2).LT.0.)PHC=-90.
  58       IF(PHC.LT.0.0) PHC=PHC+360
  59       PL=0.0
  60       PH=0.0
  61 C
  62       DO 30 IP=1,8
  63 C
  64 C           THIS IS A LOOP OVER THE 8 CORNERS.
  65 C           FIRST FIND THE LOCAL COORDINATES.
  66 C
  67       IF(ISH.EQ.28) THEN
  68 C
  69 C            General twisted trapezoid.
  70 C
  71          IL=(IP+1)/2
  72          I0=IL*4+11
  73          IS=(IP-IL*2)*2+1
  74          X(3)=PARS(1)*IS
  75          X(1)=PARS(I0)+PARS(I0+2)*X(3)
  76          X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
  77          GO TO 20
  78 C
  79       ENDIF
  80 C
  81       IP3=ISH+2
  82       IF(ISH.EQ.10) IP3=3
  83       IF(ISH.EQ.4) IP3=1
  84       X(3)=PARS(IP3)
  85       IF(IP.LE.4) X(3)=-X(3)
  86       IP2=3
  87       IF(ISH.GT.2.AND.X(3).GT.0.0) IP2=4
  88       IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
  89       IF(ISH.EQ.4) IP2=4
  90       IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
  91       X(2)=PARS(IP2)
  92       IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
  93       IP1=1
  94       IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
  95       IF(ISH.EQ.4) IP1=5
  96       IF(ISH.EQ.4.AND.X(3).GT.0.0) IP1=IP1+4
  97       IF(ISH.EQ.4.AND.X(2).GT.0.0) IP1=IP1+1
  98       X(1)=PARS(IP1)
  99       IF(MOD(IP,2).EQ.1) X(1)=-X(1)
 100 C
 101       IF(ISH.NE.10) GO TO 10
 102       X(1)=X(1)+X(2)*PARS(4)+X(3)*PARS(5)
 103       X(2)=X(2)+X(3)*PARS(6)
 104    10 CONTINUE
 105 C
 106       IF(ISH.NE.4) GO TO 20
 107       IP4=7
 108       IF(X(3).GT.0.0) IP4=11
 109       X(1)=X(1)+X(2)*PARS(IP4)+X(3)*PARS(2)
 110       X(2)=X(2)+X(3)*PARS(3)
 111    20 CONTINUE
 112 C
 113 C          ROTATE.
 114 C
 115       JROT=LQ(JROTM-IROT)
 116       XT(1)=X(1)
 117       XT(2)=X(2)
 118       XT(3)=X(3)
 119       IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
 120 C
 121       XPT=DXS+(DX(1)*XT(1)+DX(2)*XT(2))/DXS
 122       YPT=(DX(1)*XT(2)-DX(2)*XT(1))/DXS
 123 C
 124       IF(YPT.EQ.0.0.AND.XPT.EQ.0.0) GO TO 999
 125       P=ATAN2(YPT,XPT)
 126       IF(P.GT.PI) P=P-PI*2.0
 127       IF(P.LT.PL) PL=P
 128       IF(P.GT.PH) PH=P
 129 C
 130 C
 131    30 CONTINUE
 132 C
 133 C
 134       IF(PH-PL.GT.PI) GO TO 999
 135       CL=PHC+PL*RADDEG
 136       CH=PHC+PH*RADDEG
 137 C
 138 ***      SG = SIGN(1.0,CL)
 139 ***      CL = MOD( ABS(CL),360.0 )
 140 ***      IF(SG.LE.0.0) CL=360.-CL
 141 ***      SG=SIGN(1.0,CH)
 142 ***      CH = MOD( ABS(CH),360.0 )
 143 ***      IF(SG.LE.0.0) CH=360.-CH
 144 C
 145       GO TO 999
 146 C
 147    40 CONTINUE
 148       MYFLAG=0
 149       IF(ISH.EQ.11)THEN
 150        NZLAST=PARS(4)
 151        IZLAST=2+3*NZLAST
 152        CLZ=PARS(5)
 153        CHZ=PARS(IZLAST)
 154        DZ2=ABS(CHZ-CLZ)
 155        DZ=DZ2*.5
 156        TMPRAD=0.
 157        TMPMIN=100000.
 158        DO 145 I=7,IZLAST+2,3
 159          IF(PARS(I).GT.TMPRAD)TMPRAD=PARS(I)
 160          IF(PARS(I-1).LT.TMPMIN)TMPMIN=PARS(I-1)
 161  145   CONTINUE
 162        RLOW=TMPMIN
 163        PHIMIN=PARS(1)
 164        PHIMAX=PHIMIN+PARS(2)
 165        AANG=ABS(PHIMAX-PHIMIN)
 166        NANG=PARS(3)
 167        AATMAX=NANG*360./AANG
 168        LATMAX=AATMAX
 169        ALA=AATMAX-LATMAX
 170        IF(ALA.GT..5)LATMAX=LATMAX+1
 171        AFINV=1./COS(PI/LATMAX)
 172        FINV=ABS(AFINV)
 173        RM=TMPRAD*FINV
 174        IF(PARS(2).EQ.360)THEN
 175         MYFLAG=5
 176        ELSE
 177         PHIMIN=PARS(1)*DEGRAD
 178         PHIMAX=(PARS(1)+PARS(2))*DEGRAD
 179         MYFLAG=6
 180        ENDIF
 181       ELSEIF(ISH.EQ.12)THEN
 182        NZLAST=PARS(3)
 183        IZLAST=1+3*NZLAST
 184        CLZ=PARS(4)
 185        CHZ=PARS(IZLAST)
 186        DZ2=ABS(CHZ-CLZ)
 187        DZ=DZ2*.5
 188        TMPRAD=0.
 189        TMPMIN=100000.
 190        DO 146 I=6,IZLAST+2,3
 191          IF(PARS(I).GT.TMPRAD)TMPRAD=PARS(I)
 192          IF(PARS(I-1).LT.TMPMIN)TMPMIN=PARS(I-1)
 193  146   CONTINUE
 194        RM=TMPRAD
 195        RLOW=TMPMIN
 196        IF(PARS(2).EQ.360)THEN
 197         MYFLAG=5
 198        ELSE
 199         PHIMIN=PARS(1)*DEGRAD
 200         PHIMAX=(PARS(1)+PARS(2))*DEGRAD
 201         MYFLAG=6
 202        ENDIF
 203       ELSEIF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13.AND.ISH.NE.14)THEN
 204         GO TO 80
 205       ENDIF
 206 C
 207 C             TUBES AND CONES.
 208 C
 209       IERR=0
 210       CL=0.0
 211       CH=360.0
 212 C
 213 C             WHEN IN DOUBT SET TO FULL RANGE.
 214 C
 215       IF(MYFLAG.EQ.0)RM=PARS(2)
 216       IF(ISH.LE.6.OR.ISH.EQ.NSCTUB.OR.MYFLAG.NE.0) GO TO 50
 217 **
 218       IF(ISH.EQ.13) THEN
 219 **
 220 **       approxime to a cylinder whit radius
 221 **       equal to the ellipse major axis
 222 **
 223          IF(PARS(1).GT.RM) RM=PARS(1)
 224          GOTO 50
 225       ENDIF
 226       IF(ISH.EQ.14) THEN
 227         RM = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
 228         GO TO 50
 229       ENDIF
 230       RM=PARS(3)
 231       IF(PARS(5).GT.RM) RM=PARS(5)
 232 C
 233    50 CONTINUE
 234 C
 235       IF(DXS.GT.RM) GO TO 70
 236       IF(ISH.EQ.5.OR.ISH.EQ.7.OR.ISH.EQ.14.OR.MYFLAG.EQ.5) GO TO 999
 237       IF(ISH.EQ.13) GOTO 999
 238 *                 Here we treat the CONS and TUBS
 239 *                 This is the simple case, no rotation
 240 *                 Compute the position of the limits on
 241 *                 the X-Y plane.
 242       IF(MYFLAG.EQ.0.AND.ISH.LE.6)THEN
 243        PHIMIN=PARS(4)*DEGRAD
 244        PHIMAX=PARS(5)*DEGRAD
 245        RLOW=PARS(1)
 246       ELSEIF(MYFLAG.EQ.0.AND.ISH.EQ.8)THEN
 247         PHIMIN=PARS(6)*DEGRAD
 248         PHIMAX=PARS(7)*DEGRAD
 249         RLOW=MIN(PARS(2),PARS(4))
 250       ENDIF
 251       CL=PHIMIN*RADDEG
 252       CH=PHIMAX*RADDEG
 253       IF(DXS.NE.0.)THEN
 254        DDX1 = DX(1)+RM*COS(PHIMIN)
 255        DDY1 = DX(2)+RM*SIN(PHIMIN)
 256        DDX2 = DX(1)+RM*COS(PHIMAX)
 257        DDY2 = DX(2)+RM*SIN(PHIMAX)
 258        CLA = ATAN2(DDY1,DDX1)*RADDEG
 259        CHA = ATAN2(DDY2,DDX2)*RADDEG
 260        DDX1 = DX(1)+RLOW*COS(PHIMIN)
 261        DDY1 = DX(2)+RLOW*SIN(PHIMIN)
 262        DDX2 = DX(1)+RLOW*COS(PHIMAX)
 263        DDY2 = DX(2)+RLOW*SIN(PHIMAX)
 264        CLB = ATAN2(DDY1,DDX1)*RADDEG
 265        CHB = ATAN2(DDY2,DDX2)*RADDEG
 266        CL=MIN(CLA,CLB,CHA,CHB)
 267        CH=MAX(CLA,CLB,CHA,CHB)
 268        IF((CH-CL).GT.(PHIMAX-PHIMIN)*RADDEG)THEN
 269          IF(ISH.EQ.6.OR.ISH.EQ.8.OR.MYFLAG.EQ.6)THEN
 270            IF(DXS.GT.RLOW)THEN
 271               CL=0.
 272               CH=360.
 273            ENDIF
 274          ELSE
 275            CL=0.
 276            CH=360.
 277          ENDIF
 278        ENDIF
 279       ENDIF
 280 C
 281    60 CONTINUE
 282 C
 283       IF(IROT.EQ.0) GO TO 65
 284       IF(CL.EQ.0..AND.CH.EQ.360.)GOTO 65
 285       JROT=LQ(JROTM-IROT)
 286       IF(Q(JROT+15).NE.0.0.AND.Q(JROT+15).NE.180.0)THEN
 287         CL=0.
 288         CH=360.
 289         GO TO 999
 290       ENDIF
 291 C
 292       PHX=Q(JROT+12)
 293       PHY=Q(JROT+14)
 294       IF(PHY.LT.PHX) PHY=PHY+360.0
 295       ISPH=1
 296       IF(PHY-PHX.GT.180.0) ISPH=-1
 297       IF(DXS.NE.0.)THEN
 298        PHI1 = ISPH*PHIMIN+PHX*DEGRAD
 299        PHI2 = ISPH*PHIMAX+PHX*DEGRAD
 300        DDX1 = DX(1)+RM*COS(PHI1)
 301        DDY1 = DX(2)+RM*SIN(PHI1)
 302        DDX2 = DX(1)+RM*COS(PHI2)
 303        DDY2 = DX(2)+RM*SIN(PHI2)
 304        CLA = ATAN2(DDY1,DDX1)*RADDEG
 305        CHA = ATAN2(DDY2,DDX2)*RADDEG
 306        DDX1 = DX(1)+RLOW*COS(PHI1)
 307        DDY1 = DX(2)+RLOW*SIN(PHI1)
 308        DDX2 = DX(1)+RLOW*COS(PHI2)
 309        DDY2 = DX(2)+RLOW*SIN(PHI2)
 310        CLB = ATAN2(DDY1,DDX1)*RADDEG
 311        CHB = ATAN2(DDY2,DDX2)*RADDEG
 312        CL=MIN(CLA,CLB,CHA,CHB)
 313        CH=MAX(CLA,CLB,CHA,CHB)
 314       ELSE
 315        CL=ISPH*CL+PHX
 316        CH=ISPH*CH+PHX
 317       ENDIF
 318       IF(ISPH.EQ.1) GO TO 65
 319       CHT=CH
 320       CH=CL
 321       CL=CHT
 322 C
 323    65 CONTINUE
 324 C
 325 ***      SG=SIGN(1.0,CL)
 326 ***      CL = MOD( ABS(CL),360.0 )
 327 ***      IF(SG.LE.0.0) CL=360.-CL
 328 ***      SG=SIGN(1.0,CH)
 329 ***      CH = MOD( ABS(CH),360.0 )
 330 ***      IF(SG.LE.0.0) CH=360.-CH
 331 C
 332       GO TO 999
 333 C
 334    70 CONTINUE
 335 C
 336 C            DISPLACEMENT GREATER THAN MAXIMUM RADIUS SO
 337 C            ASSUME COMPLETE TUBE TO GENERATE 'WORST CASE'.
 338 C
 339       IF(MYFLAG.EQ.0)DZ=PARS(3)
 340       IF(ISH.EQ.NSCTUB) THEN
 341         S1 = (1.0-PARS(8))*(1.0+PARS(8))
 342         IF( S1 .GT. 0.0) S1 = SQRT(S1)
 343         S2 = (1.0-PARS(11))*(1.0+PARS(11))
 344         IF( S2 .GT. 0.0) S2 = SQRT(S2)
 345         IF( S2 .GT. S1 ) S1 = S2
 346         DZ = DZ+RM*S1
 347       ELSEIF(ISH.GT.6.AND.ISH.NE.13.AND.ISH.NE.14.AND.MYFLAG.EQ.0)THEN
 348         DZ=PARS(1)
 349       ENDIF
 350 C
 351       X(1)=0.0
 352       X(2)=0.0
 353       X(3)=1.0
 354 C
 355 C                    LOCAL Z AXIS.
 356 C
 357       JROT=LQ(JROTM-IROT)
 358       XT(1)=X(1)
 359       XT(2)=X(2)
 360       XT(3)=X(3)
 361       IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
 362 C
 363       COST=ABS(DX(1)*XT(1)+DX(2)*XT(2))
 364       COS2=ABS(DX(1)*XT(2)-DX(2)*XT(1))
 365       SINT=(DXS+COST)*(DXS-COST)
 366       SIN2=(DXS+COS2)*(DXS-COS2)
 367       IF(SINT.GT.0.0) SINT=SQRT(SINT)
 368       IF(SIN2.GT.0.0) SIN2=SQRT(SIN2)
 369 C
 370       XPT=DXS-(COST*DZ+SINT*RM)/DXS
 371 C
 372       IF(XPT.LE.0.0) GO TO 999
 373       YPT=(SIN2*RM+COS2*DZ)/DXS
 374       DP=ATAN(YPT/XPT)
 375 C
 376       P0=ATAN2(DX(2),DX(1))
 377       CL=(P0-DP)*RADDEG
 378       CH=(P0+DP)*RADDEG
 379 C
 380 ***      SG=SIGN(1.0,CL)
 381 ***      CL = MOD( ABS(CL),360.0 )
 382 ***      IF(SG.LE.0.0) CL=360.-CL
 383 ***      SG=SIGN(1.0,CH)
 384 ***      CH = MOD( ABS(CH),360.0 )
 385 ***      IF(SG.LE.0.0) CH=360.-CH
 386 C
 387       GO TO 999
 388 C
 389    80 CONTINUE
 390       IF(ISH.GT.9.AND.MYFLAG.EQ.0) GO TO 999
 391 C
 392 C               SPHERE.
 393 C
 394       IERR=0
 395       CL=0.0
 396       CH=360.0
 397 C
 398       IF(IROT.NE.0.OR.DXS.GT.0.0) GO TO 90
 399 C
 400 C          UNROTATED AND CENTERED.
 401 C
 402       CL=PARS(5)
 403       CH=PARS(6)
 404 C
 405 ***      SG=SIGN(1.0,CL)
 406 ***      CL = MOD( ABS(CL),360.0 )
 407 ***      IF(SG.LE.0.0) CL=360.-CL
 408 ***      SG=SIGN(1.0,CH)
 409 ***      CH = MOD( ABS(CH),360.0 )
 410 ***      IF(SG.LE.0.0) CH=360.-CH
 411 C
 412       GO TO 999
 413 C
 414    90 CONTINUE
 415 C
 416 C            ROTATED OR NOT CENTERED.
 417 C
 418       IF(DXS.LT.PARS(2)) GO TO 999
 419       P0=ATAN2(DX(2),DX(1))
 420       DP=ASIN(PARS(2)/DXS)
 421       CL=(P0-DP)*RADDEG
 422       CH=(P0+DP)*RADDEG
 423 C
 424 ***      SG=SIGN(1.0,CL)
 425 ***      CL = MOD( ABS(CL),360.0 )
 426 ***      IF(SG.LE.0.0) CL=360.-CL
 427 ***      SG=SIGN(1.0,CH)
 428 ***      CH = MOD( ABS(CH),360.0 )
 429 ***      IF(SG.LE.0.0) CH=360.-CH
 430 C
 431   999 CONTINUE
 432       END