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

*
* $Id$
*
* $Log$
* Revision 1.2  1996/02/27 09:54:45  ravndal
* Double precision for some locals, thanks to G.Poulard
*
* Revision 1.1.1.1  1995/10/24 10:20:56  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.30  by  S.Giani
*-- Author :
      SUBROUTINE GSNGTR(X,P,IACT,SNEXT,SNXT,SAFE,INSIDE)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *                                                                *
C.    *    Routine to determine the shortest distance from the point   *
C.    *    X(1-3) to the boundary of the shape of type GTRA defined    *
C.    *    by the parameters P along the vector X(4-6). If INSIDE is   *
C.    *    1 then the point is inside the shape and the distance is    *
C.    *    returned as SNEXT. If INSIDE is 0 then the point is         *
C.    *    outside the shape and if the line hits the shape then       *
C.    *    if the new distance is less than the                        *
C.    *    old value of SNEXT the new distance is returned as SNEXT.   *
C.    *                                                                *
C.    *          Called by : GNEXT, GTNEXT                             *
C.    *          A.C.McPherson   22nd April 1985.                      *
C.    *                                                                *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gconsp.inc"
#include "geant321/gcunit.inc"
C
#if !defined(CERNLIB_SINGLE)
      DOUBLE PRECISION X0,Y0,DXDZ,DYDZ,A,B,C,DISC,X1,X2,X3,SN,CP,SMALL
#endif
      PARAMETER (SMALL=1E-10)
C
      DIMENSION X(6),P(30),SN(2,5),IOUT(5),X0(4),Y0(4),DXDZ(4),DYDZ(4)
C.
C.                ---------------------------------------------
C.
C
C               Compute Safety distance
C
      DOUBLE PRECISION SL,SL1,SM,SM1
      IF(IACT.LT.3) CALL GSAGTR(X,P,SAFE,INSIDE)
      SNXT=BIG
      IF (IACT .EQ. 0) GO TO 999
      IF (IACT .EQ. 1) THEN
        IF (SNEXT .LT. SAFE) GO TO 999
      ENDIF
C
C               First compute the distance along the line to the
C               boundaries.
C
C               The distance to the planes defined by z=+/-P(1).
C
      IF(X(6).EQ.0.0) THEN
          SN(1,1)=BIG
          SN(2,1)=BIG
          GOTO 10
      ENDIF
      SN(1,1)=(-P(1)-X(3))/X(6)
      SN(2,1)=(P(1)-X(3))/X(6)
      IF(X(6).GT.0.0) GO TO 10
      ST=SN(2,1)
      SN(2,1)=SN(1,1)
      SN(1,1)=ST
   10 CONTINUE
C
C               The distance to the remaining four surfaces.
C
      DO 20 I=1,4
      X0(I)=P(I*4+11)
      Y0(I)=P(I*4+12)
      DXDZ(I)=P(I*4+13)
      DYDZ(I)=P(I*4+14)
   20 CONTINUE
C
      DO 65 I=1,4
      J=I+1
      IF(J.EQ.5) J=1
C
      A=(X(4)-DXDZ(I)*X(6))*(DYDZ(J)-DYDZ(I))*X(6) -
     +(X(5)-DYDZ(I)*X(6))*(DXDZ(J)-DXDZ(I))*X(6)
C
      B=(X(4)-DXDZ(I)*X(6))*(Y0(J)-Y0(I)+(DYDZ(J)-DYDZ(I))*X(3)) +
     +(X(1)-X0(I)-DXDZ(I)*X(3))*(DYDZ(J)-DYDZ(I))*X(6) -
     +(X(5)-DYDZ(I)*X(6))*(X0(J)-X0(I)+(DXDZ(J)-DXDZ(I))*X(3)) -
     +(X(2)-Y0(I)-DYDZ(I)*X(3))*(DXDZ(J)-DXDZ(I))*X(6)
