subroutine quadratic_fit(dxyz,fringe_or_inflector,b1,b2)
       use magfield
       use parameters_bmad
       use nr
       implicit none
       real(rp) dx
       type(magfield_struct) dxyz(-1:1)
       real(rp) b0(3), b1(3), b2(3)
!       real(rp)  A(3,3), V(3,3)
       integer ix, i 
       integer fringe_or_inflector
       dx=field_file(fringe_or_inflector)%grid_spacing

!print '(3(4x,3e12.4))', dxyz(-1)%B,dxyz(0)%B,dxyz(1)%B
!pause
!         A(1,1:3)= 1.
!         A(2,1:3) = (/-0.5,0.,0.5/)
!         A(3,1:3) = (/0.25, 0., 0.25/)
!   print '(a,3es12.4)',' dxyz(-1:1)%B(1) - Bx ', dxyz(-1:1)%B(1)
!   print '(a,3es12.4)',' dxyz(-1:1)%B(2) - By ', dxyz(-1:1)%B(2)
!   print '(a,3es12.4)',' dxyz(-1:1)%B(3) - Bz ', dxyz(-1:1)%B(3)
!do i=1,3
!     print '(a,i3,3es12.4)',' A row = ', i, A(1:3,i)  
!end do
!     print '(/,a,3es12.4)',' V before ' , V(1,1:3)

!         do i=1,3
!          V(1,1:3) = dxyz(-1:1)%B(i)
!          call gaussj(A,V)
!          B1(i) = V(1,2)
!          B2(i) = V(1,3)
!         end do 
!     print '(/,a,3es12.4)',' V after ' , V(1,1:3)

         b0(1:3) = dxyz(0)%B(1:3)
         b1(1:3) = (dxyz(1)%B(1:3) - dxyz(-1)%B(1:3))/2/dx ! derivative of B_ix wr x,y,z
         b1(1:3) = (dxyz(1)%B(1:3) - dxyz(0)%B(1:3))/dx ! derivative of B_ix wr x,y,z
         b2(1:3) = (dxyz(1)%B(1:3) + dxyz(-1)%B(1:3) - 2*b0(1:3))/2/dx**2 ! derivative of B_ix wr x,y,z


!          print '(a,i3,6e12.4)',' ix, b0,b1,b2 ',ix,b0(1:3), b1(1:3)
      return
      end subroutine quadratic_fit