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

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1999/05/18 15:55:17  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:20:56  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 21/07/94  11.48.23  by  S.Giani
*-- Author :
      SUBROUTINE GVDRAD(IAXIS,ISH,IROT,DX,PARS,CL,CH,IERR)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *    ROUTINE TO COMPUTE THE LIMITS IN R FOR THE SHAPE ISH        *
C.    *    DISPLACED BY THE VECTOR DX AND ROTATED BY THE MATRIX IROT.  *
C.    *    IF IAXIS = 4 THE R IS THE XY PLANE R, IF IAXIS = 5 IT IS    *
C.    *    THE 3 DINEMSIONAL SPACE R. THE SHAPE HAS NPAR PARAMETERS    *
C.    *    IN THE ARRAY PARS. THE LOWER LIMIT IS RETURNED IN CL AND    *
C.    *    THE HIGHER IN CH. IF THE CALCULATION CANNOT BE PERFORMED    *
C.    *    IERR IS SET TO 1 OTHERWISE IT IS SET TO 0.                  *
C.    *                                                                *
C.    *    ==>Called by : GVDLIM                                       *
C.    *         Author  S.Giani  *********                             *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gcbank.inc"
#include "geant321/gconsp.inc"
#include "geant321/gcshno.inc"
      DIMENSION DX(3),PARS(100),X(3),XT(3)
C.
C.           --------------------------------------------------
C.
      IERR=1
C
C            FIRST CALCULATE THE LENGTH OF THE DISPLACEMENT OF THE
C            ORIGIN.
C
      DXS=DX(1)*DX(1)+DX(2)*DX(2)
      IF(IAXIS.EQ.5) DXS=DXS+DX(3)*DX(3)
      IF(DXS.GT.0.0) DXS=SQRT(DXS)
C
      IF(ISH.GT.4.AND.ISH.NE.10.AND.ISH.NE.28) GO TO 40
C
C          CUBOIDS, TRAPEZOIDS, PARALLELEPIPEDS.
C
      CH=0.0
      CL=DXS
C
      DO 30 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
C          NOW COMPUTE RMIN = PROJECTED R ON DX AND RMAX = R
C          AND UPDATE LIMITS IF NECESSARY.
C
      R2=(XT(1)+DX(1))**2+(XT(2)+DX(2))**2
      IF(IAXIS.EQ.5) R2=R2+(XT(3)+DX(3))**2
      R=SQRT(R2)
      IF(R.GT.CH) CH=R
C
      IF(CL.LE.0.0) GO TO 30
C
      XPT=DX(1)*XT(1)+DX(2)*XT(2)
      IF(IAXIS.EQ.5) XPT=XPT+DX(3)*XT(3)
      IF(DXS.LE.1.0E-05) GO TO 30
      RMN=DXS+XPT/DXS
      IF(RMN.LT.CL) CL=RMN
C
   30 CONTINUE
C
      IF(CL.LE.0.0) CL=0.0
C
      IERR=0
      GO TO 999
C
   40 CONTINUE
