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