[u/mrichter/AliRoot.git] / GEANT321 / ggeom / gflthe.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:20:48  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.28  by  S.Giani
*-- Author :
      SUBROUTINE GFLTHE(ISH,IROT,DX,PARS,CL,CH,IERR)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *    ROUTINE TO COMPUTE THE THETA LIMITS FOR VOLUME OF SHAPE     *
C.    *    ISH ROTATED BY MATRIX IROT AND DISPLACED BY VECTOR DX.      *
C.    *    THE VOLUME HAS NPAR PARAMETERS IN THE ARRAY PARS. THE LOWER *
C.    *    LIMIT IS RETURNED IN CL THE HIGHER IN CH. IF THE            *
C.    *    CALCULATION CANNOT BE MADE IERR IS SET TO 1 OTHERWISE IT    *
C.    *    IS SET TO 0.                                                *
C.    *                                                                *
C.    *    ==>Called by : GFCLIM                                       *
C.    *         Author  A.McPherson  *********                         *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gcbank.inc"
#include "geant321/gconsp.inc"
#include "geant321/gcshno.inc"
      DIMENSION DX(3),PARS(11),X(3),XT(3),X1(3),X2(3),XT1(3),XT2(3)
C.
C.          ----------------------------------------------
C.
      IERR=1
C
      DXS=DX(1)*DX(1)+DX(2)*DX(2)
      IF(DXS.GT.0.0) DXS=SQRT(DXS)
C
      IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 130
C
C            CUBES, TRAPEZOIDS AND PARALLELEPIPEDS.
C
      IERR=0
      CL=0.0
      CH=180.0
C
      IF(DXS.LE.0.0) GO TO 999
C
      TH=90.
      IF(ABS(DX(3)).LT.1.0E-06)GO TO 5
      TH=ATAN(DXS/DX(3))*RADDEG
      IF(TH.LT.0.0) TH=TH+180.0
   5  TL=TH
C
      DO 50 IP=1,8
C
C           THIS IS A LOOP OVER THE 8 CORNERS.
C           FIRST FIND THE LOCAL COORDINATES.
C
      IF(ISH.EQ.28) THEN
C
C            General twisted trapezoid.
C
         IL=(IP+1)/2
         I0=IL*4+11
         IS=(IP-IL*2)*2+1
         X(3)=PARS(1)*IS
         X(1)=PARS(I0)+PARS(I0+2)*X(3)
         X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
         GO TO 20
C
      ENDIF
C
      IP3=ISH+2
      IF(ISH.EQ.10) IP3=3
      IF(ISH.EQ.4) IP3=1
      X(3)=PARS(IP3)
      IF(IP.LE.4) X(3)=-X(3)
      IP2=3
      IF(ISH.GT.2.AND.X(3).GT.0.0) IP2=4
      IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
      IF(ISH.EQ.4) IP2=4
      IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
      X(2)=PARS(IP2)
      IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
      IP1=1
      IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
      IF(ISH.EQ.4) IP1=5
      IF(ISH.EQ.4.AND.X(3).GT.0.0) IP1=IP1+4
      IF(ISH.EQ.4.AND.X(2).GT.0.0) IP1=IP1+1
      X(1)=PARS(IP1)
      IF(MOD(IP,2).EQ.1) X(1)=-X(1)
C
      IF(ISH.NE.10) GO TO 10
      X(1)=X(1)+X(2)*PARS(4)+X(3)*PARS(5)
      X(2)=X(2)+X(3)*PARS(6)
   10 CONTINUE
C
      IF(ISH.NE.4) GO TO 20
      IP4=7
      IF(X(3).GT.0.0) IP4=11
      X(1)=X(1)+X(2)*PARS(IP4)+X(3)*PARS(2)
      X(2)=X(2)+X(3)*PARS(3)
   20 CONTINUE
C
C          ROTATE.
C
      JROT=LQ(JROTM-IROT)
      XT(1)=X(1)
      XT(2)=X(2)
      XT(3)=X(3)
      IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
C
      Z=DX(3)+XT(3)
      R=(DX(1)+XT(1))**2+(DX(2)+XT(2))**2
      RT=R+(DX(3)+XT(3))**2
      IF(RT.LT.1.0E-10) GO TO 999
C
      IF(R.GT.0.0) GO TO 30
      IF(Z.GT.0.0) T=0.0
      IF(Z.LT.0.0) T=180.0
      GO TO 40
C
   30 CONTINUE
      T=90.0
      IF(ABS(Z).LT.1.0E-6) GO TO 40
C
      R=SQRT(R)
      T=ATAN(R/Z)*RADDEG
      IF(T.EQ.0.0.AND.Z.LT.0.0) T=180.0
      IF(T.LT.0.0) T=T+180.0
C
   40 CONTINUE
      IF(T.GT.TH) TH=T
      IF(T.LT.TL) TL=T
C
   50 CONTINUE
