GEANT321/ggeom/gvdrad.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:56  cernlib
   9 * Geant
  10 *
  11 *
  12 #include "geant321/pilot.h"
  13 *CMZ :  3.21/02 21/07/94  11.48.23  by  S.Giani
  14 *-- Author :
  15       SUBROUTINE GVDRAD(IAXIS,ISH,IROT,DX,PARS,CL,CH,IERR)
  16 C.
  17 C.    ******************************************************************
  18 C.    *                                                                *
  19 C.    *    ROUTINE TO COMPUTE THE LIMITS IN R FOR THE SHAPE ISH        *
  20 C.    *    DISPLACED BY THE VECTOR DX AND ROTATED BY THE MATRIX IROT.  *
  21 C.    *    IF IAXIS = 4 THE R IS THE XY PLANE R, IF IAXIS = 5 IT IS    *
  22 C.    *    THE 3 DINEMSIONAL SPACE R. THE SHAPE HAS NPAR PARAMETERS    *
  23 C.    *    IN THE ARRAY PARS. THE LOWER LIMIT IS RETURNED IN CL AND    *
  24 C.    *    THE HIGHER IN CH. IF THE CALCULATION CANNOT BE PERFORMED    *
  25 C.    *    IERR IS SET TO 1 OTHERWISE IT IS SET TO 0.                  *
  26 C.    *                                                                *
  27 C.    *    ==>Called by : GVDLIM                                       *
  28 C.    *         Author  S.Giani  *********                             *
  29 C.    *                                                                *
  30 C.    ******************************************************************
  31 C.
  32 #include "geant321/gcbank.inc"
  33 #include "geant321/gconsp.inc"
  34 #include "geant321/gcshno.inc"
  35       DIMENSION DX(3),PARS(100),X(3),XT(3)
  36 C.
  37 C.           --------------------------------------------------
  38 C.
  39       IERR=1
  40 C
  41 C            FIRST CALCULATE THE LENGTH OF THE DISPLACEMENT OF THE
  42 C            ORIGIN.
  43 C
  44       DXS=DX(1)*DX(1)+DX(2)*DX(2)
  45       IF(IAXIS.EQ.5) DXS=DXS+DX(3)*DX(3)
  46       IF(DXS.GT.0.0) DXS=SQRT(DXS)
  47 C
  48       IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 40
  49 C
  50 C          CUBOIDS, TRAPEZOIDS, PARALLELEPIPEDS.
  51 C
  52       CH=0.0
  53       CL=DXS
  54 C
  55       DO 30 IP=1,8
  56 C
  57 C           THIS IS A LOOP OVER THE 8 CORNERS.
  58 C           FIRST FIND THE LOCAL COORDINATES.
  59 C
  60       IF(ISH.EQ.28) THEN
  61 C
  62 C            General twisted trapezoid.
  63 C
  64          IL=(IP+1)/2
  65          I0=IL*4+11
  66          IS=(IP-IL*2)*2+1
  67          X(3)=PARS(1)*IS
  68          X(1)=PARS(I0)+PARS(I0+2)*X(3)
  69          X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
  70          GO TO 20
  71 C
  72       ENDIF
  73 C
  74       IP3=ISH+2
  75       IF(ISH.EQ.10) IP3=3
  76       IF(ISH.EQ.4) IP3=1
  77       X(3)=PARS(IP3)
  78       IF(IP.LE.4) X(3)=-X(3)
  79       IP2=3
  80       IF(ISH.GT.2.AND.X(3).GT.0.0) IP2=4
  81       IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
  82       IF(ISH.EQ.4) IP2=4
  83       IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
  84       X(2)=PARS(IP2)
  85       IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
  86       IP1=1
  87       IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
  88       IF(ISH.EQ.4) IP1=5
  89       IF(ISH.EQ.4.AND.X(3).GT.0.0) IP1=IP1+4
  90       IF(ISH.EQ.4.AND.X(2).GT.0.0) IP1=IP1+1
  91       X(1)=PARS(IP1)
  92       IF(MOD(IP,2).EQ.1) X(1)=-X(1)
  93 C
  94       IF(ISH.NE.10) GO TO 10
  95       X(1)=X(1)+X(2)*PARS(4)+X(3)*PARS(5)
  96       X(2)=X(2)+X(3)*PARS(6)
  97    10 CONTINUE
  98 C
  99       IF(ISH.NE.4) GO TO 20
 100       IP4=7
 101       IF(X(3).GT.0.0) IP4=11
 102       X(1)=X(1)+X(2)*PARS(IP4)+X(3)*PARS(2)
 103       X(2)=X(2)+X(3)*PARS(3)
 104    20 CONTINUE
 105 C
 106 C          ROTATE.
 107 C
 108       JROT=LQ(JROTM-IROT)
 109       XT(1)=X(1)
 110       XT(2)=X(2)
 111       XT(3)=X(3)
 112       IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
 113 C
 114 C          NOW COMPUTE RMIN = PROJECTED R ON DX AND RMAX = R
 115 C          AND UPDATE LIMITS IF NECESSARY.
 116 C
 117       R2=(XT(1)+DX(1))**2+(XT(2)+DX(2))**2
 118       IF(IAXIS.EQ.5) R2=R2+(XT(3)+DX(3))**2
 119       R=SQRT(R2)
 120       IF(R.GT.CH) CH=R
 121 C
 122       IF(CL.LE.0.0) GO TO 30
 123 C
 124       XPT=DX(1)*XT(1)+DX(2)*XT(2)
 125       IF(IAXIS.EQ.5) XPT=XPT+DX(3)*XT(3)
 126       IF(DXS.LE.1.0E-05) GO TO 30
 127       RMN=DXS+XPT/DXS
 128       IF(RMN.LT.CL) CL=RMN
 129 C
 130    30 CONTINUE
 131 C
 132       IF(CL.LE.0.0) CL=0.0
 133 C
 134       IERR=0
 135       GO TO 999
 136 C
 137    40 CONTINUE
 138 C
 139 C          POLYGONES AND POLYCONES
 140 C
 141       IF(ISH.EQ.11.AND.IAXIS.EQ.4)THEN
 142        NZLAST=PARS(4)
 143        IZLAST=2+3*NZLAST
 144        CLZ=PARS(5)
 145        CHZ=PARS(IZLAST)
 146        DZ2=ABS(CHZ-CLZ)
 147        DZ=DZ2*.5
 148        TMPRAD=0.
 149        TMPMIN=100000.
 150        DO 145 I=7,IZLAST+2,3
 151          IF(PARS(I).GT.TMPRAD)TMPRAD=PARS(I)
 152          IF(PARS(I-1).LT.TMPMIN)TMPMIN=PARS(I-1)
 153  145   CONTINUE
 154        RMN=TMPMIN
 155        PHIMIN=PARS(1)
 156        PHIMAX=PHIMIN+PARS(2)
 157        AANG=ABS(PHIMAX-PHIMIN)
 158        NANG=PARS(3)
 159        AATMAX=NANG*360./AANG
 160        LATMAX=AATMAX
 161        ALA=AATMAX-LATMAX
 162        IF(ALA.GT..5)LATMAX=LATMAX+1
 163        AFINV=1./COS(PI/LATMAX)
 164        FINV=ABS(AFINV)
 165        R=TMPRAD*FINV
 166        GOTO 50
 167       ELSEIF(ISH.EQ.12.AND.IAXIS.EQ.4)THEN
 168        NZLAST=PARS(3)
 169        IZLAST=1+3*NZLAST
 170        CLZ=PARS(4)
 171        CHZ=PARS(IZLAST)
 172        DZ2=ABS(CHZ-CLZ)
 173        DZ=DZ2*.5
 174        TMPRAD=0.
 175        TMPMIN=100000.
 176        DO 146 I=6,IZLAST+2,3
 177          IF(PARS(I).GT.TMPRAD)TMPRAD=PARS(I)
 178          IF(PARS(I-1).LT.TMPMIN)TMPMIN=PARS(I-1)
 179  146   CONTINUE
 180        RMN=TMPMIN
 181        R=TMPRAD
 182        GOTO 50
 183       ENDIF
 184       IF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13.AND.ISH.NE.14)GO TO 80
 185 C
 186 C             TUBES AND CONES
 187 C
 188       IP3=3
 189       IF(ISH.GT.6.AND.ISH.NE.NSCTUB.AND.ISH.NE.13.AND.ISH.NE.14) IP3=1
 190       DZ=PARS(IP3)
 191       R=PARS(2)
 192       IF(ISH.EQ.NSCTUB) THEN
 193         S1 = (1.0-PARS(8))*(1.0+PARS(8))
 194         IF( S1 .GT. 0.0) S1 = SQRT(S1)
 195         S2 = (1.0-PARS(11))*(1.0+PARS(11))
 196         IF( S2 .GT. 0.0) S2 = SQRT(S2)
 197         IF( S2 .GT. S1 ) S1 = S2
 198         DZ = DZ+R*S1
 199       ENDIF
 200 **
 201       IF(ISH.EQ.13) THEN
 202 **
 203 **       APPROXIME TO A CYLINDER WHIT RADIUS
 204 **       EQUAL TO THE ELLIPSE MAJOR AXIS
 205 **
 206          RMN=0.0
 207          IF(PARS(1).GT.R) R=PARS(1)
 208          GOTO 50
 209       ENDIF
 210       RMN=PARS(1)
 211 *
 212       IF(ISH.EQ.14) THEN
 213         R = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
 214         GO TO 50
 215       ENDIF
 216 C
 217       IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 50
 218 C
 219       R=PARS(3)
 220       IF(PARS(5).GT.R) R=PARS(5)
 221       RMN=PARS(2)
 222       IF(PARS(4).LT.RMN) RMN=PARS(4)
 223 C
 224    50 CONTINUE
 225 C
 226 C          ROTATE THE LOCAL Z AXIS.
 227 C
 228       X(1)=0.0
 229       X(2)=0.0
 230       X(3)=1.0
 231       JROT=LQ(JROTM-IROT)
 232       XT(1)=X(1)
 233       XT(2)=X(2)
 234       XT(3)=X(3)
 235       IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
 236 C
 237 C          COMPUTE RMIN AND RMAX ASSUMING COMPLETE TUBE HALF
 238 C          LENGTH DZ AND RADIUS R.
 239 C
 240       ST2=1.0
 241       IF(IAXIS.EQ.4) ST2=(1+XT(3))*(1-XT(3))
 242       DR=SQRT(DZ*DZ*ST2+R*R)
 243       CL=DXS-DR
 244       IF(CL.LT.0.0) CL=0.0
 245       CH=DXS+DR
 246       IF(ABS(XT(3)).EQ.1.0.AND.DXS.LT.1.0E-05) CL=RMN
 247       IERR=0
 248 C
 249       GO TO 999
 250 C
 251    80 CONTINUE
 252       IF(ISH.GT.9) GO TO 999
 253 C
 254 C           SPHERE.
 255 C
 256       CL=DXS-PARS(2)
 257       IF(CL.LT.0.0) CL=0.0
 258       CH=DXS+PARS(2)
 259       IF(IAXIS.EQ.5.AND.DXS.LT.1.0E-05) CL=PARS(1)
 260       IERR=0
 261 C
 262   999 CONTINUE
 263       END