**********************************************************************
*
*
*     routine per tracciare la particella di uno STEP 
*
      SUBROUTINE GRKUTA (CHARGE,STEP,VECT,VOUT)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *  Runge-Kutta method for tracking a particle through a magnetic *
C.    *  field. Uses Nystroem algorithm (See Handbook Nat. Bur. of     *
C.    *  Standards, procedure 25.5.20)                                 *
C.    *                                                                *
C.    *  Input parameters                                              *
C.    *       CHARGE    Particle charge                                *
C.    *       STEP      Step size                                      *
C.    *       VECT      Initial co-ords,direction cosines,momentum     *
C.    *  Output parameters                                             *
C.    *       VOUT      Output co-ords,direction cosines,momentum      *
C.    *  User routine called                                           *
C.    *       CALL GUFLD(X,F)                                          *
C.    *                                                                *
C.    *    ==>Called by : <USER>, GUSWIM                               *
C.    *       Authors    R.Brun, M.Hansroul  *********                 *
C.    *                  V.Perevoztchikov (CUT STEP implementation)    *
C.    *                                                                *
C.    *                                                                *
C.    ******************************************************************
C.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
*
      REAL VVV(3),FFF(3)
      REAL*8 CHARGE, STEP, VECT(*), VOUT(*), F(4)
      REAL*8 XYZT(3), XYZ(3), X, Y, Z, XT, YT, ZT
      DIMENSION SECXS(4),SECYS(4),SECZS(4),HXP(3)
      EQUIVALENCE (X,XYZ(1)),(Y,XYZ(2)),(Z,XYZ(3)),
     +            (XT,XYZT(1)),(YT,XYZT(2)),(ZT,XYZT(3))
*
      PARAMETER (MAXIT = 1992, MAXCUT = 11)
      PARAMETER (EC=2.9979251D-4,DLT=1D-4,DLT32=DLT/32)
      PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3)
      PARAMETER (THIRD=ONE/THREE, HALF=ONE/TWO)
      PARAMETER (PISQUA=.986960440109D+01)
      PARAMETER (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6)

*.
*.    ------------------------------------------------------------------
*.
*             This constant is for units CM,GEV/C and KGAUSS
*      
c      print *,'GRKUTA+++++ charge ',charge,' step ',step,' vect(3) ',
c     &     vect(3)
      ITER = 0
      NCUT = 0
      DO 10 J=1,7
         VOUT(J)=VECT(J)
c         print *,' grkuta j ',j,' vout ',vout(j)
c         print *,' grkuta j ',j,' vect ',vect(j)
   10 CONTINUE
      PINV   = EC * CHARGE / VECT(7)
      TL = 0.
      H      = STEP
*
*
   20 REST  = STEP-TL
      IF (DABS(H).GT.DABS(REST)) H = REST
      DO I=1,3
       VVV(I)=SNGL(VOUT(I))
c       print *,'grkuta i ',i,' vvv ',vvv(i)
      ENDDO
c      print *,'grkuta pinv ',pinv,' h ',h,' rest ',rest
      
      CALL GUFLD(VVV,FFF)
      DO I=1,3
       F(I)=DBLE(FFF(I))
c       print *,'grkuta i ',i,' f ',f(i)
      ENDDO
*
*             Start of integration
*
      X      = VOUT(1)
      Y      = VOUT(2)
      Z      = VOUT(3)
      A      = VOUT(4)
      B      = VOUT(5)
      C      = VOUT(6)
c      print *,' QUI A ',A,' B ',B,' C ',C
c      print *,' QUI x ',x,' y ',y,' z ',z
*
      H2     = HALF * H
      H4     = HALF * H2
      PH     = PINV * H
      PH2    = HALF * PH
      SECXS(1) = (B * F(3) - C * F(2)) * PH2
      SECYS(1) = (C * F(1) - A * F(3)) * PH2
      SECZS(1) = (A * F(2) - B * F(1)) * PH2
      ANG2 = (SECXS(1)**2 + SECYS(1)**2 + SECZS(1)**2)
      IF (ANG2.GT.PISQUA) GO TO 40
      DXT    = H2 * A + H4 * SECXS(1)
      DYT    = H2 * B + H4 * SECYS(1)
      DZT    = H2 * C + H4 * SECZS(1)
      XT     = X + DXT
      YT     = Y + DYT
      ZT     = Z + DZT
