GEANT321/ggeom/gvdcar.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:56  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/03 10/10/94  20.01.58  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GVDCAR(IAXIS,ISH,IROT,PARS,CL,CH,IERR)
  13 C.
  14 C.    *****************************************************************
  15 C.    *                                                               *
  16 C.    *    ROUTINE TO FIND THE LIMITS ALONG AXIS IAXIS IN CARTESIAN   *
  17 C.    *    COORDINATES FOR VOLUME OF SHAPE ISH ROTATED BY THE         *
  18 C.    *    ROTATION MATRIX IROT. THE SHAPE HAS NPAR PARAMETERS IN     *
  19 C.    *    THE ARRAY PARS. THE LOWER LIMIT IS RETURNED IN CL, THE     *
  20 C.    *    HIGHER IN CH. IF THE CALCULATION CANNOT BE MADE IERR IS    *
  21 C.    *    SET TO 1 OTHERWISE IT 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 PARS(50),X(3),XT(3)
  32 C.
  33 C.          ---------------------------------------------------
  34 C.
  35       IERR=1
  36       IF (ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 40
  37 C
  38 C           CUBOIDS, TRAPEZOIDS, PARALLELEPIPEDS.
  39 C
  40 C
  41       IERR=0
  42       CL=0
  43       CH=0
  44 C
  45       DO 30 IP=1,8
  46 C
  47 C           THIS IS A LOOP OVER THE 8 CORNERS.
  48 C           FIRST FIND THE LOCAL COORDINATES.
  49 C
  50       IF(ISH.EQ.28) THEN
  51 C
  52 C            General twisted trapezoid.
  53 C
  54          IL=(IP+1)/2
  55          I0=IL*4+11
  56          IS=(IP-IL*2)*2+1
  57          X(3)=PARS(1)*IS
  58          X(1)=PARS(I0)+PARS(I0+2)*X(3)
  59          X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
  60          GO TO 20
  61 C
  62       ENDIF
  63 C
  64       IP3=ISH+2
  65       IF(ISH.EQ.10) IP3=3
  66       IF(ISH.EQ.4) IP3=1
  67       X(3)=PARS(IP3)
  68       IF(IP.LE.4) X(3)=-X(3)
  69       IP2=3
  70       IF(ISH.GT.2.AND.X(3).GT.0.0) IP2=4
  71       IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
  72       IF(ISH.EQ.4) IP2=4
  73       IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
  74       X(2)=PARS(IP2)
  75       IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
  76       IP1=1
  77       IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
  78       IF(ISH.EQ.4) IP1=5
  79       IF(ISH.EQ.4.AND.X(3).GT.0.0) IP1=IP1+4
  80       IF(ISH.EQ.4.AND.X(2).GT.0.0) IP1=IP1+1
  81       X(1)=PARS(IP1)
  82       IF(MOD(IP,2).EQ.1) X(1)=-X(1)
  83 C
  84       IF(ISH.NE.10) GO TO 10
  85       X(1)=X(1)+X(2)*PARS(4)+X(3)*PARS(5)
  86       X(2)=X(2)+X(3)*PARS(6)
  87    10 CONTINUE
  88 C
  89       IF(ISH.NE.4) GO TO 20
  90       IP4=7
  91       IF(X(3).GT.0.0) IP4=11
  92       X(1)=X(1)+X(2)*PARS(IP4)+X(3)*PARS(2)
  93       X(2)=X(2)+X(3)*PARS(3)
  94    20 CONTINUE
  95 C
  96 C          ROTATE.
  97 C
  98       JROT=LQ(JROTM-IROT)
  99       XT(1)=X(1)
 100       XT(2)=X(2)
 101       XT(3)=X(3)
 102       IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
 103 C
 104 C          UPDATE LIMITS IF NECESSARY.
 105 C
 106       IF(XT(IAXIS).LT.CL) CL=XT(IAXIS)
 107       IF(XT(IAXIS).GT.CH) CH=XT(IAXIS)
 108 C
 109    30 CONTINUE
 110 C
 111       GO TO 999
 112 C
 113    40 CONTINUE
 114       IF(ISH.EQ.9) GO TO 90
 115 C
 116 C              TUBES , CONES, POLYGONS, POLYCONES.
 117 C              AND CUT TUBES.
 118 C
 119       MYFLAG=0
 120       IF((ISH.EQ.11.OR.ISH.EQ.12).AND.(IAXIS.LT.3))THEN
 121         MYFLAG=1
 122       ENDIF
 123       X(1)=0.0
 124       X(2)=0.0
 125       X(3)=1.0
 126       JROT=LQ(JROTM-IROT)
 127       XT(1)=X(1)
 128       XT(2)=X(2)
 129       XT(3)=X(3)
 130       IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
 131 C
 132 C          XT IS Z AXIS ROTATED.
 133 C
 134       IF(MYFLAG.EQ.0)THEN
 135        IF(ABS(XT(IAXIS)).LT.0.99) GO TO 50
 136       ELSE
 137        IF(ABS(XT(3)).LT.0.99) GO TO 50
 138       ENDIF
 139       IF(ISH.EQ.11)GO TO 45
 140       IF(ISH.EQ.12)GO TO 46
 141 C
 142 C           PARALLEL.
 143 C
 144       IP=3
 145       IF(ISH.GT.6.AND.ISH.NE.NSCTUB.AND.ISH.NE.13.AND.ISH.NE.14) IP=1
 146       CL=-PARS(IP)
 147       CH=PARS(IP)
 148       IERR=0
 149 C
 150       GO TO 999
 151   45  IF(MYFLAG.EQ.0)THEN
 152        NZLAST=PARS(4)
 153        IZLAST=2+3*NZLAST
 154        CL=PARS(5)
 155        GO TO 49
 156       ELSEIF(MYFLAG.EQ.1)THEN
 157        NZLAST=PARS(4)
 158        IZLAST=2+3*NZLAST
 159        TMPRAD=0.
 160        DO 145 I=7,IZLAST+2,3
 161          IF(PARS(I).GT.TMPRAD)TMPRAD=PARS(I)
 162  145   CONTINUE
 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        R=TMPRAD*FINV
 174        CL=-R
 175        CH= R
 176        IERR=0
 177        GOTO 999
 178       ENDIF
 179 C
 180   46  IF(MYFLAG.EQ.0)THEN
 181        NZLAST=PARS(3)
 182        IZLAST=1+3*NZLAST
 183        CL=PARS(4)
 184       ELSEIF(MYFLAG.EQ.1)THEN
 185        NZLAST=PARS(3)
 186        IZLAST=1+3*NZLAST
 187        TMPRAD=0.
 188        DO 146 I=6,IZLAST+2,3
 189          IF(PARS(I).GT.TMPRAD)TMPRAD=PARS(I)
 190  146   CONTINUE
 191        CL=-TMPRAD
 192        CH= TMPRAD
 193        IERR=0
 194        GOTO 999
 195       ENDIF
 196 C
 197   49  CH=PARS(IZLAST)
 198       IF ( ABS(XT(IAXIS)-X(IAXIS)) .GT.1.) THEN
 199          TEMP = CL
 200          CL = -CH
 201          CH = -TEMP
 202       ENDIF
 203       IERR=0
 204       GO TO 999
 205 C
 206    50 CONTINUE
 207 **
 208       IF(ISH.EQ.13) THEN
 209          CL=-PARS(IAXIS)
 210          CH=PARS(IAXIS)
 211          IERR=0
 212          GOTO 999
 213       ENDIF
 214 **
 215       IF(ISH.EQ.14) THEN
 216 C     for hyperboloid, use escribed cylinder
 217          CH = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
 218          CL = -CH
 219          IERR=0
 220          GOTO 999
 221       ENDIF
 222 **
 223       IF(ISH.GT.10.AND.ISH.NE.NSCTUB)GO TO 999
 224       IF(ABS(XT(IAXIS)).GT.0.01) GO TO 70
 225 C
 226 C         Z AXIS PERPENDICULAR TO IAXIS. ASSUME COMPLETE TUBE OR
 227 C         CONE (I.E. IGNORE PHI SEGMENTATION).
 228 C
 229       IF(ISH.GT.6.AND.ISH.NE.NSCTUB) GO TO 60
 230 C
 231       CL=-PARS(2)
 232       CH=PARS(2)
 233       IERR=0
 234       IF(ISH.EQ.6)THEN
 235         RMIN=PARS(1)
 236         RMAX=PARS(2)
 237         IF(IROT.NE.0)THEN
 238          IF(Q(JROT+15).EQ.0.)THEN
 239            PHI1=(PARS(4)+Q(JROT+12))*DEGRAD
 240            PHI2=(PARS(5)+Q(JROT+12))*DEGRAD
 241          ELSEIF(Q(JROT+15).EQ.180.)THEN
 242            PHI1=(PARS(4)+Q(JROT+12)-(PARS(5)-PARS(4)))*DEGRAD
 243            PHI2=(PARS(5)+Q(JROT+12)-(PARS(5)-PARS(4)))*DEGRAD
 244          ELSE
 245            GOTO 999
 246          ENDIF
 247         ELSE
 248          PHI1=PARS(4)*DEGRAD
 249          PHI2=PARS(5)*DEGRAD
 250         ENDIF
 251         IF(IAXIS.EQ.1)THEN
 252          IF(PHI1.GE.0..AND.PHI2.LE.PI)THEN
 253           XMIN1=RMIN*COS(PHI2)
 254           XMIN2=RMAX*COS(PHI2)
 255           CL=MIN(XMIN1,XMIN2)
 256           XMAX1=RMIN*COS(PHI1)
 257           XMAX2=RMAX*COS(PHI1)
 258           CH=MAX(XMAX1,XMAX2)
 259          ELSEIF(PHI1.GE.PI.AND.PHI2.LE.TWOPI.OR.
 260      +          PHI1.GE.-PI.AND.PHI2.LE.0.)THEN
 261           XMIN1=RMIN*COS(PHI1)
 262           XMIN2=RMAX*COS(PHI1)
 263           CL=MIN(XMIN1,XMIN2)
 264           XMAX1=RMIN*COS(PHI2)
 265           XMAX2=RMAX*COS(PHI2)
 266           CH=MAX(XMAX1,XMAX2)
 267          ELSEIF(PHI1.LT.0..AND.PHI2.GT.0..AND.
 268      +         (PHI2-PHI1).LE.PI)THEN
 269           XMIN1=RMIN*COS(PHI2)
 270           XMIN2=RMIN*COS(PHI1)
 271           CL1=MIN(XMIN1,XMIN2)
 272           XMIN3=RMAX*COS(PHI2)
 273           XMIN4=RMAX*COS(PHI1)
 274           CL2=MIN(XMIN3,XMIN4)
 275           CL=MIN(CL1,CL2)
 276           CH=RMAX
 277          ELSEIF(PHI1.LT.PI.AND.PHI2.GT.PI.AND.
 278      +         (PHI2-PHI1).LE.PI)THEN
 279           CL=-RMAX
 280           XMAX1=RMIN*COS(PHI2)
 281           XMAX2=RMIN*COS(PHI1)
 282           CH1=MAX(XMAX1,XMAX2)
 283           XMAX3=RMAX*COS(PHI2)
 284           XMAX4=RMAX*COS(PHI1)
 285           CH2=MAX(XMAX3,XMAX4)
 286           CH=MAX(CH1,CH2)
 287          ENDIF
 288         ELSEIF(IAXIS.EQ.2)THEN
 289          IF(PHI1.GE.(-PI*.5).AND.PHI2.LE.(PI*.5))THEN
 290           YMIN1=RMIN*SIN(PHI1)
 291           YMIN2=RMAX*SIN(PHI1)
 292           CL=MIN(YMIN1,YMIN2)
 293           YMAX1=RMIN*SIN(PHI2)
 294           YMAX2=RMAX*SIN(PHI2)
 295           CH=MAX(YMAX1,YMAX2)
 296          ELSEIF(PHI1.GE.(PI*.5).AND.PHI2.LE.(PI*3*.5))THEN
 297           YMIN1=RMIN*SIN(PHI2)
 298           YMIN2=RMAX*SIN(PHI2)
 299           CL=MIN(YMIN1,YMIN2)
 300           YMAX1=RMIN*SIN(PHI1)
 301           YMAX2=RMAX*SIN(PHI1)
 302           CH=MAX(YMAX1,YMAX2)
 303          ELSEIF(PHI1.LT.(PI*.5).AND.PHI2.GT.(PI*.5).AND.
 304      +         (PHI2-PHI1).LE.PI)THEN
 305           YMIN1=RMIN*SIN(PHI2)
 306           YMIN2=RMIN*SIN(PHI1)
 307           CL1=MIN(YMIN1,YMIN2)
 308           YMIN3=RMAX*SIN(PHI2)
 309           YMIN4=RMAX*SIN(PHI1)
 310           CL2=MIN(YMIN3,YMIN4)
 311           CL=MIN(CL1,CL2)
 312           CH=RMAX
 313          ELSEIF(((PHI1.LT.(PI*3*.5).AND.PHI2.GT.(PI*3*.5)).OR.
 314      +          (PHI1.LT.-(PI*.5).AND.PHI2.GT.-(PI*.5)))
 315      +          .AND.(PHI2-PHI1).LE.PI)THEN
 316           CL=-RMAX
 317           YMAX1=RMIN*SIN(PHI2)
 318           YMAX2=RMIN*SIN(PHI1)
 319           CH1=MAX(YMAX1,YMAX2)
 320           YMAX3=RMAX*SIN(PHI2)
 321           YMAX4=RMAX*SIN(PHI1)
 322           CH2=MAX(YMAX3,YMAX4)
 323           CH=MAX(CH1,CH2)
 324          ENDIF
 325         ENDIF
 326       ENDIF
 327 C
 328       GO TO 999
 329 C
 330    60 CONTINUE
 331 C
 332       RM=PARS(3)
 333       IF(PARS(5).GT.PARS(3)) RM=PARS(5)
 334 C
 335       CL=-RM
 336       CH=RM
 337       IERR=0
 338 C
 339       GO TO 999
 340 C
 341    70 CONTINUE
 342 C
 343 C           ARBITRARY ROTATION.
 344 C
 345       DZ=PARS(3)
 346       RM=PARS(2)
 347       IF(ISH.EQ.13) THEN
 348 **
 349 **       approxime to a cylinder whit radius
 350 **       equal to the ellipse major axis
 351 **
 352          IF(PARS(1).GT.RM) RM=PARS(1)
 353          GOTO 80
 354       ENDIF
 355 **
 356       IF(ISH.EQ.14) THEN
 357         RM = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
 358         GO TO 80
 359       ENDIF
 360 *
 361       IF(ISH.EQ.NSCTUB) THEN
 362         S1 = (1.0-PARS(8))*(1.0+PARS(8))
 363         IF( S1 .GT. 0.0) S1 = SQRT(S1)
 364         S2 = (1.0-PARS(11))*(1.0+PARS(11))
 365         IF( S2 .GT. 0.0) S2 = SQRT(S2)
 366         IF( S2 .GT. S1 ) S1 = S2
 367         DZ = DZ+RM*S1
 368       ENDIF
 369       IF(ISH.LE.6) GO TO 80
 370 C
 371       DZ=PARS(1)
 372       RM=PARS(3)
 373       IF(PARS(5).GT.RM) RM=PARS(5)
 374 C
 375    80 CONTINUE
 376 C
 377       COST=ABS(XT(IAXIS))
 378       SINT=(1+COST)*(1-COST)
 379       IF(SINT.GT.0.0) SINT=SQRT(SINT)
 380 C
 381       CH=COST*DZ+SINT*RM
 382       CL=-CH
 383       IERR=0
 384 C
 385       GO TO 999
 386    90 CONTINUE
 387 C
 388 C           SPHERE - ASSUME COMPLETE SPHERE, TAKE OUTER RADIUS.
 389 C
 390       IERR=0
 391       CL=-PARS(2)
 392       CH=PARS(2)
 393 C
 394   999 CONTINUE
 395       END