C
C          POLYGONES AND POLYCONES
C
      IF(ISH.EQ.11.AND.IAXIS.EQ.4)THEN
       NZLAST=PARS(4)
       IZLAST=2+3*NZLAST
       CLZ=PARS(5)
       CHZ=PARS(IZLAST)
       DZ2=ABS(CHZ-CLZ)
       DZ=DZ2*.5
       TMPRAD=0.
       TMPMIN=100000.
       DO 145 I=7,IZLAST+2,3
         IF(PARS(I).GT.TMPRAD)TMPRAD=PARS(I)
         IF(PARS(I-1).LT.TMPMIN)TMPMIN=PARS(I-1)
 145   CONTINUE
       RMN=TMPMIN
       PHIMIN=PARS(1)
       PHIMAX=PHIMIN+PARS(2)
       AANG=ABS(PHIMAX-PHIMIN)
       NANG=PARS(3)
       AATMAX=NANG*360./AANG
       LATMAX=AATMAX
       ALA=AATMAX-LATMAX
       IF(ALA.GT..5)LATMAX=LATMAX+1
       AFINV=1./COS(PI/LATMAX)
       FINV=ABS(AFINV)
       R=TMPRAD*FINV
       GOTO 50
      ELSEIF(ISH.EQ.12.AND.IAXIS.EQ.4)THEN
       NZLAST=PARS(3)
       IZLAST=1+3*NZLAST
       CLZ=PARS(4)
       CHZ=PARS(IZLAST)
       DZ2=ABS(CHZ-CLZ)
       DZ=DZ2*.5
       TMPRAD=0.
       TMPMIN=100000.
       DO 146 I=6,IZLAST+2,3
         IF(PARS(I).GT.TMPRAD)TMPRAD=PARS(I)
         IF(PARS(I-1).LT.TMPMIN)TMPMIN=PARS(I-1)
 146   CONTINUE
       RMN=TMPMIN
       R=TMPRAD
       GOTO 50
      ENDIF
      IF(ISH.GT.8.AND.ISH.NE.NSCTUB.AND.ISH.NE.13.AND.ISH.NE.14)GO TO 80
C
C             TUBES AND CONES
C
      IP3=3
      IF(ISH.GT.6.AND.ISH.NE.NSCTUB.AND.ISH.NE.13.AND.ISH.NE.14) IP3=1
      DZ=PARS(IP3)
      R=PARS(2)
      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+R*S1
      ENDIF
**
      IF(ISH.EQ.13) THEN
**
**       APPROXIME TO A CYLINDER WHIT RADIUS
**       EQUAL TO THE ELLIPSE MAJOR AXIS
**
         RMN=0.0
         IF(PARS(1).GT.R) R=PARS(1)
         GOTO 50
      ENDIF
      RMN=PARS(1)
*
      IF(ISH.EQ.14) THEN
        R = SQRT(PARS(2)**2+(PARS(3)*TAN(PARS(4)*DEGRAD))**2)
        GO TO 50
      ENDIF
C
      IF(ISH.LE.6.OR.ISH.EQ.NSCTUB) GO TO 50
C
      R=PARS(3)
      IF(PARS(5).GT.R) R=PARS(5)
      RMN=PARS(2)
      IF(PARS(4).LT.RMN) RMN=PARS(4)
C
   50 CONTINUE
C
C          ROTATE THE LOCAL Z AXIS.
C
      X(1)=0.0
      X(2)=0.0
      X(3)=1.0
      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
C          COMPUTE RMIN AND RMAX ASSUMING COMPLETE TUBE HALF
C          LENGTH DZ AND RADIUS R.
C
      ST2=1.0
      IF(IAXIS.EQ.4) ST2=(1+XT(3))*(1-XT(3))
      DR=SQRT(DZ*DZ*ST2+R*R)
      CL=DXS-DR
      IF(CL.LT.0.0) CL=0.0
      CH=DXS+DR
      IF(ABS(XT(3)).EQ.1.0.AND.DXS.LT.1.0E-05) CL=RMN
      IERR=0
C
      GO TO 999
C
   80 CONTINUE
      IF(ISH.GT.9) GO TO 999
C
C           SPHERE.
C
      CL=DXS-PARS(2)
      IF(CL.LT.0.0) CL=0.0
      CH=DXS+PARS(2)
      IF(IAXIS.EQ.5.AND.DXS.LT.1.0E-05) CL=PARS(1)
      IERR=0
C
  999 CONTINUE
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
d43b40e2	5	* Revision 1.1.1.1 1999/05/18 15:55:17 fca
	6	* AliRoot sources
	7	*
fe4da5cc	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"
d43b40e2	35	DIMENSION DX(3),PARS(100),X(3),XT(3)
fe4da5cc	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-IL2)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(DZDZST2+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