*
*              Second intermediate point
*
      EST = DABS(DXT)+DABS(DYT)+DABS(DZT)
      IF (EST.GT.H) GO TO 30

      DO I=1,3
       VVV(I)=SNGL(XYZT(I))
      ENDDO
      CALL GUFLD(VVV,FFF)
      DO I=1,3
       F(I)=DBLE(FFF(I))
c       print *,'2grkuta i ',i,' f ',f(i)
      ENDDO     
C      CALL GUFLD(XYZT,F)
      AT     = A + SECXS(1)
      BT     = B + SECYS(1)
      CT     = C + SECZS(1)
*
      SECXS(2) = (BT * F(3) - CT * F(2)) * PH2
      SECYS(2) = (CT * F(1) - AT * F(3)) * PH2
      SECZS(2) = (AT * F(2) - BT * F(1)) * PH2
c      print *,'1 at ',xt,' bt ',yt,' ct ',zt
      AT     = A + SECXS(2)
      BT     = B + SECYS(2)
      CT     = C + SECZS(2)
c      print *,'2 at ',xt,' bt ',yt,' ct ',zt
      SECXS(3) = (BT * F(3) - CT * F(2)) * PH2
      SECYS(3) = (CT * F(1) - AT * F(3)) * PH2
      SECZS(3) = (AT * F(2) - BT * F(1)) * PH2
      DXT    = H * (A + SECXS(3))
      DYT    = H * (B + SECYS(3))
      DZT    = H * (C + SECZS(3))
      XT     = X + DXT
      YT     = Y + DYT
      ZT     = Z + DZT
c      print *,' xt ',xt,' yt ',yt,' zt ',zt
c      print *,' dxt ',xt,' dyt ',yt,' dzt ',zt
c      print *,' at ',xt,' bt ',yt,' ct ',zt
      AT     = A + TWO*SECXS(3)
      BT     = B + TWO*SECYS(3)
      CT     = C + TWO*SECZS(3)
*
      EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
      IF (EST.GT.2.*ABS(H)) GO TO 30

      DO I=1,3
       VVV(I)=SNGL(XYZT(I))
c       print *,'3grkuta i ',i,' vvv ',vvv(i)
c       print *,'3grkuta i ',i,' xyzt ',xyzt(i)
      ENDDO
      CALL GUFLD(VVV,FFF)
      DO I=1,3
       F(I)=DBLE(FFF(I))
c       print *,'3grkuta i ',i,' f ',f(i)
c       print *,'3grkuta i ',i,' fff ',fff(i)
      ENDDO
C      CALL GUFLD(XYZT,F)
*
      Z      = Z + (C + (SECZS(1) + SECZS(2) + SECZS(3)) * THIRD) * H
      Y      = Y + (B + (SECYS(1) + SECYS(2) + SECYS(3)) * THIRD) * H
      X      = X + (A + (SECXS(1) + SECXS(2) + SECXS(3)) * THIRD) * H
*
      SECXS(4) = (BT*F(3) - CT*F(2))* PH2
c      print *,'secxs4 bt ',bt,' ct ',ct,' ph2 ',ph2
      SECYS(4) = (CT*F(1) - AT*F(3))* PH2
      SECZS(4) = (AT*F(2) - BT*F(1))* PH2      
      A      = A+(SECXS(1)+SECXS(4)+TWO * (SECXS(2)+SECXS(3))) * THIRD
c      print *,'=> a ',a,' secxs 1 ',SECXS(1),' 4 ',SECXS(4),' 2 ',
c     &     SECXS(2),' 3 ',SECXS(3),' third ',third
      B      = B+(SECYS(1)+SECYS(4)+TWO * (SECYS(2)+SECYS(3))) * THIRD
      C      = C+(SECZS(1)+SECZS(4)+TWO * (SECZS(2)+SECZS(3))) * THIRD
*
      EST    = ABS(SECXS(1)+SECXS(4) - (SECXS(2)+SECXS(3)))
     ++        ABS(SECYS(1)+SECYS(4) - (SECYS(2)+SECYS(3)))
     ++        ABS(SECZS(1)+SECZS(4) - (SECZS(2)+SECZS(3)))
*
      IF (EST.GT.DLT .AND. ABS(H).GT.1.E-4) GO TO 30
      ITER = ITER + 1
      NCUT = 0
*               If too many iterations, go to HELIX
      IF (ITER.GT.MAXIT) GO TO 40
