+++ /dev/null
-*
-* $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:50 cernlib
-* Geant
-*
-*
-#include "geant321/pilot.h"
-*FCA : 05/01/99 09:58:02 by Federico Carminati
-* Effectively print the message when a shape is
-* not implemented
-*CMZ : 3.21/02 29/03/94 15.41.28 by S.Giani
-*-- Author :
- SUBROUTINE GGPERP (X,U,IERR)
-C.
-C. ****************************************************************
-C. * *
-C. * This routine solves the general problem of calculating the *
-C. * unit vector normal to the surface of the current volume at *
-C. * the point X. The result is returned in the array U. X is *
-C. * assumed to be on or near a boundary of the current volume. *
-C. * The current volume is indicated by the common /GCVOLU/. *
-C. * U points from inside to outside in that neighbourhood. *
-C. * If X is equidistant to more than one boundary (in a corner) *
-C. * an arbitrary choice is made based upon the order of *
-C. * precedence implied by the IF statements below. If the *
-C. * routine fails to find the unit normal, it returns with *
-C. * IERR=1, otherwise IERR=0. *
-C. * *
-C. * Called by : GSURFP, GDSTEP *
-C. * Authors : F.Carminati, R.Jones, F.Ohlsson-Malek *
-C. * *
-C. ****************************************************************
-#include "geant321/gcvolu.inc"
-#include "geant321/gconsp.inc"
-#include "geant321/gcbank.inc"
-#include "geant321/gcshno.inc"
-#include "geant321/gctmed.inc"
-#include "geant321/gcunit.inc"
- DIMENSION X(3),U(3),XL(3),UL(3),DXL(3),PAR(100),SPAR(100),ATT(20)
- DIMENSION PERP(10)
-#if !defined(CERNLIB_SINGLE)
- DOUBLE PRECISION PERP,PMIN0
- DOUBLE PRECISION PAR,DXL,RHO,R,RINV,PHI,THE
- DOUBLE PRECISION PHI1,PHI2,THE1,THE2,XWID
- DOUBLE PRECISION GUARD,DPHI,PHI0,SPHI0,CPHI0
- DOUBLE PRECISION FACT,CALPH,SALPH,TALPH
- DOUBLE PRECISION RAT,RATL,RATH,H,BL,TL,DX,DY,DZ,DU
- DOUBLE PRECISION UU0,VV0,UU,W1,W2,W3,W4
- DOUBLE PRECISION SEW1,SEW2,SEW3,SEW4
- DOUBLE PRECISION TAN1,TAN2,TAN3,TAN4
- DOUBLE PRECISION SEC1,SEC2,SEC3,SEC4
- DOUBLE PRECISION U0,V0,U1,U1L,U2,U2L
- DOUBLE PRECISION ONE,TWO
- DOUBLE PRECISION DSECT,ZERO,FULL,FULL10,DBY2
-#endif
- LOGICAL LNTFOU
- PARAMETER (ONE=1,TWO=2)
- PARAMETER (ZERO=0.,DBY2=0.5,FULL=360.,FULL10=3600.)
-C.
-C. ------------------------------------------------------------------
-C.
- LNTFOU = .FALSE.
-*
-* *** Transform current point into local reference system
- CALL GMTOD(X,XL,1)
- DXL(1) = XL(1)
- DXL(2) = XL(2)
- DXL(3) = XL(3)
-*
-* *** Fetch the parameters of the current volume
- JVO = LQ(JVOLUM-LVOLUM(NLEVEL))
- IN = LINDEX(NLEVEL)
- IF (NLEVEL.GT.1) THEN
- JVOM = LQ(JVOLUM-LVOLUM(NLEVEL-1))
- JIN = LQ(JVOM-IN)
- ENDIF
- ISH = Q(JVO+2)
- NIN = Q(JVO+3)
- IF (NLEVEL.LT.NLDEV(NLEVEL)) THEN
- JPAR = 0
- ELSE
-* (case with structure JVOLUM locally developed)
- JPAR = LQ(LQ(JVOLUM-LVOLUM(NLDEV(NLEVEL))))
- IF (NLEVEL.EQ.NLDEV(NLEVEL)) GO TO 20
- DO 10 ILEV = NLDEV(NLEVEL), NLEVEL-1
- IF (IQ(JPAR+1).EQ.0) THEN
- JPAR = LQ(JPAR-LINDEX(ILEV+1))
- IF (JPAR.EQ.0) GO TO 20
- ELSE IF (IQ(JPAR-3).GT.1) THEN
- JPAR = LQ(JPAR-LINDEX(ILEV+1))
- ELSE
- JPAR = LQ(JPAR-1)
- ENDIF
- IF (ILEV.EQ.NLEVEL-1) THEN
- JPAR = JPAR + 5
- NPAR = IQ(JPAR)
- CALL UCOPY (Q(JPAR+1), SPAR, NPAR)
- DO 100 I=1,NPAR
- PAR(I)=SPAR(I)
- 100 CONTINUE
- ENDIF
- 10 CONTINUE
- GO TO 30
- ENDIF
-* (normal case)
- 20 CONTINUE
- CALL GFIPAR(JVO,JIN,IN,NPAR,NATT,SPAR,ATT)
- DO 101 I=1,NPAR
- PAR(I)=SPAR(I)
- 101 CONTINUE
- 30 CONTINUE
-*
-* *** Case of the BOX:
- IF (ISH.EQ.NSBOX) THEN
- PERP(1) = ABS(ABS(DXL(1))-PAR(1))
- PERP(2) = ABS(ABS(DXL(2))-PAR(2))
- PERP(3) = ABS(ABS(DXL(3))-PAR(3))
- PMIN0 = MIN(PERP(1),PERP(2),PERP(3))
- IF (PERP(1).EQ.PMIN0) THEN
- UL(1) = SIGN(ONE,DXL(1))
- UL(2) = 0.
- UL(3) = 0.
- ELSE IF (PERP(2).EQ.PMIN0) THEN
- UL(1) = 0.
- UL(2) = SIGN(ONE,DXL(2))
- UL(3) = 0.
- ELSE IF (PERP(3).EQ.PMIN0) THEN
- UL(1) = 0.
- UL(2) = 0.
- UL(3) = SIGN(ONE,DXL(3))
- ELSE
- LNTFOU=.TRUE.
- ENDIF
-*
-* *** Case of the TUBE, TUBeSection:
- ELSE IF (ISH.EQ.NSTUBE.OR.ISH.EQ.NSTUBS) THEN
- RHO = SQRT(DXL(1)**2 + DXL(2)**2)
- PERP(1) = ABS(RHO-PAR(1))
- PERP(2) = ABS(RHO-PAR(2))
- PERP(3) = ABS(ABS(DXL(3))-PAR(3))
- IF (ISH.EQ.NSTUBE) THEN
- PMIN0 = MIN(PERP(1),PERP(2),PERP(3))
- ELSE
- PHI = ATAN2(DXL(2),DXL(1))
- IF (PHI.LT.0.) PHI = PHI+TWOPI
- PHI1 = MOD(PAR(4)+FULL10,FULL)*DEGRAD
- PERP(4) = ABS(PHI-PHI1)
- IF (PERP(4).GT.PI) PERP(4) = TWOPI-PERP(4)
- PHI2 = MOD(PAR(5)+FULL10,FULL)*DEGRAD
- PERP(5) = ABS(PHI-PHI2)
- IF (PERP(5).GT.PI) PERP(5) = TWOPI-PERP(5)
- PERP(4) = PERP(4)*RHO
- PERP(5) = PERP(5)*RHO
- PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),PERP(5))
- ENDIF
- IF (PERP(1).EQ.PMIN0) THEN
- UL(1) = -DXL(1)/RHO
- UL(2) = -DXL(2)/RHO
- UL(3) = 0.
- ELSE IF (PERP(2).EQ.PMIN0) THEN
- UL(1) = DXL(1)/RHO
- UL(2) = DXL(2)/RHO
- UL(3) = 0.
- ELSE IF (PERP(3).EQ.PMIN0) THEN
- UL(1) = 0.
- UL(2) = 0.
- UL(3) = SIGN(ONE,DXL(3))
- ELSE IF (PERP(4).EQ.PMIN0) THEN
- UL(1) = SIN(PHI1)
- UL(2) = -COS(PHI1)
- UL(3) = 0.
- ELSE IF (PERP(5).EQ.PMIN0) THEN
- UL(1) = -SIN(PHI2)
- UL(2) = COS(PHI2)
- UL(3) = 0.
- ELSE
- LNTFOU=.TRUE.
- ENDIF
-*
-* *** Case of the CONE, CONeSection:
- ELSE IF (ISH.EQ.NSCONE.OR.ISH.EQ.NSCONS) THEN
- RHO = SQRT(DXL(1)**2 + DXL(2)**2)
- TAN1 = (PAR(4)-PAR(2))/(TWO*PAR(1))
- SEC1 = SQRT(ONE+TAN1**2)
- U1 = RHO-DXL(3)*TAN1
- U1L = PAR(4)-PAR(1)*TAN1
- TAN2 = (PAR(5)-PAR(3))/(TWO*PAR(1))
- SEC2 = SQRT(ONE+TAN2**2)
- U2 = RHO-DXL(3)*TAN2
- U2L = PAR(5)-PAR(1)*TAN2
- PERP(1) = ABS(ABS(DXL(3))-PAR(1))
- PERP(2) = ABS(U1-U1L)/SEC1
- PERP(3) = ABS(U2-U2L)/SEC2
- IF (ISH.EQ.NSCONE) THEN
- PMIN0 = MIN(PERP(1),PERP(2),PERP(3))
- ELSE
- PHI = ATAN2(DXL(2),DXL(1))
- IF (PHI.LT.0.) PHI = PHI+TWOPI
- PHI1 = MOD(PAR(6)+FULL10,FULL)*DEGRAD
- PERP(4) = ABS(PHI-PHI1)
- IF (PERP(4).GT.PI) PERP(4) = TWOPI-PERP(4)
- PHI2 = MOD(PAR(7)+FULL10,FULL)*DEGRAD
- PERP(5) = ABS(PHI-PHI2)
- IF (PERP(5).GT.PI) PERP(5) = TWOPI-PERP(5)
- PERP(4) = PERP(4)*RHO
- PERP(5) = PERP(5)*RHO
- PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),PERP(5))
- ENDIF
- IF (PERP(1).EQ.PMIN0) THEN
- UL(1) = 0.
- UL(2) = 0.
- UL(3) = SIGN(ONE,DXL(3))
- ELSE IF (PERP(2).EQ.PMIN0) THEN
- RHO = RHO*SEC1
- UL(1) = -DXL(1)/RHO
- UL(2) = -DXL(2)/RHO
- UL(3) = TAN1/SEC1
- ELSE IF (PERP(3).EQ.PMIN0) THEN
- RHO = RHO*SEC2
- UL(1) = DXL(1)/RHO
- UL(2) = DXL(2)/RHO
- UL(3) = -TAN2/SEC2
- ELSE IF (PERP(4).EQ.PMIN0) THEN
- UL(1) = SIN(PHI1)
- UL(2) = -COS(PHI1)
- UL(3) = 0.
- ELSE IF (PERP(5).EQ.PMIN0) THEN
- UL(1) = -SIN(PHI2)
- UL(2) = COS(PHI2)
- UL(3) = 0.
- ELSE
- LNTFOU=.TRUE.
- ENDIF
-*
-* *** Case of the PolyCONe:
- ELSE IF (ISH.EQ.NSPCON) THEN
- PERP(1) = ABS(DXL(3)-PAR(4))
- DO 400 I=7,NPAR,3
- PERP(2) = ABS(DXL(3)-PAR(I))
- IF (PERP(2).GT.PERP(1)) GOTO 401
- PERP(1) = PERP(2)
- 400 CONTINUE
- 401 I = I-3
- IF (I.GT.4) THEN
- PERP(1) = 100.
- RHO = SQRT(DXL(1)**2 + DXL(2)**2)
- DZ = PAR(I)-PAR(I-3)+1.e-10
- TAN1 = (PAR(I+1)-PAR(I-2))/DZ
- SEC1 = SQRT(ONE+TAN1**2)
- U1 = RHO-DXL(3)*TAN1
- U1L = PAR(I+1)-PAR(I)*TAN1
- TAN2 = (PAR(I+2)-PAR(I-1))/DZ
- SEC2 = SQRT(ONE+TAN2**2)
- U2 = RHO-DXL(3)*TAN2
- U2L = PAR(I+2)-PAR(I)*TAN2
- GUARD = MAX(DXL(3)-PAR(I),ZERO)
- PERP(3) = ABS(U1-U1L)/SEC1 + GUARD*SEC1
- PERP(4) = ABS(U2-U2L)/SEC2 + GUARD*SEC2
- ELSE
- PERP(3) = 100.
- PERP(4) = 100.
- ENDIF
- IF (I.LT.NPAR-2) THEN
- PERP(2) = 100.
- RHO = SQRT(DXL(1)**2 + DXL(2)**2)
- DZ = PAR(I+3)-PAR(I)+1.e-10
- TAN3 = (PAR(I+4)-PAR(I+1))/DZ
- SEC3 = SQRT(ONE+TAN3**2)
- U1 = RHO-DXL(3)*TAN3
- U1L = PAR(I+1)-PAR(I)*TAN3
- TAN4 = (PAR(I+5)-PAR(I+2))/DZ
- SEC4 = SQRT(ONE+TAN4**2)
- U2 = RHO-DXL(3)*TAN4
- U2L = PAR(I+2)-PAR(I)*TAN4
- GUARD = MAX(PAR(I)-DXL(3),ZERO)
- PERP(5) = ABS(U1-U1L)/SEC3 + GUARD*SEC3
- PERP(6) = ABS(U2-U2L)/SEC4 + GUARD*SEC4
- ELSE
- PERP(5) = 100.
- PERP(6) = 100.
- ENDIF
- PHI = ATAN2(DXL(2),DXL(1))
- IF (PHI.LT.0.) PHI = PHI+TWOPI
- PHI1 = MOD(PAR(1)+FULL10,FULL)*DEGRAD
- PERP(7) = ABS(PHI-PHI1)
- IF (PERP(7).GT.PI) PERP(7) = TWOPI-PERP(7)
- PHI2 = MOD(PAR(1)+PAR(2)+FULL10,FULL)*DEGRAD
- PERP(8) = ABS(PHI-PHI2)
- IF (PERP(8).GT.PI) PERP(8) = TWOPI-PERP(8)
- PERP(7) = PERP(7)*RHO
- PERP(8) = PERP(8)*RHO
- PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),
- + PERP(5),PERP(6),PERP(7),PERP(8))
- IF (PERP(1).EQ.PMIN0) THEN
- UL(1) = 0.
- UL(2) = 0.
- UL(3) = -1.
- ELSE IF (PERP(2).EQ.PMIN0) THEN
- UL(1) = 0.
- UL(2) = 0.
- UL(3) = 1.
- ELSE IF (PERP(3).EQ.PMIN0) THEN
- RHO = RHO*SEC1
- UL(1) = -DXL(1)/RHO
- UL(2) = -DXL(2)/RHO
- UL(3) = TAN1/SEC1
- ELSE IF (PERP(4).EQ.PMIN0) THEN
- RHO = RHO*SEC2
- UL(1) = DXL(1)/RHO
- UL(2) = DXL(2)/RHO
- UL(3) = -TAN2/SEC2
- ELSE IF (PERP(5).EQ.PMIN0) THEN
- RHO = RHO*SEC3
- UL(1) = -DXL(1)/RHO
- UL(2) = -DXL(2)/RHO
- UL(3) = TAN3/SEC3
- ELSE IF (PERP(6).EQ.PMIN0) THEN
- RHO = RHO*SEC4
- UL(1) = DXL(1)/RHO
- UL(2) = DXL(2)/RHO
- UL(3) = -TAN4/SEC4
- ELSE IF (PERP(7).EQ.PMIN0) THEN
- UL(1) = SIN(PHI1)
- UL(2) = -COS(PHI1)
- UL(3) = 0.
- ELSE IF (PERP(8).EQ.PMIN0) THEN
- UL(1) = -SIN(PHI2)
- UL(2) = COS(PHI2)
- UL(3) = 0.
- ELSE
- LNTFOU=.TRUE.
- ENDIF
-*
-* *** Case of the PolyGON:
- ELSE IF (ISH.EQ.NSPGON) THEN
- RHO = SQRT(DXL(1)**2+DXL(2)**2)
- PHI = ATAN2(DXL(2),DXL(1))
- IF (PHI.LT.0.) PHI = PHI+TWOPI
- DPHI = MOD(PHI*RADDEG-PAR(1)+FULL10,FULL)
- PDIV = PAR(2)/PAR(3)
- DSECT = INT(DPHI/PDIV + ONE)
- IF (DSECT.GT.PAR(3)) THEN
- IF (DPHI.GT.(180.+PAR(2)*DBY2)) THEN
- DSECT = ONE
- ELSE
- DSECT = PAR(3)
- ENDIF
- ENDIF
- PHI0 = MOD(PAR(1)+(DSECT-DBY2)*PDIV+FULL10,FULL)*DEGRAD
- SPHI0 = SIN(PHI0)
- CPHI0 = COS(PHI0)
- U0 = DXL(1)*CPHI0 + DXL(2)*SPHI0
- V0 = DXL(2)*CPHI0 - DXL(1)*SPHI0
- PERP(1) = ABS(DXL(3)-PAR(5))
- DO 500 I=8,NPAR,3
- PERP(2) = ABS(DXL(3)-PAR(I))
- IF (PERP(2).GT.PERP(1)) GOTO 501
- PERP(1) = PERP(2)
- 500 CONTINUE
- 501 I = I-3
- IF (I.GT.5) THEN
- PERP(1) = 100.
- DZ = PAR(I)-PAR(I-3)+1.e-10
- TAN1 = (PAR(I+1)-PAR(I-2))/DZ
- SEC1 = SQRT(ONE+TAN1**2)
- U1 = U0-DXL(3)*TAN1
- U1L = PAR(I+1)-PAR(I)*TAN1
- TAN2 = (PAR(I+2)-PAR(I-1))/DZ
- SEC2 = SQRT(ONE+TAN2**2)
- U2 = U0-DXL(3)*TAN2
- U2L = PAR(I+2)-PAR(I)*TAN2
- GUARD = MAX(DXL(3)-PAR(I),ZERO)
- PERP(3) = ABS(U1-U1L)/SEC1 + GUARD*SEC1
- PERP(4) = ABS(U2-U2L)/SEC2 + GUARD*SEC2
- ELSE
- PERP(3) = 100.
- PERP(4) = 100.
- ENDIF
- IF (I.LT.NPAR-2) THEN
- PERP(2) = 100.
- DZ = PAR(I+3)-PAR(I)+1.e-10
- TAN3 = (PAR(I+4)-PAR(I+1))/DZ
- SEC3 = SQRT(ONE+TAN3**2)
- U1 = U0-DXL(3)*TAN3
- U1L = PAR(I+1)-PAR(I)*TAN3
- TAN4 = (PAR(I+5)-PAR(I+2))/DZ
- SEC4 = SQRT(ONE+TAN4**2)
- U2 = U0-DXL(3)*TAN4
- U2L = PAR(I+2)-PAR(I)*TAN4
- GUARD = MAX(PAR(I)-DXL(3),ZERO)
- PERP(5) = ABS(U1-U1L)/SEC3 + GUARD*SEC3
- PERP(6) = ABS(U2-U2L)/SEC4 + GUARD*SEC4
- ELSE
- PERP(5) = 100.
- PERP(6) = 100.
- ENDIF
- PHI1 = MOD(PAR(1)+FULL10,FULL)*DEGRAD
- PERP(7) = ABS(PHI-PHI1)
- IF (PERP(7).GT.PI) PERP(7) = TWOPI-PERP(7)
- PHI2 = MOD(PAR(1)+PAR(2)+FULL10,FULL)*DEGRAD
- PERP(8) = ABS(PHI-PHI2)
- IF (PERP(8).GT.PI) PERP(8) = TWOPI-PERP(8)
- PERP(7) = PERP(7)*RHO
- PERP(8) = PERP(8)*RHO
- PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),
- + PERP(5),PERP(6),PERP(7),PERP(8))
- IF (PERP(1).EQ.PMIN0) THEN
- UL(1) = 0.
- UL(2) = 0.
- UL(3) = -1.
- ELSE IF (PERP(2).EQ.PMIN0) THEN
- UL(1) = 0.
- UL(2) = 0.
- UL(3) = 1.
- ELSE IF (PERP(3).EQ.PMIN0) THEN
- FACT = ONE/SEC1
- UL(1) = -CPHI0*FACT
- UL(2) = -SPHI0*FACT
- UL(3) = TAN1*FACT
- ELSE IF (PERP(4).EQ.PMIN0) THEN
- FACT = ONE/SEC2
- UL(1) = CPHI0*FACT
- UL(2) = SPHI0*FACT
- UL(3) = -TAN2*FACT
- ELSE IF (PERP(5).EQ.PMIN0) THEN
- FACT = ONE/SEC3
- UL(1) = -CPHI0*FACT
- UL(2) = -SPHI0*FACT
- UL(3) = TAN3*FACT
- ELSE IF (PERP(6).EQ.PMIN0) THEN
- FACT = ONE/SEC4
- UL(1) = CPHI0*FACT
- UL(2) = SPHI0*FACT
- UL(3) = -TAN4*FACT
- ELSE IF (PERP(7).EQ.PMIN0) THEN
- UL(1) = SIN(PHI1)
- UL(2) = -COS(PHI1)
- UL(3) = 0.
- ELSE IF (PERP(8).EQ.PMIN0) THEN
- UL(1) = -SIN(PHI2)
- UL(2) = COS(PHI2)
- UL(3) = 0.
- ELSE
- LNTFOU=.TRUE.
- ENDIF
-*
-* *** Case of the SPHEre:
- ELSE IF (ISH.EQ.NSSPHE) THEN
- R = SQRT(DXL(1)**2+DXL(2)**2+DXL(3)**2)
- RHO = SQRT(DXL(1)**2+DXL(2)**2)
- THE = ATAN2(RHO,DXL(3))
- PHI = ATAN2(DXL(2),DXL(1))
- IF (PHI.LT.0.) PHI = PHI+TWOPI
- THE1 = MOD(PAR(3)+FULL10,FULL)*DEGRAD
- THE2 = MOD(PAR(4)+FULL10,FULL)*DEGRAD
- PHI1 = MOD(PAR(5)+FULL10,FULL)*DEGRAD
- PHI2 = MOD(PAR(6)+FULL10,FULL)*DEGRAD
- PERP(1) = ABS(R-PAR(1))
- PERP(2) = ABS(R-PAR(2))
- PERP(3) = ABS(THE-THE1)*R
- PERP(4) = ABS(THE-THE2)*R
- PERP(5) = ABS(PHI-PHI1)
- IF (PERP(5).GT.PI) PERP(5) = TWOPI-PERP(5)
- PERP(5) = PERP(5)*RHO
- PERP(6) = ABS(PHI-PHI2)
- IF (PERP(6).GT.PI) PERP(6) = TWOPI-PERP(6)
- PERP(6) = PERP(6)*RHO
- PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),PERP(5),PERP(6))
- IF (PERP(1).EQ.PMIN0) THEN
- RINV = ONE/R
- UL(1) = -DXL(1)*RINV
- UL(2) = -DXL(2)*RINV
- UL(3) = -DXL(3)*RINV
- ELSE IF (PERP(2).EQ.PMIN0) THEN
- RINV = ONE/R
- UL(1) = DXL(1)*RINV
- UL(2) = DXL(2)*RINV
- UL(3) = DXL(3)*RINV
- ELSE IF (PERP(3).EQ.PMIN0) THEN
- UL(1) = -COS(THE1)*COS(PHI)
- UL(2) = -COS(THE1)*SIN(PHI)
- UL(3) = +SIN(THE1)
- ELSE IF (PERP(4).EQ.PMIN0) THEN
- UL(1) = +COS(THE2)*COS(PHI)
- UL(2) = +COS(THE2)*SIN(PHI)
- UL(3) = -SIN(THE2)
- ELSE IF (PERP(5).EQ.PMIN0) THEN
- UL(1) = +SIN(PHI1)
- UL(2) = -COS(PHI1)
- UL(3) = 0
- ELSE IF (PERP(6).EQ.PMIN0) THEN
- UL(1) = -SIN(PHI2)
- UL(2) = +COS(PHI2)
- UL(3) = 0
- ELSE
- LNTFOU=.TRUE.
- ENDIF
-*
-* *** Case of the PARAllelpiped:
-***************************************************************
-* Warning: the parameters for this shape are NOT stored in *
-* the data structure as the user supplies them. Rather, the *
-* user supplies PAR(4)=alph, PAR(5)=the, PAR(6)=phi, and the *
-* data structure contains PAR(4)=Tan(alph), PAR(5)=Tan(the)* *
-* Cos(phi), PAR(6)=Tan(the)*Sin(phi). *
-***************************************************************
- ELSE IF (ISH.EQ.NSPARA) THEN
- DX = PAR(5)
- DY = PAR(6)
- U0 = DXL(1)-DX*DXL(3)
- V0 = DXL(2)-DY*DXL(3)
- CALPH = ONE/SQRT(ONE+PAR(4)**2)
- SALPH = -CALPH*PAR(4)
- U1 = U0*CALPH+V0*SALPH
- U1L = PAR(1)*CALPH
- PERP(1) = ABS(ABS(U1)-U1L)
- PERP(2) = ABS(ABS(V0)-PAR(2))
- PERP(3) = ABS(ABS(DXL(3))-PAR(3))
- PMIN0 = MIN(PERP(1),PERP(2),PERP(3))
- IF (PERP(1).EQ.PMIN0) THEN
- DU = DX*CALPH+DY*SALPH
- FACT = SIGN(ONE/SQRT(ONE+DU**2),U1)
- UL(1) = CALPH*FACT
- UL(2) = SALPH*FACT
- UL(3) = -DU*FACT
- ELSE IF (PERP(2).EQ.PMIN0) THEN
- FACT = SIGN(ONE/SQRT(ONE+DY**2),V0)
- UL(1) = 0.
- UL(2) = FACT
- UL(3) = -DY*FACT
- ELSE IF (PERP(3).EQ.PMIN0) THEN
- UL(1) = 0.
- UL(2) = 0.
- UL(3) = SIGN(ONE,DXL(3))
- ELSE
- LNTFOU=.TRUE.
- ENDIF
-*
-* *** Case of the trapezoid TRD1
- ELSE IF (ISH.EQ.NSTRD1) THEN
- DZ = TWO*PAR(4)+1.e-10
- TAN1 = (PAR(2)-PAR(1))/DZ
- SEC1 = SQRT(ONE+TAN1**2)
- U1 = ABS(DXL(1))-DXL(3)*TAN1
- U1L = PAR(2)-PAR(4)*TAN1
- PERP(1) = ABS(U1-U1L)/SEC1
- PERP(2) = ABS(ABS(DXL(2))-PAR(3))
- PERP(3) = ABS(ABS(DXL(3))-PAR(4))
- PMIN0 = MIN(PERP(1),PERP(2),PERP(3))
- IF (PERP(1).EQ.PMIN0) THEN
- FACT = ONE/SEC1
- UL(1) = SIGN(FACT,DXL(1))
- UL(2) = 0.
- UL(3) = -TAN1*FACT
- ELSE IF (PERP(2).EQ.PMIN0) THEN
- UL(1) = 0.
- UL(2) = SIGN(ONE,DXL(2))
- UL(3) = 0.
- ELSE IF (PERP(3).EQ.PMIN0) THEN
- UL(1) = 0.
- UL(2) = 0.
- UL(3) = SIGN(ONE,DXL(3))
- ELSE
- LNTFOU=.TRUE.
- ENDIF
-*
-* *** Case of the trapezoid TRD2
- ELSE IF (ISH.EQ.NSTRD2) THEN
- DZ = TWO*PAR(5)+1.e-10
- TAN1 = (PAR(2)-PAR(1))/DZ
- SEC1 = SQRT(ONE+TAN1**2)
- U1 = ABS(DXL(1))-DXL(3)*TAN1
- U1L = PAR(2)-PAR(5)*TAN1
- TAN2 = (PAR(4)-PAR(3))/DZ
- SEC2 = SQRT(ONE+TAN2**2)
- U2 = ABS(DXL(2))-DXL(3)*TAN2
- U2L = PAR(4)-PAR(5)*TAN2
- PERP(1) = ABS(U1-U1L)/SEC1
- PERP(2) = ABS(U2-U2L)/SEC2
- PERP(3) = ABS(ABS(DXL(3))-PAR(5))
- PMIN0 = MIN(PERP(1),PERP(2),PERP(3))
- IF (PERP(1).EQ.PMIN0) THEN
- FACT = ONE/SEC1
- UL(1) = SIGN(FACT,DXL(1))
- UL(2) = 0.
- UL(3) = -TAN1*FACT
- ELSE IF (PERP(2).EQ.PMIN0) THEN
- FACT = ONE/SEC2
- UL(1) = 0.
- UL(2) = SIGN(FACT,DXL(2))
- UL(3) = -TAN2*FACT
- ELSE IF (PERP(3).EQ.PMIN0) THEN
- UL(1) = 0.
- UL(2) = 0.
- UL(3) = SIGN(ONE,DXL(3))
- ELSE
- LNTFOU=.TRUE.
- ENDIF
-*
-* *** Case of the TRAPezoid
-***************************************************************
-* Warning: the parameters for this shape are NOT stored in *
-* the data structure as the user supplies them. Rather, the *
-* user supplies PAR(2)=thet, PAR(3)=phi, PAR(7)=alp1, and *
-* PAR(11)=alp2, while the data structure contains PAR(2)= *
-* Tan(thet)*Cos(phi), PAR(3)=Tan(thet)*Sin(phi), PAR(7)= *
-* Tan(alp1), and PAR(11)=Tan(alp2). *
-***************************************************************
- ELSE IF (ISH.EQ.NSTRAP) THEN
- PERP(1) = ABS(ABS(DXL(3))-PAR(1))
- DX = PAR(2)
- DY = PAR(3)
- U0 = DX*DXL(3)
- V0 = DY*DXL(3)
- UU0 = DX*PAR(1)
- VV0 = DY*PAR(1)
- RAT = DXL(3)/PAR(1)
- RATL = (ONE-RAT)/TWO
- RATH = (ONE+RAT)/TWO
- H = PAR(4)*RATL+PAR(8)*RATH
- BL = PAR(5)*RATL+PAR(9)*RATH
- TL = PAR(6)*RATL+PAR(10)*RATH
- TALPH = PAR(7)*RATL+PAR(11)*RATH
- XWID = (TL+BL)/TWO
- TAN1 = TALPH+(TL-BL)/(TWO*H)
- SEC1 = SQRT(ONE+TAN1**2)
- U1 = DXL(1)-DXL(2)*TAN1
- U1L = U0+XWID-V0*TAN1
- TAN2 = TAN1-TWO*TALPH
- SEC2 = SQRT(ONE+TAN2**2)
- U2 = DXL(1)+DXL(2)*TAN2
- U2L = U0-XWID+V0*TAN2
- IF (DXL(3).LT.0) THEN
- DZ = PAR(1)-DXL(3)+1.e-10
- UU = UU0+(PAR(9)+PAR(10))/TWO
- W1 = (UU-VV0*TAN1-U1L)/DZ
- UU = TWO*UU0-UU
- W2 = (UU+VV0*TAN2-U2L)/DZ
- ELSE
- DZ = -PAR(1)-DXL(3)+1.e-10
- UU = -UU0+(PAR(5)+PAR(6))/TWO
- W1 = (UU+VV0*TAN1-U1L)/DZ
- UU = -TWO*UU0-UU
- W2 = (UU-VV0*TAN2-U2L)/DZ
- ENDIF
- W3 = DY+(PAR(8)-PAR(4))/(TWO*PAR(1))
- W4 = TWO*DY-W3
- SEW1 = SQRT(ONE+W1**2)
- SEW2 = SQRT(ONE+W2**2)
- SEW3 = SQRT(ONE+W3**2)
- SEW4 = SQRT(ONE+W4**2)
- PERP(2) = ABS(U1-U1L)/(SEC1*SEW1)
- PERP(3) = ABS(U2-U2L)/(SEC2*SEW2)
- PERP(4) = ABS(DXL(2)-V0-H)/SEW3
- PERP(5) = ABS(DXL(2)-V0+H)/SEW4
- PMIN0 = MIN(PERP(1),PERP(2),PERP(3),PERP(4),PERP(5))
- IF (PERP(1).EQ.PMIN0) THEN
- UL(1) = 0.
- UL(2) = 0.
- UL(3) = SIGN(ONE,DXL(3))
- ELSE IF (PERP(2).EQ.PMIN0) THEN
- FACT = ONE/(SEC1*SEW1)
- UL(1) = FACT
- UL(2) = -TAN1*FACT
- UL(3) = -W1/SEW1
- ELSE IF (PERP(3).EQ.PMIN0) THEN
- FACT = ONE/(SEC2*SEW2)
- UL(1) = -FACT
- UL(2) = -TAN2*FACT
- UL(3) = W2/SEW2
- ELSE IF (PERP(4).EQ.PMIN0) THEN
- FACT = ONE/SEW3
- UL(1) = 0.
- UL(2) = FACT
- UL(3) = -W3*FACT
- ELSE IF (PERP(5).EQ.PMIN0) THEN
- FACT = ONE/SEW4
- UL(1) = 0.
- UL(2) = -FACT
- UL(3) = W4*FACT
- ELSE
- LNTFOU=.TRUE.
- ENDIF
-*
-* *** everything else (currently NOT IMPLEMENTED)
- ELSE
- WRITE(CHMAIL,10100) ISH
- CALL GMAIL(0,0)
- IERR = 1
- GOTO 999
- ENDIF
-
- IF(LNTFOU) THEN
- WRITE(CHMAIL,10000) ISH
- CALL GMAIL(0,0)
- IERR = 1
- ELSE
-*
-* *** Transform back into the MCS
- CALL GDTOM(UL,U,2)
- IERR = 0
- ENDIF
-
-10000 FORMAT(' GGPERP - geometry check error for shape #',I2,'!')
-10100 FORMAT(' GGPERP - non implemented for shape #',I2)
- 999 END
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff