subroutine interpolate_field(x,B,b1,b2)
   use bmad
   implicit none
   real(rp) x(3), B(3),z, b1(3,3),b2(3,3)
   real(rp) xgrid(3), dx(3)
   integer ind(3),i

   ind(1:3) = (x(1:3))*2+1
   xgrid(1:3) = (ind(1:3)-1)*0.5
   dx(1:3) = x(1:3) - xgrid(1:3)
!  print '(a,6es12.4)',' dx, B ', dx,B
   forall(i=1:3) B(i) = B(i) + dot_product(b1(i,:),dx(:))    
!  print '(a,3es12.4)', 'B(1:3) ', B(1:3)
   return
  end subroutine