C
C           THETA LIMITS SET FROM THE POINTS NOW DO THE EDGES.
C
      DO 120 IL=1,12
C
C           FIND THE END POINT NUMBERS FOR EACH EDGE.
C
      IF(IL.GT.4) GO TO 60
      IPP1=1
      IF(IL.GT.2) IPP1=4
      IPP2=2
      IF(MOD(IL,2).EQ.1) IPP2=3
C
      GO TO 80
   60 CONTINUE
      IF(IL.LT.9) GO TO 70
C
      IPP1=5
      IF(IL.GT.10) IPP1=8
      IPP2=6
      IF(MOD(IL,2).EQ.1) IPP2=7
C
      GO TO 80
   70 CONTINUE
C
      IPP1=IL-4
      IPP2=IL
C
   80 CONTINUE
C
C           NOW GET THE POINTS AND ROTATE THEM.
C
      IF(ISH.EQ.28) THEN
C
C            General twisted trapezoid.
C
         ILL=IPP1
         IF(IPP1.EQ.3) ILL=4
         IF(IPP1.EQ.4) ILL=3
         I0=ILL*4+11
         X1(3)=PARS(1)
         IF(IPP1.LT.5) X1(3)=-X1(3)
         X1(1)=PARS(I0)+PARS(I0+2)*X1(3)
         X1(2)=PARS(I0+1)+PARS(I0+3)*X1(3)
         ILL=IPP2
         IF(IPP2.EQ.3) ILL=4
         IF(IPP2.EQ.4) ILL=3
         I0=ILL*4+11
         X2(3)=PARS(1)
         IF(IPP2.LT.5) X2(3)=-X2(3)
         X2(1)=PARS(I0)+PARS(I0+2)*X2(3)
         X2(2)=PARS(I0+1)+PARS(I0+3)*X2(3)
C
         GO TO 100
C
      ENDIF
C
      IP3=ISH+2
      IF(ISH.EQ.10) IP3=3
      IF(ISH.EQ.4) IP3=1
      X1(3)=PARS(IP3)
      IF(IPP1.LE.4) X1(3)=-X1(3)
      IP2=3
      IF(ISH.GT.2.AND.X1(3).GT.0.0) IP2=4
      IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
      IF(ISH.EQ.4) IP2=4
      IF(ISH.EQ.4.AND.X1(3).GT.0.0) IP2=8
      X1(2)=PARS(IP2)
      IF(MOD(IPP1+3,4).LT.2) X1(2)=-X1(2)
      IP1=1
      IF(ISH.NE.1.AND.ISH.NE.10.AND.X1(3).GT.0.0) IP1=2
      IF(ISH.EQ.4) IP1=5
      IF(ISH.EQ.4.AND.X1(3).GT.0.0) IP1=IP1+4
      IF(ISH.EQ.4.AND.X1(2).GT.0.0) IP1=IP1+1
      X1(1)=PARS(IP1)
      IF(MOD(IPP1,2).EQ.1) X1(1)=-X1(1)
C
      IP3=ISH+2
      IF(ISH.EQ.10) IP3=3
      IF(ISH.EQ.4) IP3=1
      X2(3)=PARS(IP3)
      IF(IPP2.LE.4) X2(3)=-X2(3)
      IP2=3
      IF(ISH.GT.2.AND.X2(3).GT.0.0) IP2=4
      IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
      IF(ISH.EQ.4) IP2=4
      IF(ISH.EQ.4.AND.X2(3).GT.0.0) IP2=8
      X2(2)=PARS(IP2)
      IF(MOD(IPP2+3,4).LT.2) X2(2)=-X2(2)
      IP1=1
      IF(ISH.NE.1.AND.ISH.NE.10.AND.X2(3).GT.0.0) IP1=2
      IF(ISH.EQ.4) IP1=5
      IF(ISH.EQ.4.AND.X2(3).GT.0.0) IP1=IP1+4
      IF(ISH.EQ.4.AND.X2(2).GT.0.0) IP1=IP1+1
      X2(1)=PARS(IP1)
      IF(MOD(IPP2,2).EQ.1) X2(1)=-X2(1)
C
      IF(ISH.NE.10) GO TO 90
      X1(1)=X1(1)+X1(2)*PARS(4)+X1(3)*PARS(5)
      X1(2)=X1(2)+X1(3)*PARS(6)
      X2(1)=X2(1)+X2(2)*PARS(4)+X2(3)*PARS(5)
      X2(2)=X2(2)+X2(3)*PARS(6)
   90 CONTINUE
C
      IF(ISH.NE.4) GO TO 100
      IP4=7
      IF(X1(3).GT.0.0) IP4=11
      X1(1)=X1(1)+X1(2)*PARS(IP4)+X1(3)*PARS(2)
      X1(2)=X1(2)+X1(3)*PARS(3)
      IP4=7
      IF(X2(3).GT.0.0) IP4=11
      X2(1)=X2(1)+X2(2)*PARS(IP4)+X2(3)*PARS(2)
      X2(2)=X2(2)+X2(3)*PARS(3)
  100 CONTINUE
