[u/mrichter/AliRoot.git] / GEANT321 / gphys / glisur.F

*PATCH,GEXAM1.
*CMZ :  3.21/02 29/03/94  15.41.35  by  S.Giani
*DECK,GLISUR
*CMZ :          12/09/95  11.22.49  by  S.Ravndal
*CMZ :  3.21/02 03/07/94  19.08.44  by  S.Giani
*-- Author :
      SUBROUTINE GLISUR(X0,X1,MEDI0,MEDI1,U,PDOTU,IERR)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *  This routine simulates the surface profile of a boundary      *
C.    *  between two media as seen by an approaching particle with     *
C.    *  coordinates and direction given by X0. The surface is         *
C.    *  identified by the arguments MEDI0 and MEDI1 which are the     *
C.    *  media indices of the region in which the track is presently   *
C.    *  and the one which it approaches, respectively. The input      *
C.    *  vector X1 contains the coordinates of a point on the other    *
C.    *  side of the boundary from X0, and lying within medium MEDI1.  *
C.    *  The result is a unit vector U normal to the surface of        *
C.    *  reflection at X0 and pointing into the medium from which the  *
C.    *  track is approaching. The quality of the surface finish is    *
C.    *  described by the parameter POLISH, which varies from 0 to 1.  *
C.    *  POLISH=0 means maximum roughness, with effective plane of     *
C.    *  reflection distributed as cos(alpha), where alpha is the      *
C.    *  angle between the unit normal to the effective plane of       *
C.    *  reflection and the normal to the nominal medium boundary at   *
C.    *  X0. POLISH=1 means perfect smoothness. In between the surface *
C.    *  is modeled by a bell-shaped distribution in alpha, with       *
C.    *  limits defined by                                             *
C.    *              sin(alpha) = +/- (1-POLISH)                       *
C.    *  At the interface between two media, the surface POLISH        *
C.    *  parameter is taken from a user routine GUPLSH(MEDI0,MEDI1).   *
C.    *  The indices MEDI0 and MEDI1 refer to the media declared by    *
C.    *  the user. If GGPERP returns an error, there was a precision   *
C.    *  problem with the tracking, and point X0 is on the same side   *
C.    *  of the surface as X1. In this case (or any other error        *
C.    *  condition from GGPERP) a return is made with IERR ^=0. If     *
C.    *  IERR=0 on return then U contains a valid unit vector.         *
C.    *                                                                *
C.    *    ==>Called by : <USER>, UGINIT                               *
C.    *       Author    F.Carminati, R.Jones ***********               *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gcvolu.inc"
#include "geant321/gcvol1.inc"
C
      REAL X0(6),X1(3),U(3)
      REAL D(3),RNDM(3)

*
* *** Decide which volume defines the surface
      IF (NLEVEL.GE.NLEVL1) THEN
         LVOL = 1
      ELSE
         LVOL = -1
         DO 10 I=2,NLEVEL
            IF ((NAMES1(I).NE.NAMES(I)).OR. (NUMBR1(I).NE.NUMBER(I)))
     +      LVOL = 1
   10    CONTINUE
      END IF
*
* *** Find the nominal unit normal to the surface
      IF (LVOL.EQ.1) THEN
C
C Instead of X0, now X1 is used by GGPERP, in order to calculate
C the right perpendicular even near to a corner,
C see the modifactions in GTCKOV concerning VBOU
C
         CALL GGPERP(X1,U,IERR)
      ELSE
         LVOL = NLEVEL
         NLEVEL = NLEVL1
         CALL GGPERP(X1,U,IERR)
         U(1) = -U(1)
         U(2) = -U(2)
         U(3) = -U(3)
         NLEVEL = LVOL
      END IF
      IF (IERR.NE.0) GO TO 999
      PDOTU = X0(4)*U(1) + X0(5)*U(2) + X0(6)*U(3)
      POLISH=GUPLSH(MEDI0,MEDI1)
      IF (POLISH.LT.1.) THEN
*
* *** Generate distortion vector D within a uniform sphere
         DIAM = 2.*(1.-POLISH)
         DIA2 = DIAM**2
         IROUND = 0
   20    IF(IROUND.GT.5) GO TO 999
         IROUND = IROUND+1
         CALL GRNDM(RNDM,3)
         D(1) = RNDM(1)-0.5
         D(2) = RNDM(2)-0.5
         D(3) = RNDM(3)-0.5
         D2 = D(1)**2+D(2)**2+D(3)**2
         IF (D2.GT.0.25) GO TO 20
         D(1) = D(1)*DIAM
         D(2) = D(2)*DIAM
         D(3) = D(3)*DIAM
         D2 = D2*DIA2
*
* *** Insure that V=U+D will cause reflection away from surface
         PDOTD = X0(4)*D(1) + X0(5)*D(2) + X0(6)*D(3)
         PDOTV = PDOTU+PDOTD
         UDOTD = U(1)*D(1) + U(2)*D(2) + U(3)*D(3)
         V2 = 1+D2+2*UDOTD
         IF (PDOTD*V2+PDOTV*(1.-D2).GT.0.) GO TO 20
*
* *** Normalize V and call it U
         VABS = 1./SQRT(V2)
         U(1) = (U(1)+D(1))*VABS
         U(2) = (U(2)+D(2))*VABS
         U(3) = (U(3)+D(3))*VABS
         PDOTU = PDOTV*VABS
      END IF
  999 END
