--- gpamela/gpfield/inter_B.F	2005/12/06 01:07:23	3.1
+++ gpamela/gpfield/inter_B.F	2005/12/22 16:25:03	3.2
@@ -31,7 +31,8 @@
       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
@@ -42,22 +43,90 @@
         res(ip)=0.
       enddo
 
-
 c------------------------------------------------------------------------
 c     
-c     check if the point falls inside the interpolation volume
+c     check if the point falls inside the interpolation volumes
 c     
 c------------------------------------------------------------------------
-      
-      if((x.ge.edgexmin).and.(x.le.edgexmax)
+
+*     -----------------------
+*     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
+     $     .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
 
-c     TODO: add inter_B_outer subroutine(s) for the external magnetic field map
 
       return
       end