*
* $Id$
*
* $Log$
* Revision 1.1  1995/10/24 10:20:51  cernlib
* Initial revision
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.29  by  S.Giani
*-- Author :
      SUBROUTINE GNOPCO (X, PAR, IACT, SNEXT, SNXT, SAFE)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *       COMPUTE DISTANCE UP TO INTERSECTION WITH 'PCON' VOLUME,  *
C.    *       FROM OUTSIDE POINT X(1-3) ALONG DIRECTION X(4-6)         *
C.    *                                                                *
C.    *       PAR   (input)  : volume parameters                       *
C.    *       IACT  (input)  : action flag                             *
C.    *         = 0  Compute SAFE only                                 *
C.    *         = 1  Compute SAFE, and SNXT only if SNEXT .GT.new SAFE *
C.    *         = 2  Compute both SAFE and SNXT                        *
C.    *         = 3  Compute SNXT only                                 *
C.    *       SNEXT (input)  : see IACT = 1                            *
C.    *       SNXT  (output) : distance to volume boundary             *
C.    *       SAFE  (output) : shortest distance to any boundary       *
C.    *                                                                *
C.    *    ==>Called by : GNEXT, GTNEXT                                *
C.    *         Author  A.McPherson,  P.Weidhaas  *********            *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gconsp.inc"
      DIMENSION X(6), PAR(9), PT(7), XT(6)
 
      EQUIVALENCE (PT(1),PT1), (PT(2),PT2), (PT(3),PT3), (PT(4),PT4)
      EQUIVALENCE (PT(5),PT5), (PT(6),PT6), (PT(7),PT7)
C.
C.   --------------------------------------------------------
C.
 
      SNXT = BIG
      SAFE = BIG
      R2 = X(1)*X(1) + X(2)*X(2)
      R  = SQRT (R2)
      NZ  = PAR(3)
      NSEC= NZ - 1
      PT6=PAR(1)
      PT7=PAR(1)+PAR(2)
 
      IFLAG = 2
      IF(PAR(2).EQ.360.0) IFLAG = 1
 
      IF (IACT .LT. 3) THEN
 
C       -------------------------------------------------
C       |  Compute safety-distance 'SAFE' (P.Weidhaas)  |
C       -------------------------------------------------
 
 
        SAFEZ  = 0.0
C
C......  Obtain axial distance "SAFEZ".
C
        ZMIN   = PAR(4)
        ZMAX   = PAR(3*NZ+1)
 
        IF (X(3) .LT. ZMIN) THEN
          SAFEZ  = ZMIN - X(3)
        ELSEIF (X(3) .GT. ZMAX) THEN
          SAFEZ  = X(3) - ZMAX
        ENDIF
 
 
C----------------------------------------------------
C......  Prepare parameters for PHI-segmented cone:
C----------------------------------------------------
 
      IF (IFLAG .EQ. 2) THEN
 
        PHI1 = MOD (PT6+360.0 , 360.0)
        PHI2 = MOD (PT7+360.0 , 360.0)
        IF ( X(1).NE.0.0 .OR. X(2).NE.0.0 ) THEN
           PHI = ATAN2( X(2), X(1) ) * RADDEG
        ELSE
           PHI = 0.0
        ENDIF
        PHI  = MOD (PHI+360.0 , 360.0)
 
        SINPH1 = SIN(PHI1*DEGRAD)
        COSPH1 = COS(PHI1*DEGRAD)
        SINPH2 = SIN(PHI2*DEGRAD)
        COSPH2 = COS(PHI2*DEGRAD)
 
 
C......  Set flag "IOPENP" if point (X(1),X(2),X(3)) lies in the
C......  PHI-opening.
 
        IOPENP = 0
        IF (PHI2 .GT. PHI1) THEN
          IF (PHI.GT.PHI2 .OR. PHI.LT.PHI1) IOPENP = 1
        ELSE
          IF (PHI.GT.PHI2 .AND. PHI.LT.PHI1) IOPENP = 1
        ENDIF
      ENDIF
 
 
 
C------------------   Start of loop over Z-sections  --------------
 
        IPZ2=4
 
        DO 150 IS=1,NSEC
 
        IPZ1=IPZ2
        IPZ2=IPZ1+3
 
        SAFSEG = 0.0
        SAFER  = 0.0
 
        XT3=X(3) - 0.5 * (PAR(IPZ1)+PAR(IPZ2))
        DZ   = 0.5 * (PAR(IPZ2)-PAR(IPZ1))
 
        PT2  = PAR(IPZ1+1)
        PT3  = PAR(IPZ1+2)
        PT4  = PAR(IPZ2+1)
        PT5  = PAR(IPZ2+2)