*
      TL = TL + H
      IF (EST.LT.(DLT32)) THEN
         H = H*TWO
      ENDIF
      CBA    = ONE/ SQRT(A*A + B*B + C*C)
      VOUT(1) = X
      VOUT(2) = Y
      VOUT(3) = Z
      VOUT(4) = CBA*A
      VOUT(5) = CBA*B
      VOUT(6) = CBA*C
      REST = STEP - TL
      IF (STEP.LT.0.) REST = -REST
      IF (REST .GT. 1.E-5*DABS(STEP)) GO TO 20
*
c      print *,' x ',x,' y ',y,' z ',z,' cba ',cba,
c     &     ' a ',a,' b ',b,' c ',c,' step ',step,
c     &     ' tl ',tl,' rest ',rest,' est ',est,
c     &     ' h ',h,' two ',two,' one ',one
      GO TO 999
*
**              CUT STEP
   30 NCUT = NCUT + 1
*               If too many cuts , go to HELIX
      IF (NCUT.GT.MAXCUT)       GO TO 40
      H = H*HALF
      GO TO 20
*
**              ANGLE TOO BIG, USE HELIX
   40 F1  = F(1)
      F2  = F(2)
      F3  = F(3)
      F4  = DSQRT(F1**2+F2**2+F3**2)
      RHO = -F4*PINV
      TET = RHO * STEP
      IF(TET.NE.0.) THEN
         HNORM = ONE/F4
         F1 = F1*HNORM
         F2 = F2*HNORM
         F3 = F3*HNORM
*
         HXP(1) = F2*VECT(IPZ) - F3*VECT(IPY)
         HXP(2) = F3*VECT(IPX) - F1*VECT(IPZ)
         HXP(3) = F1*VECT(IPY) - F2*VECT(IPX)

         HP = F1*VECT(IPX) + F2*VECT(IPY) + F3*VECT(IPZ)
*
         RHO1 = ONE/RHO
         SINT = DSIN(TET)
         COST = TWO*DSIN(HALF*TET)**2
*
         G1 = SINT*RHO1
         G2 = COST*RHO1
         G3 = (TET-SINT) * HP*RHO1
         G4 = -COST
         G5 = SINT
         G6 = COST * HP

         VOUT(IX) = VECT(IX) + (G1*VECT(IPX) + G2*HXP(1) + G3*F1)
         VOUT(IY) = VECT(IY) + (G1*VECT(IPY) + G2*HXP(2) + G3*F2)
         VOUT(IZ) = VECT(IZ) + (G1*VECT(IPZ) + G2*HXP(3) + G3*F3)
c         print *,' iz ',iz,' vect ',vect(iz),' g1 ',g1,' ipz ',ipz,
c     &        ' vect ',vect(ipz),' g2 ',g2,' hxp ',hxp(3),' g3 ',g3,
c     &        ' f3 ',f3

         VOUT(IPX) = VECT(IPX) + (G4*VECT(IPX) + G5*HXP(1) + G6*F1)
         VOUT(IPY) = VECT(IPY) + (G4*VECT(IPY) + G5*HXP(2) + G6*F2)
         VOUT(IPZ) = VECT(IPZ) + (G4*VECT(IPZ) + G5*HXP(3) + G6*F3)
*
      ELSE
         VOUT(IX) = VECT(IX) + STEP*VECT(IPX)
         VOUT(IY) = VECT(IY) + STEP*VECT(IPY)
         VOUT(IZ) = VECT(IZ) + STEP*VECT(IPZ)
c         print *,' iz ',iz,' vect ',vect(iz),' step ',step,' ipz ',ipz,
c     &        ' vect ',vect(ipz)
*     
      ENDIF
*
  999 END
*     
*     



**********************************************************************
*     
*     gives the value of the magnetic field in the tracking point
*     
**********************************************************************

      subroutine  gufld(v,f)    !coordinates in cm, B field in kGauss

      real v(3),f(3)            !coordinates in cm, B field in kGauss, error in kGauss

      real*8 vv(3),ff(3)        !inter_B.f works in double precision


        do i=1,3
          vv(i)=v(i)/100.       !inter_B.f works in meters
c          print *,'IN gufld i ',i,' v ',v(i)
        enddo
c       inter_B: coordinates in m, B field in Tesla
        call inter_B(vv(1),vv(2),vv(3),ff) 
        do i=1,3                !change back the field in kGauss 
          f(i)=ff(i)*10.
        enddo
c        do i=1,3
c          print *,'OUT gufld i ',i,' v ',v(i)
c        enddo
      return
      end

