subroutine compute_derivatives(map)

   use magfield
   use magfield_interface

   implicit none
   type(magfield_struct), allocatable :: map(:,:,:)
   type(magfield_struct) x0, xmax, xmin, dxyz(-1:1)
   integer i,j,k,l, ntot(3), ix,jx
   real(rp)  b1(3), b2(3)


! next compute dB(1:3)/dxyz(1:3)

!print *
ntot(1) = size(map,1)
ntot(2) = size(map,2)
ntot(3) = size(map,3)

print '(a,3i4)',' ntot ', ntot(1:3)
do j=1,ntot(1)
  do k=1,ntot(2)
    do l=1,ntot(3)
     do ix=1,3
      do jx=-1,1,1

        dxyz(jx)%B(1:3) = map(j,k,l)%B(1:3)
        if(ix == 1 .and. j+jx >= 1 .and. j+jx <= ntot(ix))then
          dxyz(jx)%B(1:3) = map(j+jx,k,l)%B(1:3)
         elseif(ix == 2 .and. k+jx >= 1 .and. k+jx <= ntot(ix))then
          dxyz(jx)%B(1:3) = map(j,k+jx,l)%B(1:3)
         elseif(ix == 3 .and. l+jx >= 1 .and. l+jx <= ntot(ix))then
          dxyz(jx)%B(1:3) = map(j,k,l+jx)%B(1:3)
         else
          cycle
        endif
      end do


!      print '(a,i,3es12.4)','2  ix, dxyz', ix,  dxyz(-1:1)%B(2)
      call quadratic_fit(dxyz, b1,b2)
!      print '(a,3es12.4)','b1(1:3) ', b1
      map(j,k,l)%dB(1:3,ix) = b1(1:3) !first derivative at dB_1:3/dx_1:3
      map(j,k,l)%d2B(1:3,ix) = b2(1:3) !second derivative at d2B_1:3/d2x_1:3 /2
     end do
!      write(11,'(6F12.4)') map(j,k,l)%x(1:3),map(j,k,l)%B(1:3)
   end do
  end do
end do

return
end subroutine