SUBROUTINE field_errors(ele,param,s_rel,orb, field_error)
  USE bmad
  use parameters_bmad
  implicit none

interface
  subroutine read_multipoles(npoints,nterms,PP,AA,BB)

   use parameters_bmad
   use sim_utils
   implicit none
   integer npoints, nterms
   real(rp),allocatable :: PP(:), AA(:,:), BB(:,:)
 end subroutine read_multipoles
end interface

  interface
 subroutine generalize_cyl_exp(azimuthal_exp_z, azimuthal_exp_r, r,phi,z, field)
  use bmad
  implicit none
  type(em_field_struct) field
  real(rp) r,phi,z
  character *120 azimuthal_exp_z, azimuthal_exp_r
end subroutine
end interface

  type(ele_struct), target :: ele
  type(lat_param_struct) param
  type(coord_struct) orb
  type(em_field_struct)  field_error, field
  real(rp) s_rel
  real(rp) phi, phase, br

  real(rp), allocatable, save :: PP(:), AA(:,:), BB(:,:)
  integer i
  integer n, ix_phi
  integer,  save :: nterms, npoints
  real(rp) A(0:4)/68.9569,21.6858,30.4212,0.,0./, B(0:4)/0,-28.6175,26.6923,0.,0./
  real(rp) d/68.9569/
  real(rp), save :: r0
  real(rp) psi,r
  real(rp) z
  real(rp) Bfield(3)
  real(rp) b_ref, dphi, delta,theta
  logical first/.true./
  character *120 azimuthal_exp_z/'BzFourier20170628_LogID983.dat'/, azimuthal_exp_r/'BrFourier2016.dat'/

!  Br= A(0)r0/r + sum  [A(n)*cos(n phi)+B(n)*sin(n phi)](r0/r)**(n+1)
!  psi= A(0)r0ln(r) - 1/n sum  [A(n)*cos(n phi)+B(n)*sin(n phi)] r0**(n+1)/r**n

  psi=0
  field_error%B=0
  Bfield = 0
  phi = (ele%s_start+s_rel)/param%total_length * twopi
  !r0 = ele%value(rho$)

  z = orb%vec(3)

  if(multipole_field_file(1:1)  /= '')then
     if(first)then
        call read_multipoles(npoints,nterms,PP,AA,BB)
!        PP(:)=0
        AA(:,:)=0
        BB(:,:)=0
        r0 = 0.045
        first=.false.
     endif


   b_ref = ele%value(p0c$)/ele%value(rho$)/c_light
   phi = (ele%s_start + s_rel)/param%total_length * 360
   ix_phi = phi /360 * npoints
   dphi = phi - PP(ix_phi)
   if(ix_phi >= npoints)then
!      print *,' ix_phi =', ix_phi
   endif
   
   delta = PP(ix_phi+1)-PP(ix_phi)
   do n=0,nterms
      if(delta > 0)then
         A(n) = (AA(ix_phi,n) * (phi-PP(ix_phi)) +AA(ix_phi+1,n)*(PP(ix_phi+1)-phi))/delta
         B(n) = (BB(ix_phi,n) * (phi-PP(ix_phi)) +BB(ix_phi+1,n)*(PP(ix_phi+1)-phi))/delta
      else
         A(n) = AA(ix_phi,n) 
         B(n) = BB(ix_phi,n) 
      endif
      
   end do
   theta = atan2(orb%vec(3),orb%vec(1))

      if(multipole_field_file(1:1) == '')then
         r = ele%value(rho$) + orb%vec(1)
         Bfield(1) = A(0) * ele%value(rho$)/r
         psi = A(0) * ele%value(rho$) *log(r)
         phi = (ele%s_start + s_rel)/param%total_length * twopi
         do n=0,nterms
           psi = psi - (1/n)*(A(n)*cos(n*phi)+B(n)*sin(n*phi))* ele%value(rho$)**(n+1)/r**n
           Bfield(1) = Bfield(1) +   (A(n)*cos(n*phi)+B(n)*sin(n*phi))*(ele%value(rho$)/r)**(n+1)  !1/r d psi/d r
           Bfield(3) = Bfield(3) +  (A(n)*sin(n*phi)-B(n)*cos(n*phi))*(ele%value(rho$)/r)**(n+1)   !1/r d psi/d phi
         end do    
      else
         do n=0,nterms
            Bfield(2) = Bfield(2) + (orb%vec(1)/r0)**n * (B(n)*cos(n*theta) + A(n) * sin(n*theta))
         end do
      endif
      
   
   field_error%B(1:3) = Bfield(1:3)*b_ref*1.e-6

 else !3d
 
 call generalize_cyl_exp(azimuthal_exp_z, azimuthal_exp_r, r,phi,z, field)
 field_error%B(1:3) = field%B(1:3)*b_ref*1.e-6
 field_error%E=0



 endif
!  write(97,'(a,2es12.4,5es12.4)')ele%name, ele%s_start, s_rel,phi,theta,field_error%B
 return
 end subroutine

 subroutine read_multipoles(npoints,nterms,PP,AA,BB)

   use parameters_bmad
   use sim_utils
   implicit none
   integer lun
   integer npoints, nterms
   integer reason
   real(rp),allocatable :: PP(:), AA(:,:), BB(:,:)
   character*200 string
   npoints = 0
   lun = lunget()
   open(unit = lun, file = trim(multipole_field_file))
   do while(.true.)
      read(lun,'(a)', IOSTAT = reason)string
      if(reason < 0)exit
      if(index(string,'#')/= 0 .or. string(1:2) == '  '.or. index(string,'angle')/=0)cycle
      npoints = npoints + 1
   end do
   allocate(PP(0:npoints+1), BB(0:npoints+1,0:4),AA(0:npoints+1,0:4))
   rewind(lun)
   npoints = 0
   nterms = 4
   do while(.true.)
      read(lun,'(a)', IOSTAT = reason)string
      if(reason < 0)exit
      if(index(string,'#')/= 0 .or. string(1:2) == '  '.or. index(string,'angle')/=0)cycle
      npoints = npoints + 1
      read(string,*)PP(npoints),BB(npoints,0),BB(npoints,1),AA(npoints,1),BB(npoints,2),AA(npoints,2),BB(npoints,3),AA(npoints,3),BB(npoints,4),AA(npoints,4)
!      print '(i10,10es12.4)', npoints, PP(npoints),BB(npoints,0),BB(npoints,1),AA(npoints,1),BB(npoints,2),AA(npoints,2),BB(npoints,3),AA(npoints,3),BB(npoints,4),AA(npoints,4)
      
   end do
   PP(npoints+1) = PP(npoints)
  print '(a27,a)',' Read multipoles from file: ', trim(multipole_field_file)
   return
 end subroutine read_multipoles