C
C          ROTATE.
C
      JROT=LQ(JROTM-IROT)
      XT1(1)=X1(1)
      XT1(2)=X1(2)
      XT1(3)=X1(3)
      IF(IROT.NE.0) CALL GINROT(X1,Q(JROT+1),XT1)
      XT2(1)=X2(1)
      XT2(2)=X2(2)
      XT2(3)=X2(3)
      IF(IROT.NE.0) CALL GINROT(X2,Q(JROT+1),XT2)
C
C           NOW WE HAVE THE END POINTS IN THE MASTER SYSTEM.
C           FIND THE TANGENT POINT TO THE SET OF CONES OF CONSTANT
C           THETA.
C
      DS=(XT2(1)-XT1(1))**2+(XT2(2)-XT1(2))**2+(XT2(3)-XT1(3))**2
      IF(DS.LE.0.0) GO TO 120
C
      DS=SQRT(DS)
C
      X0=(XT2(1)+XT1(1))/2.0+DX(1)
      Y0=(XT2(2)+XT1(2))/2.0+DX(2)
      Z0=(XT2(3)+XT1(3))/2.0+DX(3)
      ALX=(XT2(1)-XT1(1))/DS
      ALY=(XT2(2)-XT1(2))/DS
      ALZ=(XT2(3)-XT1(3))/DS
C
      T=90.0
      IF(Z0.EQ.0.0.AND.ALZ.EQ.0.0) GO TO 110
C
      IF(ALX.EQ.0.0.AND.ALY.EQ.0.0) GO TO 120
C             THIS CHECKS WHETHER THE LINE IS PARALLEL TO THE
C             Z AXIS IN WHICH CASE THERE IS NO TANGENT AND
C             THE END POINTS DETERMINE THE THETA RANGE.
C
      SNUM=(X0*Z0*ALX+Y0*Z0*ALY+X0*X0*ALZ+Y0*Y0*ALZ)
      SDEN=(Z0*ALX*ALX-X0*ALX*ALZ+Z0*ALY*ALY-Y0*ALY*ALZ)
C
      IF(ABS(SNUM).GT.0.5*DS*ABS(SDEN)) GO TO 120
C
C           TANGENT EXIST BETWEEN THE TWO ENDS.
C
      S = SNUM/SDEN
      X0=X0+S*ALX
      Y0=Y0+S*ALY
      Z0=Z0+S*ALZ
C
      R=X0*X0+Y0*Y0
      RT=R+Z0*Z0
C
      IF(RT.LT.1.0E-10) GO TO 999
C
      IF(R.GT.0.0) R=SQRT(R)
      T=90.0
      IF(ABS(Z0).LT.1.0E-06) GO TO 110
C
      T=ATAN(R/Z0)*RADDEG
      IF(T.EQ.0.0.AND.Z0.LT.0.0) T=180.0
      IF(T.LT.0.0) T=T+180.0
C
  110 CONTINUE
C
      IF(T.LT.TL) TL=T
      IF(T.GT.TH) TH=T
C
C        CHECK FOR THE POSSIBILITY OF STRADDLING T=0.0 &/OR 180.0.
C
      C=X0*DX(1)+Y0*DX(2)
      IF(C.GT.0.0) GO TO 120
C
C          CHECK IF SAME SIGN OF Z.
C
      IF(Z0*DX(3).LT.0.0) GO TO 999
C
      T=0.0
      IF(Z0.LT.0.0) T=180.0
      IF(T.LT.TL) TL=T
      IF(T.GT.TH) TH=T
C
  120 CONTINUE
C
C            DONE SET THE LIMITS.
C
      CL=TL
      CH=TH
C
      GO TO 999
C
  130 CONTINUE
C
C          TUBES, CONES ETC.
C
      IF(IROT.NE.0.OR.DXS.GT.1.0E-05) GO TO 180
C
C          UNROTATED AND CENTERED.
C
      IF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13
     +           .AND.ISH.NE.14)GO TO 170
C
C               TUBES AND CONES.
C
      IERR=0
      DZ=PARS(3)
      RMN=PARS(1)
      RMX=PARS(2)
      IF(ISH.EQ.13) THEN
**
**         approxime to a cylinder whit radius
**         equal to the ellipse major axis
**
        RMN=0.0
        IF(PARS(1).GT.RMX) RMX=PARS(1)
        GOTO 140
      ENDIF
      IF(ISH.EQ.14) THEN
