GEANT321/ggeom/gflphi.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 GFLPHI(ISH,IROT,DX,PARS,CL,CH,IERR)
  13 C.
  14 C.    **********************************************************
  15 C.    *                                                        *
  16 C.    *    ROUTINE TO FIND THE PHI LIMITS OF THE OBJECT SHAPE  *
  17 C.    *    ISH ROTATED BY MATRIX IROT AND DISPLACED BY VECTOR  *
  18 C.    *    DX. IT HAS NPAR PARAMTERS IN THE ARRAY PARS. THE    *
  19 C.    *    LOWER LIMIT IS RETURNED IN CL AND THE HIGHER IN CH. *
  20 C.    *    NOTE THE OBJECT IS CONTAINED IN THE RANGE OF        *
  21 C.    *    INCREASING PHI FROM CL TO CH THOUGH CL AND CH ARE   *
  22 C.    *    FORCED TO LIE IN THE RANGE 0.0 TO 360.0 SO THAT THE *
  23 C.    *    VALUE OF CL CAN BE HIGHER THAN THAT OF CH.          *
  24 C.    *                                                        *
  25 C.    *    ==>Called by : GFCLIM                               *
  26 C.    *         Author  A.McPherson  *********                 *
  27 C.    *                                                        *
  28 C.    **********************************************************
  29 C.
  30 #include "geant321/gcbank.inc"
  31 #include "geant321/gconsp.inc"
  32 #include "geant321/gcshno.inc"
  33       DIMENSION DX(3),PARS(11),X(3),XT(3)
  34 C.
  35 C.          -------------------------------------------
  36 C.
  37       IERR=1
  38 C
  39       DXS=DX(1)*DX(1)+DX(2)*DX(2)
  40       IF(DXS.GT.0.0) DXS=SQRT(DXS)
  41 C
  42       IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 40
  43 C
  44 C
  45 C           CUBOIDS, TRAPEZOIDS, PARALLELEPIPEDS.
  46 C
  47       IERR=0
  48       CL=0.0
  49       CH=360.0
  50 C
  51 C           IF IN DOUBT SET IT TO FULL RANGE.
  52 C
  53       IF(DXS.LE.0.0) GO TO 999
  54 C
  55       PHC=90.
  56       IF(DX(1).NE.0.)PHC=ATAN2(DX(2),DX(1))*RADDEG
  57       IF(PHC.LT.0.0) PHC=PHC+360
  58       PL=0.0
  59       PH=0.0
  60 C
  61       DO 30 IP=1,8
  62 C
  63 C           THIS IS A LOOP OVER THE 8 CORNERS.
  64 C           FIRST FIND THE LOCAL COORDINATES.
  65 C
  66       IF(ISH.EQ.28) THEN
  67 C
  68 C            General twisted trapezoid.
  69 C
  70          IL=(IP+1)/2
  71          I0=IL*4+11
  72          IS=(IP-IL*2)*2+1
  73          X(3)=PARS(1)*IS
  74          X(1)=PARS(I0)+PARS(I0+2)*X(3)
  75          X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
  76          GO TO 20
  77 C
  78       ENDIF
  79 C
  80       IP3=ISH+2
  81       IF(ISH.EQ.10) IP3=3
  82       IF(ISH.EQ.4) IP3=1
  83       X(3)=PARS(IP3)
  84       IF(IP.LE.4) X(3)=-X(3)
  85       IP2=3
  86       IF(ISH.GT.2.AND.X(3).GT.0.0) IP2=4
  87       IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
  88       IF(ISH.EQ.4) IP2=4
  89       IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
  90       X(2)=PARS(IP2)
  91       IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
  92       IP1=1
  93       IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
  94       IF(ISH.EQ.4) IP1=5
  95       IF(ISH.EQ.4.AND.X(3).GT.0.0) IP1=IP1+4
  96       IF(ISH.EQ.4.AND.X(2).GT.0.0) IP1=IP1+1
  97       X(1)=PARS(IP1)
  98       IF(MOD(IP,2).EQ.1) X(1)=-X(1)
  99 C
 100       IF(ISH.NE.10) GO TO 10
 101       X(1)=X(1)+X(2)*PARS(4)+X(3)*PARS(5)
 102       X(2)=X(2)+X(3)*PARS(6)
 103    10 CONTINUE
 104 C
 105       IF(ISH.NE.4) GO TO 20
 106       IP4=7
 107       IF(X(3).GT.0.0) IP4=11
 108       X(1)=X(1)+X(2)*PARS(IP4)+X(3)*PARS(2)
 109       X(2)=X(2)+X(3)*PARS(3)
 110    20 CONTINUE
 111 C
 112 C          ROTATE.
 113 C
 114       JROT=LQ(JROTM-IROT)
 115       XT(1)=X(1)
 116       XT(2)=X(2)
 117       XT(3)=X(3)
 118       IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
 119 C
 120       XPT=DXS+(DX(1)*XT(1)+DX(2)*XT(2))/DXS
 121       YPT=(DX(1)*XT(2)-DX(2)*XT(1))/DXS
 122 C
 123       IF(YPT.EQ.0.0.AND.XPT.EQ.0.0) GO TO 999
 124       P=ATAN2(YPT,XPT)
 125       IF(P.GT.PI) P=P-PI*2.0
 126       IF(P.LT.PL) PL=P
 127       IF(P.GT.PH) PH=P
 128 C
 129 C
 130    30 CONTINUE
 131 C
 132 C
 133       IF(PH-PL.GT.PI) GO TO 999
 134       CL=PHC+PL*RADDEG
 135       CH=PHC+PH*RADDEG
 136 C
 137       SG = SIGN(1.0,CL)
 138       CL = MOD( ABS(CL),360.0 )
 139       IF(SG.LE.0.0) CL=360.-CL
 140       SG=SIGN(1.0,CH)
 141       CH = MOD( ABS(CH),360.0 )
 142       IF(SG.LE.0.0) CH=360.-CH
 143 C
 144       GO TO 999
 145 C
 146    40 CONTINUE
 147       IF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13.AND.ISH.NE.14)GO TO 80
 148 C
 149 C             TUBES AND CONES.
 150 C
 151       IERR=0
 152       CL=0.0
 153       CH=360.0
 154 C
 155 C             WHEN IN DOUBT SET TO FULL RANGE.
 156 C
 157       RM=PARS(2)
 158       IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 50
 159 **
 160       IF(ISH.EQ.13) THEN
 161 **
 162 **       approxime to a cylinder whit radius
 163 **       equal to the ellipse major axis
 164 **
 165          IF(PARS(1).GT.RM) RM=PARS(1)
 166          GOTO 50
 167       ENDIF
 168       IF(ISH.EQ.14) THEN
 169         RM = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
 170         GO TO 50
 171       ENDIF
 172       RM=PARS(3)
 173       IF(PARS(5).GT.RM) RM=PARS(5)
 174 C
 175    50 CONTINUE
 176 C
 177       IF(DXS.GT.RM) GO TO 70
 178       IF(ISH.EQ.5.OR.ISH.EQ.7.OR.ISH.EQ.14) GO TO 999
 179       IF(ISH.EQ.13) GOTO 999
 180 *                 Here we treat the CONS
 181       IF(IROT.EQ.0) THEN
 182 *                 This is the simple case, no rotation
 183 *                 Compute the position of the limits on
 184 *                 the X-Y plane.
 185       PHIMIN=PARS(4)*DEGRAD
 186       PHIMAX=PARS(5)*DEGRAD
 187       DDX1 = DX(1)+RM*COS(PHIMIN)
 188       DDY1 = DX(2)+RM*SIN(PHIMIN)
 189       DDX2 = DX(1)+RM*COS(PHIMAX)
 190       DDY2 = DX(2)+RM*SIN(PHIMAX)
 191       CL = ATAN2(DDY1,DDX1)*RADDEG
 192       CH = ATAN2(DDY2,DDX2)*RADDEG
 193       ELSE
 194 *                 Rotated tubes might be more difficult
 195 *                 Just leave it for later
 196       CONTINUE
 197       ENDIF
 198       IF(ISH.LE.6) GO TO 60
 199       CL=PARS(6)
 200       CH=PARS(7)
 201    60 CONTINUE
 202 C
 203       IF(IROT.EQ.0) GO TO 65
 204       JROT=LQ(JROTM-IROT)
 205       IF(Q(JROT+15).NE.0.0.AND.Q(JROT+15).NE.180.0) GO TO 999
 206 C
 207       PHX=Q(JROT+12)
 208       PHY=Q(JROT+14)
 209       IF(PHY.LT.PHX) PHY=PHY+360.0
 210       ISPH=1
 211       IF(PHY-PHX.GT.180.0) ISPH=-1
 212       CL=ISPH*CL+PHX
 213       CH=ISPH*CH+PHX
 214       IF(ISPH.EQ.1) GO TO 65
 215       CHT=CH
 216       CH=CL
 217       CL=CHT
 218 C
 219    65 CONTINUE
 220 C
 221       SG=SIGN(1.0,CL)
 222       CL = MOD( ABS(CL),360.0 )
 223       IF(SG.LE.0.0) CL=360.-CL
 224       SG=SIGN(1.0,CH)
 225       CH = MOD( ABS(CH),360.0 )
 226       IF(SG.LE.0.0) CH=360.-CH
 227 C
 228       GO TO 999
 229 C
 230    70 CONTINUE
 231 C
 232 C            DISPLACEMENT GREATER THAN MAXIMUM RADIUS SO
 233 C            ASSUME COMPLETE TUBE TO GENERATE 'WORST CASE'.
 234 C
 235       DZ=PARS(3)
 236       IF(ISH.EQ.NSCTUB) THEN
 237         S1 = (1.0-PARS(8))*(1.0+PARS(8))
 238         IF( S1 .GT. 0.0) S1 = SQRT(S1)
 239         S2 = (1.0-PARS(11))*(1.0+PARS(11))
 240         IF( S2 .GT. 0.0) S2 = SQRT(S2)
 241         IF( S2 .GT. S1 ) S1 = S2
 242         DZ = DZ+RM*S1
 243       ELSEIF(ISH.GT.6.AND.ISH.NE.13.AND.ISH.NE.14) THEN
 244         DZ=PARS(1)
 245       ENDIF
 246 C
 247       X(1)=0.0
 248       X(2)=0.0
 249       X(3)=1.0
 250 C
 251 C                    LOCAL Z AXIS.
 252 C
 253       JROT=LQ(JROTM-IROT)
 254       XT(1)=X(1)
 255       XT(2)=X(2)
 256       XT(3)=X(3)
 257       IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
 258 C
 259       COST=ABS(DX(1)*XT(1)+DX(2)*XT(2))
 260       COS2=ABS(DX(1)*XT(2)-DX(2)*XT(1))
 261       SINT=(DXS+COST)*(DXS-COST)
 262       SIN2=(DXS+COS2)*(DXS-COS2)
 263       IF(SINT.GT.0.0) SINT=SQRT(SINT)
 264       IF(SIN2.GT.0.0) SIN2=SQRT(SIN2)
 265 C
 266       XPT=DXS-(COST*DZ+SINT*RM)/DXS
 267 C
 268       IF(XPT.LE.0.0) GO TO 999
 269       YPT=(SIN2*RM+COS2*DZ)/DXS
 270       DP=ATAN(YPT/XPT)
 271 C
 272       P0=ATAN2(DX(2),DX(1))
 273       CL=(P0-DP)*RADDEG
 274       CH=(P0+DP)*RADDEG
 275 C
 276       SG=SIGN(1.0,CL)
 277       CL = MOD( ABS(CL),360.0 )
 278       IF(SG.LE.0.0) CL=360.-CL
 279       SG=SIGN(1.0,CH)
 280       CH = MOD( ABS(CH),360.0 )
 281       IF(SG.LE.0.0) CH=360.-CH
 282 C
 283       GO TO 999
 284 C
 285    80 CONTINUE
 286       IF(ISH.GT.9) GO TO 999
 287 C
 288 C               SPHERE.
 289 C
 290       IERR=0
 291       CL=0.0
 292       CH=360.0
 293 C
 294       IF(IROT.NE.0.OR.DXS.GT.0.0) GO TO 90
 295 C
 296 C          UNROTATED AND CENTERED.
 297 C
 298       CL=PARS(5)
 299       CH=PARS(6)
 300 C
 301       SG=SIGN(1.0,CL)
 302       CL = MOD( ABS(CL),360.0 )
 303       IF(SG.LE.0.0) CL=360.-CL
 304       SG=SIGN(1.0,CH)
 305       CH = MOD( ABS(CH),360.0 )
 306       IF(SG.LE.0.0) CH=360.-CH
 307 C
 308       GO TO 999
 309 C
 310    90 CONTINUE
 311 C
 312 C            ROTATED OR NOT CENTERED.
 313 C
 314       IF(DXS.LT.PARS(2)) GO TO 999
 315       P0=ATAN2(DX(2),DX(1))
 316       DP=ASIN(PARS(2)/DXS)
 317       CL=(P0-DP)*RADDEG
 318       CH=(P0+DP)*RADDEG
 319 C
 320       SG=SIGN(1.0,CL)
 321       CL = MOD( ABS(CL),360.0 )
 322       IF(SG.LE.0.0) CL=360.-CL
 323       SG=SIGN(1.0,CH)
 324       CH = MOD( ABS(CH),360.0 )
 325       IF(SG.LE.0.0) CH=360.-CH
 326 C
 327   999 CONTINUE
 328       END