*************************************************************************
*     
*     Subroutine inter_B.f
*     
*     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_B_inner.f common file for the inner magnetic field map
*     - ./inter_B_inner.f common file for the inner magnetic field map
*     - common_B_outer.f common file for the outer magnetic field map
*     - ./inter_B_outer.f common file for the outer 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 'common_B.f'


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


*************************************************************************
*     
*     Subroutine inter_B_inner.f
*     
*     it computes the magnetic field in a chosen point x,y,z inside the
*     magnetic cavity, using a trilinear interpolation of
*     B field measurements (read before by means of ./read_B.f)
*     the value is computed for two different inner maps and then averaged
*     
*     needs:
*     - ../common/common_B_inner.f
*     
*     input: coordinates in m
*     output: magnetic field in T
*     
*************************************************************************

      subroutine inter_B_inner(x,y,z,res) !coordinates in m, magnetic field in T

      implicit double precision (a-h,o-z)
      include 'common_B.f'


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

      real*8 x,y,z              !point of interpolation
      real*8 res(3)             !interpolated B components: res = (Bx, By, Bz)
      real*8 res1(3),res2(3)    !interpolated B components for the two maps

      integer ic                !index for B components:
                                ! ic=1 ---> Bx
                                ! ic=2 ---> By
                                ! ic=3 ---> Bz

      integer cube(3)           !vector of indexes identifying the cube
                                ! containing the point of interpolation
                                ! (see later...)
      
      real*8 xl,xh,yl,yh,zl,zh  !cube vertexes coordinates

      real*8 xr,yr,zr           !reduced variables (coordinates of the 
                                ! point of interpolation inside the cube)

      real*8 Bp(8)              !vector of values of B component 
                                ! being computed, on the eight cube vertexes


c------------------------------------------------------------------------
c     
c     *** FIRST MAP ***
c     
c------------------------------------------------------------------------

      do ic=1,3                 !loops on the three B components

c------------------------------------------------------------------------
c     
c     chooses the coordinates interval containing the input point
c     
c------------------------------------------------------------------------
c     e.g.:
c     
c     x1    x2    x3    x4    x5...
c     |-----|-+---|-----|-----|----
c     ~~~~~~~~x
c     
c     in this case the right interval is identified by indexes 2-3, so the 
c     value assigned to cube variable is 2

        cube(1)=INT((nx-1)*(x-px1min(ic))/(px1max(ic)-px1min(ic))) +1
        cube(2)=INT((ny-1)*(y-py1min(ic))/(py1max(ic)-py1min(ic))) +1
        cube(3)=INT((nz-1)*(z-pz1min(ic))/(pz1max(ic)-pz1min(ic))) +1
        
c------------------------------------------------------------------------
c     
c     if the point falls beyond the extremes of the grid...
c     
c------------------------------------------------------------------------
c     
c     ~~~~~~~~~~x1    x2    x3...
c     - - + - - |-----|-----|----
c     ~~~~x
c     
c     in the case cube = 1
c     
c     
c     ...nx-2  nx-1  nx
c     ----|-----|-----| - - - + - -
c     ~~~~~~~~~~~~~~~~~~~~~~~~x
c     
c     in this case cube = nx-1

        if (cube(1).le.0) cube(1) = 1
        if (cube(2).le.0) cube(2) = 1
        if (cube(3).le.0) cube(3) = 1
        if (cube(1).ge.nx) cube(1) = nx-1
        if (cube(2).ge.ny) cube(2) = ny-1
        if (cube(3).ge.nz) cube(3) = nz-1


c------------------------------------------------------------------------
c     
c     temporary variables definition for field computation
c     
c------------------------------------------------------------------------

        xl = px1(cube(1),ic)    !X coordinates of cube vertexes
        xh = px1(cube(1)+1,ic)

        yl = py1(cube(2),ic)    !Y coordinates of cube vertexes
        yh = py1(cube(2)+1,ic)

        zl = pz1(cube(3),ic)    !Z coordinates of cube vertexes
        zh = pz1(cube(3)+1,ic)

        xr = (x-xl) / (xh-xl)   !reduced variables
        yr = (y-yl) / (yh-yl)
        zr = (z-zl) / (zh-zl)

        Bp(1) = b1(cube(1)  ,cube(2)  ,cube(3)  ,ic) !ic-th component of B
        Bp(2) = b1(cube(1)+1,cube(2)  ,cube(3)  ,ic) ! on the eight cube
        Bp(3) = b1(cube(1)  ,cube(2)+1,cube(3)  ,ic) ! vertexes
        Bp(4) = b1(cube(1)+1,cube(2)+1,cube(3)  ,ic)
        Bp(5) = b1(cube(1)  ,cube(2)  ,cube(3)+1,ic)
        Bp(6) = b1(cube(1)+1,cube(2)  ,cube(3)+1,ic)
        Bp(7) = b1(cube(1)  ,cube(2)+1,cube(3)+1,ic)
        Bp(8) = b1(cube(1)+1,cube(2)+1,cube(3)+1,ic)