C
      C=(X(1)-X0(I)-DXDZ(I)*X(3))*(Y0(J)-Y0(I)+(DYDZ(J)-DYDZ(I))*X(3))
     + - (X(2)-Y0(I)-DYDZ(I)*X(3))*(X0(J)-X0(I)+(DXDZ(J)-DXDZ(I))*X(3))
C
      IOUT(I+1)=0
      IF(C.GT.0.0) IOUT(I+1)=1
C
C             The solutions are in the normal form:
C             s = (-B+/-SQRT(B*B-4.0*A*C))*0.5/A
C
      SN(1,I+1)=BIG
      SN(2,I+1)=BIG
      IF(ABS(A).GT.1.0E-10) GO TO 30
C
C             A = 0 only one solution.
C
      IF(ABS(B).LT.1.0E-10) GO TO 60
C
      SN(1,I+1)=-C/B
      GO TO 60
C
   30 CONTINUE
      IF(ABS(C).GT.1.0E-10) GO TO 40
      SN(1,I+1)=0.0
      SN(2,I+1)=0.0
      IF(ABS(B).LT.1.0E-10) GO TO 60
      SN(1,I+1)=-C/B
      IF(C.EQ.0.0) SN(1,I+1)=SIGN(SMALL,B)
      SN(2,I+1)=-B/A
      GO TO 50
C
   40 CONTINUE
      DISC=B*B-A*C*4.0
      IF(DISC.LT.0.0) GO TO 60
      IF(DISC.GT.0.0) DISC=SQRT(DISC)
      SN(1,I+1)=(-B-DISC)*0.5/A
      SN(2,I+1)=(-B+DISC)*0.5/A
C
   50 CONTINUE
      IF(SN(2,I+1).GT.SN(1,I+1)) GO TO 60
      ST=SN(2,I+1)
      SN(2,I+1)=SN(1,I+1)
      SN(1,I+1)=ST
C
   60 CONTINUE
C
      DO 65 K=1,2
      IF(ABS(SN(K,I+1)).GT.1.0E+05.OR.ABS(SN(K,I+1)).LT.1.0E-05)
     +GO TO 65
C
      X1=X(1)+SN(K,I+1)*X(4)
      X2=X(2)+SN(K,I+1)*X(5)
      X3=X(3)+SN(K,I+1)*X(6)
      CP=(X1-X0(I)-DXDZ(I)*X3)*(Y0(J)-Y0(I)+(DYDZ(J)-DYDZ(I))*X3)
     + - (X2-Y0(I)-DYDZ(I)*X3)*(X0(J)-X0(I)+(DXDZ(J)-DXDZ(I))*X3)
      CP=CP/SQRT((X0(J)-X0(I)+(DXDZ(J)-DXDZ(I))*X3)**2+
     +   (Y0(J)-Y0(I)+(DYDZ(J)-DYDZ(I))*X3)**2)
C
      IF(ABS(CP).LT.0.0001) GO TO 65
      IF(ABS(CP/SN(K,I+1)).LT.1.0E-06) GO TO 65
      WRITE(CHMAIL,1020) I,K,SN(K,I+1)
      CALL GMAIL(0,0)
      WRITE(CHMAIL,1021) X
      CALL GMAIL(0,0)
      WRITE(CHMAIL,1022) X1,X2,X3,CP
      CALL GMAIL(0,0)
      WRITE(CHMAIL,1023) A,B,C,DISC
      CALL GMAIL(0,0)
      WRITE(CHMAIL,1024) INSIDE
      CALL GMAIL(0,0)
      WRITE(CHMAIL,1025) X0
      CALL GMAIL(0,0)
      WRITE(CHMAIL,1026) Y0
      CALL GMAIL(0,0)
      WRITE(CHMAIL,1027) DXDZ
      CALL GMAIL(0,0)
      WRITE(CHMAIL,1028) DYDZ
      CALL GMAIL(0,0)
 1020 FORMAT('0 GSNGTR ERROR - I,K =',2I2,' SN =',E13.5)
 1021 FORMAT(' X =',6F15.6)
 1022 FORMAT(' X1,X2,X3 =',3F15.6,' CP =',E15.6)
 1023 FORMAT(' A =',E15.6,' B =',E15.6,' C =',E15.6,' DISC =',E15.6)
 1024 FORMAT(' INSIDE =',I3)
 1025 FORMAT('   X0 =',4E15.6)
 1026 FORMAT('   Y0 =',4E15.6)
 1027 FORMAT(' DXDZ =',4E15.6)
 1028 FORMAT(' DYDZ =',4E15.6)
