*************************************************************************
*     
*     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 "gpfield.inc"


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