c------------------------------------------------------------------------
c     
c     computes interpolated ic-th component of B in (x,y,z)
c     
c------------------------------------------------------------------------

        res1(ic) = 
     +       Bp(1)*(1-xr)*(1-yr)*(1-zr) + 
     +       Bp(2)*xr*(1-yr)*(1-zr) +
     +       Bp(3)*(1-xr)*yr*(1-zr) + 
     +       Bp(4)*xr*yr*(1-zr) +
     +       Bp(5)*(1-xr)*(1-yr)*zr + 
     +       Bp(6)*xr*(1-yr)*zr +
     +       Bp(7)*(1-xr)*yr*zr +
     +       Bp(8)*xr*yr*zr


      enddo
      
c------------------------------------------------------------------------
c     
c     *** SECOND MAP ***
c     
c------------------------------------------------------------------------

c     second map is rotated by 180 degree along the Z axis. so change sign
c     of x and y coordinates and then change sign to Bx and By components
c     to obtain the correct result

      x=-x
      y=-y

      do ic=1,3

        cube(1)=INT((nx-1)*(x-px2min(ic))/(px2max(ic)-px2min(ic))) +1
        cube(2)=INT((ny-1)*(y-py2min(ic))/(py2max(ic)-py2min(ic))) +1
        cube(3)=INT((nz-1)*(z-pz2min(ic))/(pz2max(ic)-pz2min(ic))) +1
        
        if (cube(1).le.0) cube(1) = 1
        if (cube(2).le.0) cube(2) = 1
        if (cube(3).le.0) cube(3) = 1
        if (cube(1).ge.nx) cube(1) = nx-1
        if (cube(2).ge.ny) cube(2) = ny-1
        if (cube(3).ge.nz) cube(3) = nz-1

        xl = px2(cube(1),ic)
        xh = px2(cube(1)+1,ic)

        yl = py2(cube(2),ic)
        yh = py2(cube(2)+1,ic)

        zl = pz2(cube(3),ic)
        zh = pz2(cube(3)+1,ic)

        xr = (x-xl) / (xh-xl)
        yr = (y-yl) / (yh-yl)
        zr = (z-zl) / (zh-zl)

        Bp(1) = b2(cube(1)  ,cube(2)  ,cube(3)  ,ic)
        Bp(2) = b2(cube(1)+1,cube(2)  ,cube(3)  ,ic)
        Bp(3) = b2(cube(1)  ,cube(2)+1,cube(3)  ,ic)
        Bp(4) = b2(cube(1)+1,cube(2)+1,cube(3)  ,ic)
        Bp(5) = b2(cube(1)  ,cube(2)  ,cube(3)+1,ic)
        Bp(6) = b2(cube(1)+1,cube(2)  ,cube(3)+1,ic)
        Bp(7) = b2(cube(1)  ,cube(2)+1,cube(3)+1,ic)
        Bp(8) = b2(cube(1)+1,cube(2)+1,cube(3)+1,ic)

        res2(ic) = 
     +       Bp(1)*(1-xr)*(1-yr)*(1-zr) + 
     +       Bp(2)*xr*(1-yr)*(1-zr) +
     +       Bp(3)*(1-xr)*yr*(1-zr) + 
     +       Bp(4)*xr*yr*(1-zr) +
     +       Bp(5)*(1-xr)*(1-yr)*zr + 
     +       Bp(6)*xr*(1-yr)*zr +
     +       Bp(7)*(1-xr)*yr*zr +
     +       Bp(8)*xr*yr*zr

      enddo

c     change Bx and By components sign
      res2(1)=-res2(1)
      res2(2)=-res2(2)

c     change back the x and y coordinate signs
      x=-x
      y=-y


c------------------------------------------------------------------------
c     
c     average the two maps results
c     
c------------------------------------------------------------------------

      do ic=1,3
        res(ic)=(res1(ic)+res2(ic))/2
      enddo


      return
      end
*************************************************************************
*     
*     Subroutine inter_B_outer.f
*     
*     it computes the magnetic field in a chosen point x,y,z OUTSIDE the
*     magnetic cavity, using a trilinear interpolation of
*     B field measurements (read before by means of ./read_B.f)
*     the value is computed for the outer map
*     
*     needs:
*     - ../common/common_B_outer.f
*     
*     input: coordinates in m
*     output: magnetic field in T
*     
*************************************************************************

      subroutine inter_B_outer(x,y,z,res) !coordinates in m, magnetic field in T

      implicit double precision (a-h,o-z)
      include 'common_B.f'


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

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

      integer ic                
c     !index for B components:
c     ! ic=1 ---> Bx
c     ! ic=2 ---> By
c     ! ic=3 ---> Bz
      
      integer cube(3)           