C   not really sure of the function of these... keep RMN=PARS(1)
        RMX = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
        GO TO 140
      ENDIF
      IF(ISH.EQ.NSCTUB) THEN
        S1 = (1.0-PARS(8))*(1.0+PARS(8))
        IF( S1 .GT. 0.0) S1 = SQRT(S1)
        S2 = (1.0-PARS(11))*(1.0+PARS(11))
        IF( S2 .GT. 0.0) S2 = SQRT(S2)
        IF( S2 .GT. S1 ) S1 = S2
        DZ = DZ+RMX*S1
      ENDIF
      IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 140
      DZ=PARS(1)
      RMN=PARS(2)
      RMX=PARS(3)
      IF(PARS(4).LT.RMN) RMN=PARS(4)
      IF(PARS(5).GT.RMX) RMX=PARS(5)
  140 CONTINUE
C
      IF(DZ.GT.ABS(DX(3))) GO TO 160
C
C          ALL FORWARD OR ALL BACK.
C
      DSH=DX(3)-DZ
      DLN=DX(3)+DZ
      IF(DX(3).GT.0.0) GO TO 150
      DSS=DSH
      DSH=DLN
      DLN=DSS
  150 CONTINUE
C
      CL=90.0
      CH=90.0
      IF(DLN.NE.0.0) CL=ATAN(RMN/DLN)*RADDEG
      IF(DSH.NE.0.0) CH=ATAN(RMX/DSH)*RADDEG
      IF(DX(3).GT.0.0) GO TO 999
      CS=CL
      CL=CH
      CH=CS
      IF(CH.EQ.0.0) CH=180.0
      IF(CL.LT.0.0) CL=CL+180.0
      IF(CH.LT.0.0) CH=CH+180.0
C
      GO TO 999
C
  160 CONTINUE
C
      CL=90.0
      CH=90.0
      IF(DZ+DX(3).NE.0.0) CL=ATAN(RMN/(DZ+DX(3)))*RADDEG
      IF(-DZ+DX(3).NE.0.0) CH=ATAN(RMN/(-DZ+DX(3)))*RADDEG
      IF(CH.LE.0.0) CH=CH+180.0
C
      GO TO 999
C
  170 CONTINUE
      IF(ISH.GT.9) GO TO 999
C
C           SPHERE.
C
      IERR=0
      CL=PARS(3)
      CH=PARS(4)
C
      GO TO 999
  180 CONTINUE
C
      IF(ISH.EQ.11.OR.ISH.EQ.12) GOTO 999
**
      RM=PARS(2)
      IF(ISH.EQ.9) GO TO 200
C
      DZ=PARS(3)
      IF(ISH.EQ.NSCTUB) THEN
        S1 = (1.0-PARS(8))*(1.0+PARS(8))
        IF( S1 .GT. 0.0) S1 = SQRT(S1)
        S2 = (1.0-PARS(11))*(1.0+PARS(11))
        IF( S2 .GT. 0.0) S2 = SQRT(S2)
        IF( S2 .GT. S1 ) S1 = S2
        DZ = DZ+RM*S1
      ENDIF
**
      IF(ISH.EQ.13) THEN
        IF(PARS(1).GT.RM) RM=PARS(1)
        GOTO 190
      ENDIF
**
      IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 190
C
      DZ=PARS(1)
      RM=PARS(3)
      IF(PARS(5).GT.RM) RM=PARS(5)
C
  190 CONTINUE
C
      RM=SQRT(RM**2+DZ**2)
C
  200 CONTINUE
C
      CL=0.0
      CH=180.0
      IERR=0
      RC=DXS**2+DX(3)**2
      IF(RC.GT.0.0) RC=SQRT(RC)
      IF(RM.GE.RC) GO TO 999
C
      TC=90.0
      IF(ABS(DX(3)).GT.0.0) TC=ATAN(DXS/DX(3))*RADDEG
      IF(TC.LT.0.0) TC=TC+180.0
C
      DT=ASIN(RM/RC)*RADDEG
      CL=TC-DT
      CH=TC+DT
      IF(CL.LT.0.0) CL=0.0
      IF(CH.GT.180.0) CH=180.0
C
  999 CONTINUE
      END
