***   
*     subroutine to compute the Integral(B dl) along the trajectory
***   


***   
*     
*     created by paolo papini, 11 oct 2005
*     
*     c$$$ max:  modified by massimo bongi, 11 oct 2005 (later...)
*     
***


      subroutine CalcBdL(Nstep,BdL,IFAIL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      real BdL,vvv(3),fff(3)
      real*8 f(3)
      include 'commontracker.f' !tracker general common
      include 'common_mini_2.f' !common for the tracking procedure
*     
*     set parameters for GRKUTA from ZINI to the first plane
*     
c$$$  max    IF(AL(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
c$$$  max    IF(AL(5).NE.0) CHARGE=AL_P(5)/DABS(AL_P(5))
c$$$      IF(AL(5).NE.0) CHARGE=AL(5)/DABS(AL(5))
      IF(AL(5).NE.0) CHARGE=AL(5)/DABS(AL(5))
      IF(AL(5).EQ.0) CHARGE=1.
      VOUT(1)=AL(1)
      VOUT(2)=AL(2)
      VOUT(3)=ZINI              !  DBLE(Z0)-DBLE(ZSPEC)
      VOUT(4)=AL(3)*DCOS(AL(4))
      VOUT(5)=AL(3)*DSIN(AL(4))
      VOUT(6)=-1.*DSQRT(1.-AL(3)**2)
      IF(AL(5).NE.0.) VOUT(7)=DABS(1./AL(5))
      IF(AL(5).EQ.0.) VOUT(7)=1.E8
*     
*     first integrate the track from reference plane to first trk plane
*
*      zout=zm(1)     
      zout=zv(1)
      step=vout(3)-zout
 10   DO J=1,7
         VECT(J)=VOUT(J)
         VECTINI(J)=VOUT(J)
      ENDDO
 11   continue
      CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
      IF(VOUT(3).GT.VECT(3)) THEN
         IFAIL=1
         PRINT *,'BdL (grkuta): WARNING ===> backward track!!'
         print*,'charge',charge
         print*,'vect',vect
         print*,'vout',vout
         print*,'step',step
         RETURN
      ENDIF
      Z=VOUT(3)
      IF(Z.LE.Zout+TOLL.AND.Z.GE.Zout-TOLL) GOTO 100
      IF(Z.GT.Zout+TOLL) GOTO 10      
      IF(Z.LE.Zout-TOLL) THEN
         STEP=STEP*(Zout-VECT(3))/(Z-VECT(3))
         DO J=1,7
            VECT(J)=VECTINI(J)
         ENDDO
         GOTO 11
      ENDIF
*
*     then integrate B o dL along the track, from the first to the last plane
*
 100  continue
      BdL=0.
*      dz=(zm(1)-zm(6))/float(Nstep)
      dz = (zv(1)-zv(6))/float(Nstep)
      zout = zv(6)
      
      do i=1,10*Nstep
         BdLstep=0.
         vvv(1)=sngl(vout(1))
         vvv(2)=sngl(vout(2))
         vvv(3)=sngl(vout(3))
         CALL GUFLD(VVV,FFF)
         do j=1,3
            f(j)=dble(fff(j))
         enddo
c         zout=vout(3)-dz
         step = dz
 20      DO J=1,7
            VECT(J)=VOUT(J)
            VECTINI(J)=VOUT(J)
         ENDDO
 21      continue
         CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
c$$$         BdLstep=BdLstep+
c$$$     $        (vout(1)-vect(1))*f(1)+
c$$$     $        (vout(2)-vect(2))*f(2)+
c$$$     $        (vout(3)-vect(3))*f(3)
         BdLstep=BdLstep+
     $        sqrt(
     $        ((vout(2)-vect(2))*f(3)-(vout(3)-vect(3))*f(2))**2+
     $        ((vout(3)-vect(3))*f(1)-(vout(1)-vect(1))*f(3))**2+
     $        ((vout(1)-vect(1))*f(2)-(vout(2)-vect(2))*f(1))**2
     $        )
         IF(VOUT(3).GT.VECT(3)) THEN
            IFAIL=1
            PRINT *,'BdL (grkuta): WARNING ===> backward track!!'
            print*,'charge',charge
            print*,'vect',vect
            print*,'vout',vout
            print*,'step',step
            RETURN
         ENDIF
 200     continue
         BdL=BdL+sngl(BdLstep*1.D-3) !Tesla * meters
c         print*,i,' #### ',BdL

         Z=VOUT(3)
         IF(Z.LE.Zout+TOLL.AND.Z.GE.Zout-TOLL) GOTO 2000 !end
         IF(Z.LE.Zout-TOLL) THEN
            STEP=STEP*(Zout-VECT(3))/(Z-VECT(3))
            DO J=1,7
               VECT(J)=VECTINI(J)
            ENDDO
            BdL=BdL-sngl(BdLstep*1.D-3) !Tesla * meters            
            BdLstep=0.
            GOTO 21
         ENDIF
      enddo
 2000 continue
c$$$      print*,' #### ',BdL
c$$$      print*,' #### ',xgood
c$$$      print*,' #### ',ygood
c$$$      print*,'========================================='
      return
      end