c     !vector of indexes identifying the cube
c     ! containing the point of interpolation
c     ! (see later...)
      
      real*8 xl,xh,yl,yh,zl,zh  !cube vertexes coordinates

      real*8 xr,yr,zr           
c     !reduced variables (coordinates of the 
c     ! point of interpolation inside the cube)

      real*8 Bp(8)              
c     !vector of values of B component 
c     ! being computed, on the eight cube vertexes


c     LOWER MAP
c     ---> up/down simmetry
      zin=z
      if(zin.le.edgelzmax)z=-z

c------------------------------------------------------------------------
c     
c     *** MAP ***
c     
c------------------------------------------------------------------------

      do ic=1,3                 !loops on the three B components

c------------------------------------------------------------------------
c     
c     chooses the coordinates interval containing the input point
c     
c------------------------------------------------------------------------
c     e.g.:
c     
c     x1    x2    x3    x4    x5...               xN
c     |-----|-+---|-----|-----|---- ... ----|-----|
c     ~~~~~~~~x
c     
c     in this case the right interval is identified by indexes 2-3, so the 
c     value assigned to cube variable is 2

        cube(1)=INT((nox-1)*(x-poxmin(ic))/(poxmax(ic)-poxmin(ic))) +1
        cube(2)=INT((noy-1)*(y-poymin(ic))/(poymax(ic)-poymin(ic))) +1
        cube(3)=INT((noz-1)*(z-pozmin(ic))/(pozmax(ic)-pozmin(ic))) +1
        
c------------------------------------------------------------------------
c     
c     if the point falls beyond the extremes of the grid...
c     
c------------------------------------------------------------------------
c     
c     ~~~~~~~~~~x1    x2    x3...
c     - - + - - |-----|-----|----
c     ~~~~x
c     
c     in the case cube = 1
c     
c     
c     ...nx-2  nx-1  nx
c     ----|-----|-----| - - - + - -
c     ~~~~~~~~~~~~~~~~~~~~~~~~x
c     
c     in this case cube = nx-1

        if (cube(1).le.0) cube(1) = 1
        if (cube(2).le.0) cube(2) = 1
        if (cube(3).le.0) cube(3) = 1
        if (cube(1).ge.nox) cube(1) = nox-1
        if (cube(2).ge.noy) cube(2) = noy-1
        if (cube(3).ge.noz) cube(3) = noz-1


c------------------------------------------------------------------------
c     
c     temporary variables definition for field computation
c     
c------------------------------------------------------------------------

        xl = pox(cube(1),ic)    !X coordinates of cube vertexes
        xh = pox(cube(1)+1,ic)

        yl = poy(cube(2),ic)    !Y coordinates of cube vertexes
        yh = poy(cube(2)+1,ic)

        zl = poz(cube(3),ic)    !Z coordinates of cube vertexes
        zh = poz(cube(3)+1,ic)

        xr = (x-xl) / (xh-xl)   !reduced variables
        yr = (y-yl) / (yh-yl)
        zr = (z-zl) / (zh-zl)

        Bp(1) = bo(cube(1)  ,cube(2)  ,cube(3)  ,ic) !ic-th component of B
        Bp(2) = bo(cube(1)+1,cube(2)  ,cube(3)  ,ic) ! on the eight cube
        Bp(3) = bo(cube(1)  ,cube(2)+1,cube(3)  ,ic) ! vertexes
        Bp(4) = bo(cube(1)+1,cube(2)+1,cube(3)  ,ic)
        Bp(5) = bo(cube(1)  ,cube(2)  ,cube(3)+1,ic)
        Bp(6) = bo(cube(1)+1,cube(2)  ,cube(3)+1,ic)
        Bp(7) = bo(cube(1)  ,cube(2)+1,cube(3)+1,ic)
        Bp(8) = bo(cube(1)+1,cube(2)+1,cube(3)+1,ic)

c------------------------------------------------------------------------
c     
c     computes interpolated ic-th component of B in (x,y,z)
c     
c------------------------------------------------------------------------

        res(ic) = 
     +       Bp(1)*(1-xr)*(1-yr)*(1-zr) + 
     +       Bp(2)*xr*(1-yr)*(1-zr) +
     +       Bp(3)*(1-xr)*yr*(1-zr) + 
     +       Bp(4)*xr*yr*(1-zr) +
     +       Bp(5)*(1-xr)*(1-yr)*zr + 
     +       Bp(6)*xr*(1-yr)*zr +
     +       Bp(7)*(1-xr)*yr*zr +
     +       Bp(8)*xr*yr*zr


      enddo     

c     LOWER MAP
c     ---> up/down simmetry 
      if(zin.le.edgelzmax)then
         z=-z                   !back to initial ccoordinate
         res(3)=-res(3)         !invert BZ component
      endif

      return
      end