GEANT321/ggeom/gflrad.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:48  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.28  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GFLRAD(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 : GFCLIM                                       *
  25 C.    *         Author  A.McPherson  *********                         *
  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(11),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       IF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13.AND.ISH.NE.14)GO TO 80
 136 C
 137 C             TUBES AND CONES.
 138 C
 139       IP3=3
 140       IF(ISH.GT.6.AND.ISH.NE.NSCTUB.AND.ISH.NE.13.AND.ISH.NE.14) IP3=1
 141       DZ=PARS(IP3)
 142       R=PARS(2)
 143       IF(ISH.EQ.NSCTUB) THEN
 144         S1 = (1.0-PARS(8))*(1.0+PARS(8))
 145         IF( S1 .GT. 0.0) S1 = SQRT(S1)
 146         S2 = (1.0-PARS(11))*(1.0+PARS(11))
 147         IF( S2 .GT. 0.0) S2 = SQRT(S2)
 148         IF( S2 .GT. S1 ) S1 = S2
 149         DZ = DZ+R*S1
 150       ENDIF
 151 **
 152       IF(ISH.EQ.13) THEN
 153 **
 154 **       APPROXIME TO A CYLINDER WHIT RADIUS
 155 **       EQUAL TO THE ELLIPSE MAJOR AXIS
 156 **
 157          RMN=0.0
 158          IF(PARS(1).GT.R) R=PARS(1)
 159          GOTO 50
 160       ENDIF
 161       RMN=PARS(1)
 162 *
 163       IF(ISH.EQ.14) THEN
 164         R = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
 165         GO TO 50
 166       ENDIF
 167 C
 168       IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 50
 169 C
 170       R=PARS(3)
 171       IF(PARS(5).GT.R) R=PARS(5)
 172       RMN=PARS(2)
 173       IF(PARS(4).LT.RMN) RMN=PARS(4)
 174 C
 175    50 CONTINUE
 176 C
 177 C          ROTATE THE LOCAL Z AXIS.
 178 C
 179       X(1)=0.0
 180       X(2)=0.0
 181       X(3)=1.0
 182       JROT=LQ(JROTM-IROT)
 183       XT(1)=X(1)
 184       XT(2)=X(2)
 185       XT(3)=X(3)
 186       IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
 187 C
 188 C          COMPUTE RMIN AND RMAX ASSUMING COMPLETE TUBE HALF
 189 C          LENGTH DZ AND RADIUS R.
 190 C
 191       ST2=1.0
 192       IF(IAXIS.EQ.4) ST2=(1+XT(3))*(1-XT(3))
 193       DR=SQRT(DZ*DZ*ST2+R*R)
 194       CL=DXS-DR
 195       IF(CL.LT.0.0) CL=0.0
 196       CH=DXS+DR
 197       IF(IROT.EQ.0.AND.DXS.LT.1.0E-05) CL=RMN
 198       IERR=0
 199 C
 200       GO TO 999
 201 C
 202    80 CONTINUE
 203       IF(ISH.GT.9) GO TO 999
 204 C
 205 C           SPHERE.
 206 C
 207       CL=DXS-PARS(2)
 208       IF(CL.LT.0.0) CL=0.0
 209       CH=DXS+PARS(2)
 210       IF(IAXIS.EQ.5.AND.DXS.LT.1.0E-05) CL=PARS(1)
 211       IERR=0
 212 C
 213   999 CONTINUE
 214       END