subroutine compute_derivatives(fringe_or_inflector,map)

   use magfield
   use parameters_bmad
   use magfield_interface, dummy => compute_derivatives
   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)
   integer fringe_or_inflector

! next compute dB(1:3)/dxyz(1:3)
if(field_file(fringe_or_inflector)%type == 4)then
  print '(a)', ' Field file name = ', field_file(fringe_or_inflector)%name,' type = ',field_file(fringe_or_inflector)%type
  print '(a)', 'No derivatives computed for uniform field'
 return
endif
!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 ! change in x,y, or z

     !at each grid point calculate B(1:3) displaced by +1,-1,0 grid points from central value

       !if ix ==1,  dxyz(+-1,0)%B(1:3) is the B(1:3) displaced in x by +1,-1,0 grid point
       !if ix == 2  dxyz(+-1,0)%B(1:3) is the B(1:3) displaced in y by +1,-1,0 grid point
       !if ix == 3  dxyz(+-1,0)%B(1:3) is the B(1:3) displaced in z by +1,-1,0 grid point
      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

      call quadratic_fit(dxyz,fringe_or_inflector, b1,b2)   !b1(1:3) is the first derivative in ix=1,2 or 3 direction

      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

   end do
  end do
end do

return
end subroutine