Commit	Line	Data
fe4da5cc	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 GFLTHE(ISH,IROT,DX,PARS,CL,CH,IERR)
	13	C.
	14	C. ******************************************************************
	15	C. * *
	16	C. * ROUTINE TO COMPUTE THE THETA LIMITS FOR VOLUME OF SHAPE *
	17	C. * ISH ROTATED BY MATRIX IROT AND DISPLACED BY VECTOR DX. *
	18	C. * THE VOLUME HAS NPAR PARAMETERS IN THE ARRAY PARS. THE LOWER *
	19	C. * LIMIT IS RETURNED IN CL THE HIGHER IN CH. IF THE *
	20	C. * CALCULATION CANNOT BE MADE IERR IS SET TO 1 OTHERWISE IT *
	21	C. * IS SET TO 0. *
	22	C. * *
	23	C. * ==>Called by : GFCLIM *
	24	C. * Author A.McPherson ********* *
	25	C. * *
	26	C. ******************************************************************
	27	C.
	28	#include "geant321/gcbank.inc"
	29	#include "geant321/gconsp.inc"
	30	#include "geant321/gcshno.inc"
	31	DIMENSION DX(3),PARS(11),X(3),XT(3),X1(3),X2(3),XT1(3),XT2(3)
	32	C.
	33	C. ----------------------------------------------
	34	C.
	35	IERR=1
	36	C
	37	DXS=DX(1)DX(1)+DX(2)DX(2)
	38	IF(DXS.GT.0.0) DXS=SQRT(DXS)
	39	C
	40	IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 130
	41	C
	42	C CUBES, TRAPEZOIDS AND PARALLELEPIPEDS.
	43	C
	44	IERR=0
	45	CL=0.0
	46	CH=180.0
	47	C
	48	IF(DXS.LE.0.0) GO TO 999
	49	C
	50	TH=90.
	51	IF(ABS(DX(3)).LT.1.0E-06)GO TO 5
	52	TH=ATAN(DXS/DX(3))*RADDEG
	53	IF(TH.LT.0.0) TH=TH+180.0
	54	5 TL=TH
	55	C
	56	DO 50 IP=1,8
	57	C
	58	C THIS IS A LOOP OVER THE 8 CORNERS.
	59	C FIRST FIND THE LOCAL COORDINATES.
	60	C
	61	IF(ISH.EQ.28) THEN
	62	C
	63	C General twisted trapezoid.
	64	C