C
   65 CONTINUE
C
C
C             Have computed the two distances for the z planes and
C             the four surfaces. Combine them accordingly as to
C             whether the point is inside or outside the shape.
C
      IF(INSIDE.EQ.0) GO TO 80
C
C             Point is inside shape.
C
      DO 70 I=1,5
      DO 70 J=1,2
      IF(SN(J,I).GT.0.0.AND.SN(J,I).LT.SNXT) SNXT=SN(J,I)
   70 CONTINUE
      GO TO 999
C
   80 CONTINUE
C
C             Point is outside shape.
C
      IOUT(1)=0
      IF(ABS(X(3)).GT.P(1)) IOUT(1)=1
C
C             For each of five sets of SN and IOUT, IOUT(I) equal to 1
C             indicates that the point is outside the shape by the Ith
C             test, SN(1,I) is the distance to the first change in the
C             test and SN(2,I) is the distance to the second change.
C             The remaining logic just attempts to find a distance when
C             the line is inside by all five tests, bearing in mind that
C             for some tests the line can start inside, leave and return
C             inside.
C
      SL=-1.0
      SM=BIG
      SM1=BIG
      DO 100 I=1,5
      IF(IOUT(I).EQ.0) GO TO 90
      IF(SN(2,I).LT.0.0) GO TO 999
      IF(SN(1,I).LT.0.0.AND.SN(2,I).GT.SL) SL=SN(2,I)
      IF(SN(1,I).GT.SL) SL=SN(1,I)
      IF(SN(1,I).GE.0.0.AND.SN(2,I).LT.SM) SM=SN(2,I)
      GO TO 100
   90 CONTINUE
      IF(SN(1,I).LT.0.0.AND.SN(2,I).GE.0.0.AND.SN(2,I).LT.SM) SM=SN(2,I)
      IF(SN(1,I).LT.0.0.OR.SN(1,I).GT.SM1) GO TO 100
      IF(SN(1,I).GE.SN(2,I)) GO TO 100
      SM1=SN(1,I)
      SL1=SN(2,I)
  100 CONTINUE
C
C             SL is the largest of the five distances to the first
C             time the line is inside. SM is the smallest to the
C             last time the point is inside. SM1 is the smallest
C             distance to when the line is temporarily outside
C             one of the tests.
C
      IF(SM.LE.SL) GO TO 999
      IF(SM1.GT.SL) GO TO 130
C
  110 CONTINUE
C
C             In this loop SL is updated by the return after SM1
C             if SM1 is less than SL.
C
      SL=SL1
      IF(SM.LE.SL) GO TO 999
      SM1=SM
C
      DO 120 I=1,5
      IF(IOUT(I).EQ.1) GO TO 120
      IF(SN(2,I).LE.SL.OR.SN(1,I).GT.SM1) GO TO 120
      IF(SN(1,I).GE.SN(2,I)) GO TO 120
      SM1=SN(1,I)
      SL1=SN(2,I)
  120 CONTINUE
C
      IF(SM1.GT.SL) GO TO 130
C
      GO TO 110
  130 CONTINUE
C
      IF(SL.LT.SNXT) SNXT=SL
C
  999 CONTINUE
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.2 1996/02/27 09:54:45 ravndal
	6	* Double precision for some locals, thanks to G.Poulard
	7	*
	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 29/03/94 15.41.30 by S.Giani
	14	*-- Author :
	15	SUBROUTINE GSNGTR(X,P,IACT,SNEXT,SNXT,SAFE,INSIDE)
	16	C.
	17	C. ******************************************************************
	18	C. * *
	19	C. * *
	20	C. * Routine to determine the shortest distance from the point *
	21	C. * X(1-3) to the boundary of the shape of type GTRA defined *
	22	C. * by the parameters P along the vector X(4-6). If INSIDE is *
	23	C. * 1 then the point is inside the shape and the distance is *
	24	C. * returned as SNEXT. If INSIDE is 0 then the point is *
	25	C. * outside the shape and if the line hits the shape then *
	26	C. * if the new distance is less than the *
	27	C. * old value of SNEXT the new distance is returned as SNEXT. *
	28	C. * *
	29	C. * Called by : GNEXT, GTNEXT *
	30	C. * A.C.McPherson 22nd April 1985. *
	31	C. * *
	32	C. * *
	33	C. ******************************************************************
	34	C.
	35	#include "geant321/gconsp.inc"
	36	#include "geant321/gcunit.inc"
	37	C
	38	#if !defined(CERNLIB_SINGLE)
	39	DOUBLE PRECISION X0,Y0,DXDZ,DYDZ,A,B,C,DISC,X1,X2,X3,SN,CP,SMALL
	40	#endif
	41	PARAMETER (SMALL=1E-10)
	42	C
	43	DIMENSION X(6),P(30),SN(2,5),IOUT(5),X0(4),Y0(4),DXDZ(4),DYDZ(4)
	44	C.
	45	C. ---------------------------------------------
	46	C.
	47	C
	48	C Compute Safety distance
	49	C
	50	DOUBLE PRECISION SL,SL1,SM,SM1
	51	IF(IACT.LT.3) CALL GSAGTR(X,P,SAFE,INSIDE)
	52	SNXT=BIG
	53	IF (IACT .EQ. 0) GO TO 999
	54	IF (IACT .EQ. 1) THEN
	55	IF (SNEXT .LT. SAFE) GO TO 999
	56	ENDIF
	57	C
	58	C First compute the distance along the line to the
	59	C boundaries.
	60	C
	61	C The distance to the planes defined by z=+/-P(1).
	62	C
	63	IF(X(6).EQ.0.0) THEN
	64	SN(1,1)=BIG
