*************************************************************
*
*     Routine to compute the NPOINT track intersection points
*     with planes of z-coordinate given by ZIN
*     given the track parameters
*
*     The routine is based on GRKUTA, which computes the 
*     trajectory of a charged particle in a magnetic field
*     by solving the equatoins of motion with Runge-Kuta method
*
*     Variables that have to be assigned when the subroutine 
*     is called:
*
*     ZIN(1-NPOINT) ----> z coordinates of the planes
*     AL_P(1-5)     ----> track-parameter vector
*
*     NB !!! 
*     The routine works properly only if the 
*     planes are numbered in descending order starting from the
*     reference plane (ZINI)
*     
**************************************************************

      SUBROUTINE TRACK(NPOINT,ZIN,XOUT,YOUT,AL_P,IFAIL)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      DIMENSION VECT(7),VECTINI(7),VOUT(7)
      DATA TOLL/1.d-8/          
*     tolerance in reaching the next plane during the tracking procedure
*     -----------------------------------------------
*     I/O parameters
      PARAMETER (NPOINT_MAX=100)
      DIMENSION ZIN(NPOINT_MAX)
      DIMENSION XOUT(NPOINT_MAX)
      DIMENSION YOUT(NPOINT_MAX)
      DIMENSION AL_P(5)
      DOUBLE PRECISION Z
      DOUBLE PRECISION STEP
      DOUBLE PRECISION CHARGE
*     -----------------------------------------------      
      DATA ZINI/23.5/           !z coordinate of the reference plane
      
*     ==================================================================
*     divide the track in two parts: below and above the reference plane
*     ==================================================================      
      Z = 0.
      IUPDOWN=0
      DO I=1,NPOINT   
         IF(ZIN(I).LT.ZINI)THEN
            if(i.ne.1)IUPDOWN=I
            GOTO 88
         ENDIF
         IUPDOWN=NPOINT+1
      ENDDO
 88   CONTINUE
c      DO I=1,NPOINT   
c         print *,' TRACK >>>>>>>>> ',I,' ZIN  ',ZIN(i)
c      enddo
c      print *,' iupdown = ',iupdown


*     ==================================================================
*     track from reference plane DOWN
*     ==================================================================
*     parameters for GRKUTA
      IF(AL_P(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
      IF(AL_P(5).EQ.0) CHARGE=1.
      VOUT(1)=AL_P(1)
      VOUT(2)=AL_P(2)
      VOUT(3)=ZINI             
      VOUT(4)=AL_P(3)*DCOS(AL_P(4))
      VOUT(5)=AL_P(3)*DSIN(AL_P(4))
      VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
      IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
      IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
c
      DO I=MAX(IUPDOWN,1),NPOINT 
         step=vout(3)-zin(i)
c         print*,'DOWN ',i,' - Track from ',
c     $        vout(3),' to ',zin(i),' step ',step
 10      DO J=1,7
            VECT(J)=VOUT(J)
            VECTINI(J)=VOUT(J)
c            print *,'vect ',j,' ',vect(j)
         ENDDO
 11      continue
         CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
         IF(VOUT(3).GT.VECT(3)) THEN
            IFAIL=1
            PRINT *,'=== WARNING ===> inverted tracking (DOWN)'
            print*,'charge',charge
            print*,'vect',vect
            print*,'vout',vout
            print*,'step',step
            RETURN
         ENDIF
         Z=VOUT(3)
c         if(Z.ne.Z) goto 11
         IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 100
         IF(Z.GT.ZIN(I)+TOLL) GOTO 10      
         IF(Z.LE.ZIN(I)-TOLL) THEN
            STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
            DO J=1,7
               VECT(J)=VECTINI(J)
            ENDDO
            GOTO 11
         ENDIF
 100     XOUT(I)=VOUT(1)
         YOUT(I)=VOUT(2)
         ZIN(I)=VOUT(3)
c         print *,'vout3 ',VOUT(3)
      ENDDO

*     ==================================================================
*     track from refernce plane UP
*     ==================================================================
*     parameters for GRKUTA:
*     -opposite charge
*     -opposite momentum direction
      IF(AL_P(5).NE.0) CHARGE=-AL_P(5)/DABS(AL_P(5))
      IF(AL_P(5).EQ.0) CHARGE=-1.
      VOUT(1)=AL_P(1)
      VOUT(2)=AL_P(2)
      VOUT(3)=ZINI             
      VOUT(4)=-AL_P(3)*DCOS(AL_P(4))
      VOUT(5)=-AL_P(3)*DSIN(AL_P(4))
      VOUT(6)=1.*DSQRT(1.-AL_P(3)**2)
      IF(AL_P(5).NE.0.) VOUT(7)=DABS(1./AL_P(5))
      IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
      DO I=MIN((IUPDOWN-1),NPOINT),1,-1
         step=vout(3)-zin(i)
         step = -step
c         print*,'UP ',i,' - Track from ',
c     $        vout(3),' to ',zin(i),' step ',step
 20      DO J=1,7
            VECT(J)=VOUT(J)
            VECTINI(J)=VOUT(J)
         ENDDO
 22      continue
         CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
         IF(VOUT(3).LT.VECT(3)) THEN
            IFAIL=1
            PRINT *,'=== WARNING ===> inverted tracking (UP)'
            print*,'charge',charge
            print*,'vect',vect
            print*,'vout',vout
            print*,'step',step
            RETURN
         ENDIF
         Z=VOUT(3)
         IF(Z.LE.ZIN(I)+TOLL.AND.Z.GE.ZIN(I)-TOLL) GOTO 200
         IF(Z.LT.ZIN(I)-TOLL) GOTO 20      
         IF(Z.GE.ZIN(I)+TOLL) THEN
            STEP=STEP*(ZIN(I)-VECT(3))/(Z-VECT(3))
            DO J=1,7
               VECT(J)=VECTINI(J)
            ENDDO
            GOTO 22
         ENDIF
 200     XOUT(I)=VOUT(1)
         YOUT(I)=VOUT(2)
         ZIN(I)=VOUT(3)
      ENDDO

      RETURN
      END

      include './grkuta.for'
      include './inter_B.for'