C**** check DZ=0 segments
        IF (DZ.LE.0.) THEN
           IF ((R-PT2)*(R-PT4).LE.0. .OR. (R-PT3)*(R-PT5).LE.0.) THEN
              SAFER = ABS(XT3)
           ELSE
              GO TO 150
           ENDIF
           GO TO 100
        ENDIF
 
        IF (PT2 .NE. PT4) GO TO 50
        IF (PT3 .NE. PT5) GO TO 50
 
C**********************************************************
C
C......  The segment is a tube;  invoke the algorithm
C......  from "GNOTUB" inline to get "SAFER" and "SAFSEG".
C
C**********************************************************
 
        IF (R .LT. PT2) SAFER = PT2 - R
        IF (R .GT. PT3) SAFER = R - PT3
 
        IF (IFLAG .EQ. 2)  THEN
 
C********************************************************************
C......  Handle the case in which we have a PHI-segment of a tube.
C......  In addition to the radial distance (SAFER) and the
C......  axial distance (SAFEZ) we compute here the distance (SAFSEG)
C......  to the PHI-segment boundary that is closest to the point:
C
C......  SAFSEG is only calculated if PHI lies outside the interval
C......  [PHI1, PHI2]. Here PHI is the angle to the given point
C......  (thus we only consider SAFSEG if the point is outside the
C......  PHI-segment).
C
C......  Algorithm to find SAFSEG (same as in routine "GNTUBE"):
C
C......  For each PHI-boundary we find the distance from the given
C......  point to the outer (at RMAX) point of the segment boundary
C......  (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
C......  SAFSEG to be the distance to segment PHI1; else we set
C......  SAFSEG to be the distance to segment PHI2.
C********************************************************************
 
 
C......  Next eliminate those points whose angle PHI places them
C......  inside the given PHI-segment (IOPENP = 0).
 
          IF (IOPENP .EQ. 0) GO TO 100
 
 
C......  Get coordinates of outer endpoints (at RMAX) of both PHI-segments.
 
          XS1  = PT3 * COSPH1
          YS1  = PT3 * SINPH1
          XS2  = PT3 * COSPH2
          YS2  = PT3 * SINPH2
 
C......  Get distances (squared) from the given point to each endpoint.
 
          DISTS1 = (X(1) - XS1)**2  +  (X(2) - YS1)**2
          DISTS2 = (X(1) - XS2)**2  +  (X(2) - YS2)**2
 
C......  Get distance to that PHI-segment whose endpoint
C......  is closest to the given point.
 
          IF (DISTS1 .LE. DISTS2) THEN
            SAFSEG = ABS(X(1) * SINPH1 - X(2) * COSPH1)
          ELSE
            SAFSEG = ABS(X(1) * SINPH2 - X(2) * COSPH2)
          ENDIF
 
        ENDIF
 
        GO TO 100
 
 
   50   CONTINUE
 
C*********************************************************
C
C......  The segment is a cone; invoke the algorithm
C......  from "GNOCON" inline to get "SAFER" and "SAFSEG".
C
C*********************************************************
 
        ZLENI = 0.5 / DZ
 
        FACT1  = (PT4 - PT2) * ZLENI
        FACT2  = (PT5 - PT3) * ZLENI
        RIN   = PT2 + FACT1 * (DZ + XT3)
        ROUT  = PT3 + FACT2 * (DZ + XT3)
 
        IF (R .LT. RIN) THEN
          SAFER = (RIN - R) / SQRT(1.0 + FACT1 * FACT1)
        ELSE
          IF (R .GT. ROUT) THEN
            SAFER = (R - ROUT) / SQRT(1.0 + FACT2 * FACT2)
          ENDIF
        ENDIF
 
 
        IF (IFLAG .EQ. 2)  THEN
 
C********************************************************************
C......  Handle the case in which we have a PHI-segment of a cone.
C......  In addition to the radial distance (SAFER) and the
C......  axial distance (SAFEZ) we compute here the distance (SAFSEG)
C......  to the PHI-segment boundary that is closest to the point:
C
C......  SAFSEG is only calculated if PHI lies outside the interval
C......  [PHI1, PHI2]. Here PHI is the angle to the given point
C......  (thus we only consider SAFSEG if the point is outside the
C......  PHI-segment).
C
C......  Algorithm to find SAFSEG (same as in routine "GNTUBE"):
C
C......  For each PHI-boundary we find the distance from the given
C......  point to the outer (at ROUT) point of the segment boundary
C......  (DISTS1 and DISTS2, resp.). If DISTS1 < DISTS2, we define
C......  SAFSEG to be the distance to segment PHI1; else we set
C......  SAFSEG to be the distance to segment PHI2.
C********************************************************************
 
 
C......  Next eliminate those points whose angle PHI places them
C......  inside the given PHI-segment (IOPENP = 0).
 
           IF (IOPENP .EQ. 0) GO TO 100
 
C......  Get coordinates of outer endpoints (at ROUT) of both PHI-segments.
 
          IF (XT3 .LT. -DZ) THEN
            ROUT = PT3
          ELSEIF (XT3 .GT. DZ) THEN
            ROUT = PT5
          ENDIF
 
          XS1  = ROUT * COSPH1
          YS1  = ROUT * SINPH1
          XS2  = ROUT * COSPH2
          YS2  = ROUT * SINPH2
 
C......  Get distances (squared) from the given point to each endpoint.
 
          DISTS1 = (X(1) - XS1)**2  +  (X(2) - YS1)**2
          DISTS2 = (X(1) - XS2)**2  +  (X(2) - YS2)**2
 
C......  Obtain distance to that PHI-segment whose endpoint
C......  is closest to the given point.
 
          IF (DISTS1 .LE. DISTS2) THEN
            SAFSEG = ABS(X(1) * SINPH1 - X(2) * COSPH1)
          ELSE
            SAFSEG = ABS(X(1) * SINPH2 - X(2) * COSPH2)
          ENDIF
 
        ENDIF
 
 
  100   CONTINUE
        TSAFE  = MAX (SAFEZ, SAFER, SAFSEG)
 
        IF (TSAFE .GT. 0.0) THEN
          IF (TSAFE .LT. SAFE) SAFE  = TSAFE
        ENDIF
 
        IF (TSAFE .EQ. 0.0) THEN
          IF (X(3) .LT. PAR(IPZ1)) GO TO 200
        ENDIF
  150   CONTINUE
 
  200   CONTINUE
        IF (IACT .EQ. 0) GO TO 999
        IF (IACT .EQ. 1) THEN
          IF (SNEXT .LT. SAFE) GO TO 999
        ENDIF
      ENDIF
 
 
C     ------------------------------------------------
C     |  Compute vector-distance 'SNXT' (McPherson)  |
C     ------------------------------------------------
 
      TSNXT = BIG
 
      DO 210 I=1, 6
        XT(I) = X(I)
  210 CONTINUE
 
      IPZ2 = 4
 
      DO 300  IS=1, NSEC
        IPZ1  = IPZ2
        IPZ2  = IPZ1 + 3
 
        XT(3) = X(3) - 0.5 * (PAR(IPZ1) + PAR(IPZ2))
 
        PT1 = 0.5 * (PAR(IPZ2) - PAR(IPZ1))
        IF (PT1 .LE. 0.0) GO TO 300
 
        IF (PAR(IPZ1+1) .NE. PAR(IPZ2+1)) GO TO 250
        IF (PAR(IPZ1+2) .NE. PAR(IPZ2+2)) GO TO 250
 
C
C......  This Z-section is a tube.
C
        PT3 = PT1
        PT1 = PAR(IPZ1+1)
        PT2 = PAR(IPZ1+2)
        PT4 = PT6
        PT5 = PT7
 
        CALL GNOTUB (XT, PT, 3, IFLAG, SNEXT, TSNXT, TSAFE)
        GO TO 280
 
 
  250   CONTINUE
C
C......   This Z-section is a cone.
C
        PT2 = PAR(IPZ1+1)
        PT3 = PAR(IPZ1+2)
        PT4 = PAR(IPZ2+1)
        PT5 = PAR(IPZ2+2)
 
        CALL GNOCON (XT, PT, 3, IFLAG, SNEXT, TSNXT, TSAFE)
 
  280 CONTINUE
      IF (TSNXT .LT. SNXT)  SNXT = TSNXT
  300 CONTINUE
 
  999 CONTINUE
      END