65	SN(2,1)=BIG
66	GOTO 10
67	ENDIF
68	SN(1,1)=(-P(1)-X(3))/X(6)
69	SN(2,1)=(P(1)-X(3))/X(6)
70	IF(X(6).GT.0.0) GO TO 10
71	ST=SN(2,1)
72	SN(2,1)=SN(1,1)
73	SN(1,1)=ST
74	10 CONTINUE
75	C
76	C The distance to the remaining four surfaces.
77	C
78	DO 20 I=1,4
79	X0(I)=P(I*4+11)
80	Y0(I)=P(I*4+12)
81	DXDZ(I)=P(I*4+13)
82	DYDZ(I)=P(I*4+14)
83	20 CONTINUE
84	C
85	DO 65 I=1,4
86	J=I+1
87	IF(J.EQ.5) J=1
88	C
89	A=(X(4)-DXDZ(I)X(6))(DYDZ(J)-DYDZ(I))*X(6) -
90	+(X(5)-DYDZ(I)X(6))(DXDZ(J)-DXDZ(I))*X(6)
91	C
92	B=(X(4)-DXDZ(I)X(6))(Y0(J)-Y0(I)+(DYDZ(J)-DYDZ(I))*X(3)) +
93	+(X(1)-X0(I)-DXDZ(I)X(3))(DYDZ(J)-DYDZ(I))*X(6) -
94	+(X(5)-DYDZ(I)X(6))(X0(J)-X0(I)+(DXDZ(J)-DXDZ(I))*X(3)) -
95	+(X(2)-Y0(I)-DYDZ(I)X(3))(DXDZ(J)-DXDZ(I))*X(6)
96	C
97	C=(X(1)-X0(I)-DXDZ(I)X(3))(Y0(J)-Y0(I)+(DYDZ(J)-DYDZ(I))*X(3))
98	+ - (X(2)-Y0(I)-DYDZ(I)X(3))(X0(J)-X0(I)+(DXDZ(J)-DXDZ(I))*X(3))
99	C
100	IOUT(I+1)=0
101	IF(C.GT.0.0) IOUT(I+1)=1
102	C
103	C The solutions are in the normal form:
104	C s = (-B+/-SQRT(BB-4.0AC))0.5/A
105	C
106	SN(1,I+1)=BIG
107	SN(2,I+1)=BIG
108	IF(ABS(A).GT.1.0E-10) GO TO 30
109	C
110	C A = 0 only one solution.
111	C
112	IF(ABS(B).LT.1.0E-10) GO TO 60
113	C
114	SN(1,I+1)=-C/B
115	GO TO 60
116	C
117	30 CONTINUE
118	IF(ABS(C).GT.1.0E-10) GO TO 40
119	SN(1,I+1)=0.0
120	SN(2,I+1)=0.0
121	IF(ABS(B).LT.1.0E-10) GO TO 60
122	SN(1,I+1)=-C/B
123	IF(C.EQ.0.0) SN(1,I+1)=SIGN(SMALL,B)
124	SN(2,I+1)=-B/A
125	GO TO 50
126	C
127	40 CONTINUE
128	DISC=BB-AC*4.0
129	IF(DISC.LT.0.0) GO TO 60
130	IF(DISC.GT.0.0) DISC=SQRT(DISC)
131	SN(1,I+1)=(-B-DISC)*0.5/A
132	SN(2,I+1)=(-B+DISC)*0.5/A
133	C
134	50 CONTINUE
135	IF(SN(2,I+1).GT.SN(1,I+1)) GO TO 60
136	ST=SN(2,I+1)
137	SN(2,I+1)=SN(1,I+1)
138	SN(1,I+1)=ST
139	C
140	60 CONTINUE
141	C
142	DO 65 K=1,2
143	IF(ABS(SN(K,I+1)).GT.1.0E+05.OR.ABS(SN(K,I+1)).LT.1.0E-05)
144	+GO TO 65
145	C
146	X1=X(1)+SN(K,I+1)*X(4)
147	X2=X(2)+SN(K,I+1)*X(5)
148	X3=X(3)+SN(K,I+1)*X(6)
149	CP=(X1-X0(I)-DXDZ(I)X3)(Y0(J)-Y0(I)+(DYDZ(J)-DYDZ(I))*X3)
150	+ - (X2-Y0(I)-DYDZ(I)X3)(X0(J)-X0(I)+(DXDZ(J)-DXDZ(I))*X3)
151	CP=CP/SQRT((X0(J)-X0(I)+(DXDZ(J)-DXDZ(I))X3)*2+
152	+ (Y0(J)-Y0(I)+(DYDZ(J)-DYDZ(I))X3)*2)
153	C
154	IF(ABS(CP).LT.0.0001) GO TO 65
155	IF(ABS(CP/SN(K,I+1)).LT.1.0E-06) GO TO 65
156	WRITE(CHMAIL,1020) I,K,SN(K,I+1)
157	CALL GMAIL(0,0)
158	WRITE(CHMAIL,1021) X
159	CALL GMAIL(0,0)
160	WRITE(CHMAIL,1022) X1,X2,X3,CP
161	CALL GMAIL(0,0)
162	WRITE(CHMAIL,1023) A,B,C,DISC
163	CALL GMAIL(0,0)
164	WRITE(CHMAIL,1024) INSIDE
165	CALL GMAIL(0,0)
166	WRITE(CHMAIL,1025) X0
167	CALL GMAIL(0,0)
168	WRITE(CHMAIL,1026) Y0
169	CALL GMAIL(0,0)
170	WRITE(CHMAIL,1027) DXDZ
171	CALL GMAIL(0,0)
172	WRITE(CHMAIL,1028) DYDZ
173	CALL GMAIL(0,0)
174	1020 FORMAT('0 GSNGTR ERROR - I,K =',2I2,' SN =',E13.5)
175	1021 FORMAT(' X =',6F15.6)
176	1022 FORMAT(' X1,X2,X3 =',3F15.6,' CP =',E15.6)
177	1023 FORMAT(' A =',E15.6,' B =',E15.6,' C =',E15.6,' DISC =',E15.6)
178	1024 FORMAT(' INSIDE =',I3)
179	1025 FORMAT(' X0 =',4E15.6)
180	1026 FORMAT(' Y0 =',4E15.6)
181	1027 FORMAT(' DXDZ =',4E15.6)
182	1028 FORMAT(' DYDZ =',4E15.6)
183	C
184	65 CONTINUE
185	C
186	C
187	C Have computed the two distances for the z planes and
188	C the four surfaces. Combine them accordingly as to
189	C whether the point is inside or outside the shape.
190	C
191	IF(INSIDE.EQ.0) GO TO 80
192	C
193	C Point is inside shape.
194	C
195	DO 70 I=1,5
196	DO 70 J=1,2
197	IF(SN(J,I).GT.0.0.AND.SN(J,I).LT.SNXT) SNXT=SN(J,I)
198	70 CONTINUE
199	GO TO 999
200	C
201	80 CONTINUE
202	C
203	C Point is outside shape.
204	C
205	IOUT(1)=0
206	IF(ABS(X(3)).GT.P(1)) IOUT(1)=1
207	C
208	C For each of five sets of SN and IOUT, IOUT(I) equal to 1
209	C indicates that the point is outside the shape by the Ith
210	C test, SN(1,I) is the distance to the first change in the
211	C test and SN(2,I) is the distance to the second change.
212	C The remaining logic just attempts to find a distance when
213	C the line is inside by all five tests, bearing in mind that
214	C for some tests the line can start inside, leave and return
215	C inside.
216	C
217	SL=-1.0
218	SM=BIG
219	SM1=BIG
220	DO 100 I=1,5
221	IF(IOUT(I).EQ.0) GO TO 90
222	IF(SN(2,I).LT.0.0) GO TO 999
223	IF(SN(1,I).LT.0.0.AND.SN(2,I).GT.SL) SL=SN(2,I)
224	IF(SN(1,I).GT.SL) SL=SN(1,I)
225	IF(SN(1,I).GE.0.0.AND.SN(2,I).LT.SM) SM=SN(2,I)
226	GO TO 100
227	90 CONTINUE
228	IF(SN(1,I).LT.0.0.AND.SN(2,I).GE.0.0.AND.SN(2,I).LT.SM) SM=SN(2,I)
229	IF(SN(1,I).LT.0.0.OR.SN(1,I).GT.SM1) GO TO 100
230	IF(SN(1,I).GE.SN(2,I)) GO TO 100
231	SM1=SN(1,I)
232	SL1=SN(2,I)
233	100 CONTINUE
234	C
235	C SL is the largest of the five distances to the first
236	C time the line is inside. SM is the smallest to the
237	C last time the point is inside. SM1 is the smallest
238	C distance to when the line is temporarily outside
239	C one of the tests.
240	C
241	IF(SM.LE.SL) GO TO 999
242	IF(SM1.GT.SL) GO TO 130
243	C
244	110 CONTINUE
245	C
246	C In this loop SL is updated by the return after SM1
247	C if SM1 is less than SL.
248	C
249	SL=SL1
250	IF(SM.LE.SL) GO TO 999
251	SM1=SM
252	C
253	DO 120 I=1,5
254	IF(IOUT(I).EQ.1) GO TO 120
255	IF(SN(2,I).LE.SL.OR.SN(1,I).GT.SM1) GO TO 120
256	IF(SN(1,I).GE.SN(2,I)) GO TO 120
257	SM1=SN(1,I)
258	SL1=SN(2,I)
259	120 CONTINUE
260	C
261	IF(SM1.GT.SL) GO TO 130
262	C
263	GO TO 110
264	130 CONTINUE
265	C
266	IF(SL.LT.SNXT) SNXT=SL
267	C
268	999 CONTINUE
269	END