subroutine test_quad(lat,nbranch,ix_quad)
  use bmad
  use muon_mod
  use parameters_bmad
  use muon_interface
  implicit none
  
  type (lat_struct), target::lat
  type (branch_struct), pointer:: branch ! ring
  type (coord_struct) co
  type (ele_struct) ele
  type (em_field_struct) field

  integer i
  integer nbranch
  integer ix_quad
  integer lun
  logical err_flag, local_ref_frame/.false./

  real(rp) x, delta, s_pos

  lun=lunget()
  open(unit=lun, file=trim(directory)//'/quad_field.dat')

  branch =>lat%branch(nbranch)
  ele = branch%ele(ix_quad)
  print *,' element = ', ele%name, ele%s

  s_pos = ele%s/2
  if(associated(ele%b_pole_elec)) s_pos=0
  delta= 0.001
  call init_coord(co)
   write(lun,'(a12,a36)')'x',' field%E '
   x=-0.045
   do while(x <0.045)
    co%vec=0
     x=x+delta
     co%vec(1) = x

    call em_field_calc(ele, lat%param, s_pos, co,local_ref_frame,field,err_flag=err_flag)
    write(lun, '(4es12.4)')x,field%E
   end do

  if(associated(ele%b_pole_elec)) then
   do i=0, n_pole_maxx
!    print '(i10,1x,2es12.4)',i, ele%b_pole_elec(i),ele%a_pole_elec(i)
    ele%a_pole_elec(i) =0
    if(i /= 1)ele%b_pole_elec(i)=0
   end do

   write(lun,'(/,a,/)')'#    use only m=2'
   write(lun,'(a12,a36)')'x',' field%E '
   x=-0.045
   do while(x <0.045)
    co%vec=0
     x=x+delta
     co%vec(1) = x

    call em_field_calc(ele, lat%param, s_pos, co,local_ref_frame,field,err_flag=err_flag)
    write(lun, '(4es12.4)')x,field%E
   end do
 endif
 return
end