Commit	Line	Data
fe4da5cc	1	*PATCH,GEXAM1.
	2	*CMZ : 3.21/02 29/03/94 15.41.35 by S.Giani
	3	*DECK,GLISUR
	4	*CMZ : 12/09/95 11.22.49 by S.Ravndal
	5	*CMZ : 3.21/02 03/07/94 19.08.44 by S.Giani
	6	*-- Author :
	7	SUBROUTINE GLISUR(X0,X1,MEDI0,MEDI1,U,PDOTU,IERR)
	8	C.
	9	C. ******************************************************************
	10	C. * *
	11	C. * This routine simulates the surface profile of a boundary *
	12	C. * between two media as seen by an approaching particle with *
	13	C. * coordinates and direction given by X0. The surface is *
	14	C. * identified by the arguments MEDI0 and MEDI1 which are the *
	15	C. * media indices of the region in which the track is presently *
	16	C. * and the one which it approaches, respectively. The input *
	17	C. * vector X1 contains the coordinates of a point on the other *
	18	C. * side of the boundary from X0, and lying within medium MEDI1. *
	19	C. * The result is a unit vector U normal to the surface of *
	20	C. * reflection at X0 and pointing into the medium from which the *
	21	C. * track is approaching. The quality of the surface finish is *
	22	C. * described by the parameter POLISH, which varies from 0 to 1. *
	23	C. * POLISH=0 means maximum roughness, with effective plane of *
	24	C. * reflection distributed as cos(alpha), where alpha is the *
	25	C. * angle between the unit normal to the effective plane of *
	26	C. * reflection and the normal to the nominal medium boundary at *
	27	C. * X0. POLISH=1 means perfect smoothness. In between the surface *
	28	C. * is modeled by a bell-shaped distribution in alpha, with *
	29	C. * limits defined by *
	30	C. * sin(alpha) = +/- (1-POLISH) *
	31	C. * At the interface between two media, the surface POLISH *
	32	C. * parameter is taken from a user routine GUPLSH(MEDI0,MEDI1). *
	33	C. * The indices MEDI0 and MEDI1 refer to the media declared by *
	34	C. * the user. If GGPERP returns an error, there was a precision *
	35	C. * problem with the tracking, and point X0 is on the same side *
	36	C. * of the surface as X1. In this case (or any other error *
	37	C. * condition from GGPERP) a return is made with IERR ^=0. If *
	38	C. * IERR=0 on return then U contains a valid unit vector. *
	39	C. * *
	40	C. * ==>Called by : <USER>, UGINIT *
	41	C. * Author F.Carminati, R.Jones *********** *
	42	C. * *
	43	C. ******************************************************************
	44	C.
	45	#include "geant321/gcvolu.inc"
	46	#include "geant321/gcvol1.inc"
	47	C
	48	REAL X0(6),X1(3),U(3)
	49	REAL D(3),RNDM(3)
	50
	51	*
	52	* *** Decide which volume defines the surface
	53	IF (NLEVEL.GE.NLEVL1) THEN
	54	LVOL = 1
	55	ELSE
	56	LVOL = -1
	57	DO 10 I=2,NLEVEL
	58	IF ((NAMES1(I).NE.NAMES(I)).OR. (NUMBR1(I).NE.NUMBER(I)))
	59	+ LVOL = 1
	60	10 CONTINUE
	61	END IF
	62	*
	63	* *** Find the nominal unit normal to the surface
	64	IF (LVOL.EQ.1) THEN
65	C
66	C Instead of X0, now X1 is used by GGPERP, in order to calculate
67	C the right perpendicular even near to a corner,
68	C see the modifactions in GTCKOV concerning VBOU
69	C
70	CALL GGPERP(X1,U,IERR)
71	ELSE
72	LVOL = NLEVEL
73	NLEVEL = NLEVL1
74	CALL GGPERP(X1,U,IERR)
75	U(1) = -U(1)
76	U(2) = -U(2)
77	U(3) = -U(3)
78	NLEVEL = LVOL
79	END IF
80	IF (IERR.NE.0) GO TO 999
81	PDOTU = X0(4)U(1) + X0(5)U(2) + X0(6)*U(3)
82	POLISH=GUPLSH(MEDI0,MEDI1)
83	IF (POLISH.LT.1.) THEN
84	*
85	* *** Generate distortion vector D within a uniform sphere
86	DIAM = 2.*(1.-POLISH)
87	DIA2 = DIAM**2
88	IROUND = 0
89	20 IF(IROUND.GT.5) GO TO 999
90	IROUND = IROUND+1
91	CALL GRNDM(RNDM,3)
92	D(1) = RNDM(1)-0.5
93	D(2) = RNDM(2)-0.5
94	D(3) = RNDM(3)-0.5
95	D2 = D(1)2+D(2)2+D(3)**2
96	IF (D2.GT.0.25) GO TO 20
97	D(1) = D(1)*DIAM
98	D(2) = D(2)*DIAM
99	D(3) = D(3)*DIAM
100	D2 = D2*DIA2
101	*
102	* *** Insure that V=U+D will cause reflection away from surface
103	PDOTD = X0(4)D(1) + X0(5)D(2) + X0(6)*D(3)
104	PDOTV = PDOTU+PDOTD
105	UDOTD = U(1)D(1) + U(2)D(2) + U(3)*D(3)
106	V2 = 1+D2+2*UDOTD
107	IF (PDOTDV2+PDOTV(1.-D2).GT.0.) GO TO 20
108	*
109	* *** Normalize V and call it U
110	VABS = 1./SQRT(V2)
111	U(1) = (U(1)+D(1))*VABS
112	U(2) = (U(2)+D(2))*VABS
113	U(3) = (U(3)+D(3))*VABS
114	PDOTU = PDOTV*VABS
115	END IF
116	999 END