65	IL=(IP+1)/2
66	I0=IL*4+11
67	IS=(IP-IL2)2+1
68	X(3)=PARS(1)*IS
69	X(1)=PARS(I0)+PARS(I0+2)*X(3)
70	X(2)=PARS(I0+1)+PARS(I0+3)*X(3)
71	GO TO 20
72	C
73	ENDIF
74	C
75	IP3=ISH+2
76	IF(ISH.EQ.10) IP3=3
77	IF(ISH.EQ.4) IP3=1
78	X(3)=PARS(IP3)
79	IF(IP.LE.4) X(3)=-X(3)
80	IP2=3
81	IF(ISH.GT.2.AND.X(3).GT.0.0) IP2=4
82	IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
83	IF(ISH.EQ.4) IP2=4
84	IF(ISH.EQ.4.AND.X(3).GT.0.0) IP2=8
85	X(2)=PARS(IP2)
86	IF(MOD(IP+3,4).LT.2) X(2)=-X(2)
87	IP1=1
88	IF(ISH.NE.1.AND.ISH.NE.10.AND.X(3).GT.0.0) IP1=2
89	IF(ISH.EQ.4) IP1=5
90	IF(ISH.EQ.4.AND.X(3).GT.0.0) IP1=IP1+4
91	IF(ISH.EQ.4.AND.X(2).GT.0.0) IP1=IP1+1
92	X(1)=PARS(IP1)
93	IF(MOD(IP,2).EQ.1) X(1)=-X(1)
94	C
95	IF(ISH.NE.10) GO TO 10
96	X(1)=X(1)+X(2)PARS(4)+X(3)PARS(5)
97	X(2)=X(2)+X(3)*PARS(6)
98	10 CONTINUE
99	C
100	IF(ISH.NE.4) GO TO 20
101	IP4=7
102	IF(X(3).GT.0.0) IP4=11
103	X(1)=X(1)+X(2)PARS(IP4)+X(3)PARS(2)
104	X(2)=X(2)+X(3)*PARS(3)
105	20 CONTINUE
106	C
107	C ROTATE.
108	C
109	JROT=LQ(JROTM-IROT)
110	XT(1)=X(1)
111	XT(2)=X(2)
112	XT(3)=X(3)
113	IF(IROT.NE.0) CALL GINROT(X,Q(JROT+1),XT)
114	C
115	Z=DX(3)+XT(3)
116	R=(DX(1)+XT(1))2+(DX(2)+XT(2))2
117	RT=R+(DX(3)+XT(3))**2
118	IF(RT.LT.1.0E-10) GO TO 999
119	C
120	IF(R.GT.0.0) GO TO 30
121	IF(Z.GT.0.0) T=0.0
122	IF(Z.LT.0.0) T=180.0
123	GO TO 40
124	C
125	30 CONTINUE
126	T=90.0
127	IF(ABS(Z).LT.1.0E-6) GO TO 40
128	C
129	R=SQRT(R)
130	T=ATAN(R/Z)*RADDEG
131	IF(T.EQ.0.0.AND.Z.LT.0.0) T=180.0
132	IF(T.LT.0.0) T=T+180.0
133	C
134	40 CONTINUE
135	IF(T.GT.TH) TH=T
136	IF(T.LT.TL) TL=T
137	C
138	50 CONTINUE
139	C
140	C THETA LIMITS SET FROM THE POINTS NOW DO THE EDGES.
141	C
142	DO 120 IL=1,12
143	C
144	C FIND THE END POINT NUMBERS FOR EACH EDGE.
145	C
146	IF(IL.GT.4) GO TO 60
147	IPP1=1
148	IF(IL.GT.2) IPP1=4
149	IPP2=2
150	IF(MOD(IL,2).EQ.1) IPP2=3
151	C
152	GO TO 80
153	60 CONTINUE
154	IF(IL.LT.9) GO TO 70
155	C
156	IPP1=5
157	IF(IL.GT.10) IPP1=8
158	IPP2=6
159	IF(MOD(IL,2).EQ.1) IPP2=7
160	C
161	GO TO 80
162	70 CONTINUE
163	C
164	IPP1=IL-4
165	IPP2=IL
166	C
167	80 CONTINUE
168	C
169	C NOW GET THE POINTS AND ROTATE THEM.
170	C
171	IF(ISH.EQ.28) THEN
172	C
173	C General twisted trapezoid.
174	C
175	ILL=IPP1
176	IF(IPP1.EQ.3) ILL=4
177	IF(IPP1.EQ.4) ILL=3
178	I0=ILL*4+11
179	X1(3)=PARS(1)
180	IF(IPP1.LT.5) X1(3)=-X1(3)
181	X1(1)=PARS(I0)+PARS(I0+2)*X1(3)
182	X1(2)=PARS(I0+1)+PARS(I0+3)*X1(3)
183	ILL=IPP2
184	IF(IPP2.EQ.3) ILL=4
185	IF(IPP2.EQ.4) ILL=3
186	I0=ILL*4+11
187	X2(3)=PARS(1)
188	IF(IPP2.LT.5) X2(3)=-X2(3)
189	X2(1)=PARS(I0)+PARS(I0+2)*X2(3)
190	X2(2)=PARS(I0+1)+PARS(I0+3)*X2(3)
191	C
192	GO TO 100
193	C
194	ENDIF
195	C
196	IP3=ISH+2
197	IF(ISH.EQ.10) IP3=3
198	IF(ISH.EQ.4) IP3=1
199	X1(3)=PARS(IP3)
200	IF(IPP1.LE.4) X1(3)=-X1(3)
201	IP2=3
202	IF(ISH.GT.2.AND.X1(3).GT.0.0) IP2=4
203	IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
204	IF(ISH.EQ.4) IP2=4
205	IF(ISH.EQ.4.AND.X1(3).GT.0.0) IP2=8
206	X1(2)=PARS(IP2)
207	IF(MOD(IPP1+3,4).LT.2) X1(2)=-X1(2)
208	IP1=1
209	IF(ISH.NE.1.AND.ISH.NE.10.AND.X1(3).GT.0.0) IP1=2
210	IF(ISH.EQ.4) IP1=5
211	IF(ISH.EQ.4.AND.X1(3).GT.0.0) IP1=IP1+4
212	IF(ISH.EQ.4.AND.X1(2).GT.0.0) IP1=IP1+1
213	X1(1)=PARS(IP1)
214	IF(MOD(IPP1,2).EQ.1) X1(1)=-X1(1)
215	C
216	IP3=ISH+2
217	IF(ISH.EQ.10) IP3=3
218	IF(ISH.EQ.4) IP3=1
219	X2(3)=PARS(IP3)
220	IF(IPP2.LE.4) X2(3)=-X2(3)
221	IP2=3
222	IF(ISH.GT.2.AND.X2(3).GT.0.0) IP2=4
223	IF(ISH.EQ.1.OR.ISH.EQ.10) IP2=2
224	IF(ISH.EQ.4) IP2=4
225	IF(ISH.EQ.4.AND.X2(3).GT.0.0) IP2=8
226	X2(2)=PARS(IP2)
227	IF(MOD(IPP2+3,4).LT.2) X2(2)=-X2(2)
228	IP1=1
229	IF(ISH.NE.1.AND.ISH.NE.10.AND.X2(3).GT.0.0) IP1=2
230	IF(ISH.EQ.4) IP1=5
231	IF(ISH.EQ.4.AND.X2(3).GT.0.0) IP1=IP1+4
232	IF(ISH.EQ.4.AND.X2(2).GT.0.0) IP1=IP1+1
233	X2(1)=PARS(IP1)
234	IF(MOD(IPP2,2).EQ.1) X2(1)=-X2(1)
235	C
236	IF(ISH.NE.10) GO TO 90
237	X1(1)=X1(1)+X1(2)PARS(4)+X1(3)PARS(5)
238	X1(2)=X1(2)+X1(3)*PARS(6)
239	X2(1)=X2(1)+X2(2)PARS(4)+X2(3)PARS(5)
240	X2(2)=X2(2)+X2(3)*PARS(6)
241	90 CONTINUE
242	C
243	IF(ISH.NE.4) GO TO 100
244	IP4=7
245	IF(X1(3).GT.0.0) IP4=11
246	X1(1)=X1(1)+X1(2)PARS(IP4)+X1(3)PARS(2)
247	X1(2)=X1(2)+X1(3)*PARS(3)
248	IP4=7
249	IF(X2(3).GT.0.0) IP4=11
250	X2(1)=X2(1)+X2(2)PARS(IP4)+X2(3)PARS(2)
251	X2(2)=X2(2)+X2(3)*PARS(3)
252	100 CONTINUE
253	C
254	C ROTATE.
255	C
256	JROT=LQ(JROTM-IROT)
257	XT1(1)=X1(1)
258	XT1(2)=X1(2)
259	XT1(3)=X1(3)
260	IF(IROT.NE.0) CALL GINROT(X1,Q(JROT+1),XT1)
261	XT2(1)=X2(1)
262	XT2(2)=X2(2)
263	XT2(3)=X2(3)
264	IF(IROT.NE.0) CALL GINROT(X2,Q(JROT+1),XT2)
265	C
266	C NOW WE HAVE THE END POINTS IN THE MASTER SYSTEM.
267	C FIND THE TANGENT POINT TO THE SET OF CONES OF CONSTANT
268	C THETA.
269	C
270	DS=(XT2(1)-XT1(1))2+(XT2(2)-XT1(2))2+(XT2(3)-XT1(3))**2
271	IF(DS.LE.0.0) GO TO 120
272	C
273	DS=SQRT(DS)
274	C
275	X0=(XT2(1)+XT1(1))/2.0+DX(1)
276	Y0=(XT2(2)+XT1(2))/2.0+DX(2)
277	Z0=(XT2(3)+XT1(3))/2.0+DX(3)
278	ALX=(XT2(1)-XT1(1))/DS
279	ALY=(XT2(2)-XT1(2))/DS
280	ALZ=(XT2(3)-XT1(3))/DS
281	C
282	T=90.0
283	IF(Z0.EQ.0.0.AND.ALZ.EQ.0.0) GO TO 110
284	C
285	IF(ALX.EQ.0.0.AND.ALY.EQ.0.0) GO TO 120
286	C THIS CHECKS WHETHER THE LINE IS PARALLEL TO THE
287	C Z AXIS IN WHICH CASE THERE IS NO TANGENT AND
288	C THE END POINTS DETERMINE THE THETA RANGE.
289	C
290	SNUM=(X0Z0ALX+Y0Z0ALY+X0X0ALZ+Y0Y0ALZ)
291	SDEN=(Z0ALXALX-X0ALXALZ+Z0ALYALY-Y0ALYALZ)
292	C
293	IF(ABS(SNUM).GT.0.5DSABS(SDEN)) GO TO 120
294	C
295	C TANGENT EXIST BETWEEN THE TWO ENDS.
296	C
297	S = SNUM/SDEN
298	X0=X0+S*ALX
299	Y0=Y0+S*ALY
300	Z0=Z0+S*ALZ
301	C
302	R=X0X0+Y0Y0
303	RT=R+Z0*Z0
304	C
305	IF(RT.LT.1.0E-10) GO TO 999
306	C
307	IF(R.GT.0.0) R=SQRT(R)
308	T=90.0
309	IF(ABS(Z0).LT.1.0E-06) GO TO 110
310	C
311	T=ATAN(R/Z0)*RADDEG
312	IF(T.EQ.0.0.AND.Z0.LT.0.0) T=180.0
313	IF(T.LT.0.0) T=T+180.0
314	C
315	110 CONTINUE
316	C
317	IF(T.LT.TL) TL=T
318	IF(T.GT.TH) TH=T
319	C
320	C CHECK FOR THE POSSIBILITY OF STRADDLING T=0.0 &/OR 180.0.
321	C
322	C=X0DX(1)+Y0DX(2)
323	IF(C.GT.0.0) GO TO 120
324	C
325	C CHECK IF SAME SIGN OF Z.
326	C
327	IF(Z0*DX(3).LT.0.0) GO TO 999
328	C
329	T=0.0
330	IF(Z0.LT.0.0) T=180.0
331	IF(T.LT.TL) TL=T
332	IF(T.GT.TH) TH=T
333	C
334	120 CONTINUE
335	C
336	C DONE SET THE LIMITS.
337	C
338	CL=TL
339	CH=TH
340	C
341	GO TO 999
342	C
343	130 CONTINUE
344	C
345	C TUBES, CONES ETC.
346	C
347	IF(IROT.NE.0.OR.DXS.GT.1.0E-05) GO TO 180
348	C
349	C UNROTATED AND CENTERED.
350	C
351	IF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13
352	+ .AND.ISH.NE.14)GO TO 170
353	C
354	C TUBES AND CONES.
355	C
356	IERR=0
357	DZ=PARS(3)
358	RMN=PARS(1)
359	RMX=PARS(2)
360	IF(ISH.EQ.13) THEN
361	**
362	** approxime to a cylinder whit radius
363	** equal to the ellipse major axis
364	**
365	RMN=0.0
366	IF(PARS(1).GT.RMX) RMX=PARS(1)
367	GOTO 140
368	ENDIF
369	IF(ISH.EQ.14) THEN
370	C not really sure of the function of these... keep RMN=PARS(1)
371	RMX = SQRT(PARS(2)*2+(PARS(3)TAN(PARS(4)DEGRAD))*2)
372	GO TO 140
373	ENDIF
374	IF(ISH.EQ.NSCTUB) THEN
375	S1 = (1.0-PARS(8))*(1.0+PARS(8))
376	IF( S1 .GT. 0.0) S1 = SQRT(S1)
377	S2 = (1.0-PARS(11))*(1.0+PARS(11))
378	IF( S2 .GT. 0.0) S2 = SQRT(S2)
379	IF( S2 .GT. S1 ) S1 = S2
380	DZ = DZ+RMX*S1
381	ENDIF
382	IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 140
383	DZ=PARS(1)
384	RMN=PARS(2)
385	RMX=PARS(3)
386	IF(PARS(4).LT.RMN) RMN=PARS(4)
387	IF(PARS(5).GT.RMX) RMX=PARS(5)
388	140 CONTINUE
389	C
390	IF(DZ.GT.ABS(DX(3))) GO TO 160
391	C
392	C ALL FORWARD OR ALL BACK.
393	C
394	DSH=DX(3)-DZ
395	DLN=DX(3)+DZ
396	IF(DX(3).GT.0.0) GO TO 150
397	DSS=DSH
398	DSH=DLN
399	DLN=DSS
400	150 CONTINUE
401	C
402	CL=90.0
403	CH=90.0
404	IF(DLN.NE.0.0) CL=ATAN(RMN/DLN)*RADDEG
405	IF(DSH.NE.0.0) CH=ATAN(RMX/DSH)*RADDEG
406	IF(DX(3).GT.0.0) GO TO 999
407	CS=CL
408	CL=CH
409	CH=CS
410	IF(CH.EQ.0.0) CH=180.0
411	IF(CL.LT.0.0) CL=CL+180.0
412	IF(CH.LT.0.0) CH=CH+180.0
413	C
414	GO TO 999
415	C
416	160 CONTINUE
417	C
418	CL=90.0
419	CH=90.0
420	IF(DZ+DX(3).NE.0.0) CL=ATAN(RMN/(DZ+DX(3)))*RADDEG
421	IF(-DZ+DX(3).NE.0.0) CH=ATAN(RMN/(-DZ+DX(3)))*RADDEG
422	IF(CH.LE.0.0) CH=CH+180.0
423	C
424	GO TO 999
425	C
426	170 CONTINUE
427	IF(ISH.GT.9) GO TO 999
428	C
429	C SPHERE.
430	C
431	IERR=0
432	CL=PARS(3)
433	CH=PARS(4)
434	C
435	GO TO 999
436	180 CONTINUE
437	C
438	IF(ISH.EQ.11.OR.ISH.EQ.12) GOTO 999
439	**
440	RM=PARS(2)
441	IF(ISH.EQ.9) GO TO 200
442	C
443	DZ=PARS(3)
444	IF(ISH.EQ.NSCTUB) THEN
445	S1 = (1.0-PARS(8))*(1.0+PARS(8))
446	IF( S1 .GT. 0.0) S1 = SQRT(S1)
447	S2 = (1.0-PARS(11))*(1.0+PARS(11))
448	IF( S2 .GT. 0.0) S2 = SQRT(S2)
449	IF( S2 .GT. S1 ) S1 = S2
450	DZ = DZ+RM*S1
451	ENDIF
452	**
453	IF(ISH.EQ.13) THEN
454	IF(PARS(1).GT.RM) RM=PARS(1)
455	GOTO 190
456	ENDIF
457	**
458	IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 190
459	C
460	DZ=PARS(1)
461	RM=PARS(3)
462	IF(PARS(5).GT.RM) RM=PARS(5)
463	C
464	190 CONTINUE
465	C
466	RM=SQRT(RM2+DZ2)
467	C
468	200 CONTINUE
469	C
470	CL=0.0
471	CH=180.0
472	IERR=0
473	RC=DXS2+DX(3)2
474	IF(RC.GT.0.0) RC=SQRT(RC)
475	IF(RM.GE.RC) GO TO 999
476	C
477	TC=90.0
478	IF(ABS(DX(3)).GT.0.0) TC=ATAN(DXS/DX(3))*RADDEG
479	IF(TC.LT.0.0) TC=TC+180.0
480	C
481	DT=ASIN(RM/RC)*RADDEG
482	CL=TC-DT
483	CH=TC+DT
484	IF(CL.LT.0.0) CL=0.0
485	IF(CH.GT.180.0) CH=180.0
486	C
487	999 CONTINUE
488	END