GEANT321/ggeom/gvdphi.F

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