*************************************************************************
*     
*     Subroutine inter_B.f (from tracker software analysis)
*     
*     it computes the magnetic field in a chosen point x,y,z inside or 
*     outside the magnetic cavity, using a trilinear interpolation of
*     B field measurements (read before by means of ./read_B.f)
*     if the point falls outside the interpolation volume, set the field to 0
*     
*     needs:
*     - ../common/common_B_inner.f common file for the inner magnetic field map
*     - ./inter_B_inner.f common file for the inner magnetic field map
*     
*     to be called after ./read_B.f (magnetic field map reading subroutine)
*     
*     input: coordinates in m
*     output: magnetic field in T
*     
*************************************************************************

      subroutine inter_B(x,y,z,res) !coordinates in m, magnetic field in T
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "gpfield.inc"

c------------------------------------------------------------------------
c     
c     local variables
c     
c------------------------------------------------------------------------

      real*8 x,y,z              !point of interest
      real*8 res(3)             !interpolated B components: res = (Bx, By, Bz)

      real*8 zl,zu
      real*8 resu(3),resl(3)
c------------------------------------------------------------------------
c     
c     set the field outside the interpolation volume to be 0
c     
c------------------------------------------------------------------------

      do ip=1,3
        res(ip)=0.
      enddo

c------------------------------------------------------------------------
c     
c     check if the point falls inside the interpolation volumes
c     
c------------------------------------------------------------------------

*     -----------------------
*     INNER MAP
*     -----------------------
      if(
     $          (x.ge.edgexmin).and.(x.le.edgexmax)
     $     .and.(y.ge.edgeymin).and.(y.le.edgeymax)
     $     .and.(z.ge.edgezmin).and.(z.le.edgezmax)
     $     ) then


        call inter_B_inner(x,y,z,res)
c        print*,'INNER - ',z,res

      endif

*     -----------------------
*     OUTER MAP
*     -----------------------
      if(
     $         ((x.ge.edgeuxmin).and.(x.le.edgeuxmax)
     $     .and.(y.ge.edgeuymin).and.(y.le.edgeuymax)
     $     .and.(z.ge.edgeuzmin).and.(z.le.edgeuzmax))
     $     .or.
     $         ((x.ge.edgelxmin).and.(x.le.edgelxmax)
     $     .and.(y.ge.edgelymin).and.(y.le.edgelymax)
     $     .and.(z.ge.edgelzmin).and.(z.le.edgelzmax))
     $     ) then
         

        call inter_B_outer(x,y,z,res)
c        res(2)=res(2)*10
c        print*,'OUTER - ',z,res

      endif

*     --------------------------------
*     GAP between INNER and OUTER MAPS
*     --------------------------------
      if(
     $          (x.gt.edgexmin).and.(x.lt.edgexmax)
     $     .and.(y.gt.edgeymin).and.(y.lt.edgeymax)
     $     .and.(z.gt.edgezmax).and.(z.lt.edgeuzmin)
     $     )then
         
         zu = edgeuzmin
         zl = edgezmax
         call inter_B_inner(x,y,zl,resu)
         call inter_B_outer(x,y,zu,resl)

         do i=1,3
            res(i) = z  * ((resu(i)-resl(i))/(zu-zl))
     $           + resu(i)
     $           -   zu * ((resu(i)-resl(i))/(zu-zl))
         enddo
c        print*,'GAP U - ',z,res

      elseif(
     $          (x.gt.edgexmin).and.(x.lt.edgexmax)
     $     .and.(y.gt.edgeymin).and.(y.lt.edgeymax)
     $     .and.(z.gt.edgelzmax).and.(z.lt.edgezmin)
     $     ) then
         
        
         zu = edgezmin
         zl = edgelzmax
         call inter_B_inner(x,y,zu,resu)
         call inter_B_outer(x,y,zl,resl)

         
         do i=1,3
            res(i) = z  * ((resu(i)-resl(i))/(zu-zl))
     $           + resu(i)
     $           -   zu * ((resu(i)-resl(i))/(zu-zl))
         enddo
c        print*,'GAP D - ',z,res

      endif


      return
      end