recursive SUBROUTINE em_field_custom(ele, param, s_rel, orb, local_ref, field, calcd, err_flag, &
                                         calc_potential, use_overlap, grid_allow_s_out_of_bounds, rf_time, used_eles)
  USE bmad !, except_dummy => em_field_custom
  use parameters_bmad
  use quad_scrape_parameters
  use em_field_interface
  use muon_mod

  IMPLICIT NONE

  type(ele_struct), target :: ele
  type(lat_param_struct) param
  type(coord_struct), intent(in) :: orb
  type(coord_struct) orb0
  type(em_field_struct) field, field_error
  type (ele_struct), pointer :: lord
  type (ele_pointer_struct), allocatable, optional :: used_eles(:)
  real(rp), intent(in) :: s_rel
  real(rp) t_rel
  real(rp) scale
 real(rp), optional :: rf_time
  logical, optional :: calcd, err_flag, calc_potential, grid_allow_s_out_of_bounds, use_overlap
  logical  out_of_range,local_ref
  logical grid, mpoles, err
  logical :: local_ref_frame = .false.
  logical  file_opened/.false./
  logical itsopen
  character(32) :: r_name = 'em_field_custom'
  character*3 fitn
  integer i
  integer, save :: n
  integer quad_number, param_number/0/
  integer, save :: lun

t_rel = orb%t 

mpoles = .false.
grid = .false.
if(ele%key == sbend$)then
 ! first check to see if the orbit is outside the boundaries of the map
  if(abs(orb%vec(1)) < 0.069 .or. index(ele%name,'QUAD') == 0)then  !inside the map
 ! then check to see if this element has a grid map
  if ((associated(ele%grid_field)))then
     grid = .true.
    elseif(ele%field_calc == refer_to_lords$)then
         do i=1,ele%n_lord
          lord =>pointer_to_lord(ele,i)
          if(associated(lord%grid_field))then
           grid = .true.
           exit
          endif
         end do
   endif
! next check for multipoles
   if(associated(ele%b_pole_elec) .or. associated(ele%a_pole_elec))mpoles = .true.
   if(grid .or. mpoles)then
    if(grid)ele%field_calc = fieldmap$
    if(mpoles)ele%field_calc = bmad_standard$ !multipole$

!  rf quad
!    if(index(ele%name,'QUAD') /= 0 .and. t_rel < max(maxval(t0_h_max(:)), maxval(t0_v_max(:))) .and. grid) then
    if(index(ele%name,'QUAD') /= 0 .and. grid) then       
       read(ele%name(5:5),'(i1)')quad_number
       if(ele%name(7:7) =='S')quad_number = -quad_number 

       call time_dep_quad_field(orb, ele, quad_number,  err_flag)
       orb0=orb
       orb0%vec(1:4)=0
       if(em_field_calc_verbose)then
          if(.not. file_opened .or. rfquad_params_reset)then
             lun=92
             inquire(unit=lun, opened=itsopen)
             if(itsopen)close(lun)
             param_number=param_number+1
             if(param_number < 10) write(fitn,'(i1)')param_number
             if(param_number >= 10) write(fitn,'(i2)')param_number           
             open(unit=lun, file=trim(directory)//'/quad_fields_vs_time_'//trim(fitn)//'.dat')
             write(lun,'(a)')' Fields on axis'
             write(lun, '(11a12)')'quad','time','x','xp','y','yp','z','delta','Ex','Ey','Ez'
             file_opened = .true.
             rfquad_params_reset = .false.             
          endif
          call em_field_calc(ele, param, s_rel, orb0, local_ref_frame, field)
          write(lun,'(i12,10es12.4)')quad_number, orb%t,orb%vec,field%E
       endif
       

!      scale = rf_quad(2)%amp_h*sin(rf_quad(2)%freq_h * t_rel + rf_quad(2)%phi_h) + rf_quad(2)%amp_v*sin(rf_quad(2)%freq_v * t_rel + rf_quad(2)%phi_v)
!      if(index(ele%name,'SHORT')/=0) then
!                 quad_params%short_quad_field_index(2)  = quad_params_0%short_quad_field_index(2)* (1.+ scale)
!                 call set_ele_real_attribute (ele, 'FIELD_INDEX',  quad_params%short_quad_field_index(2), err)
!                 call set_a_quad(ele,.false.)
  
!      endif
    endif
!

    call em_field_calc(ele, param, s_rel, orb, local_ref_frame, field, calcd, err_flag, &
         calc_potential, use_overlap, grid_allow_s_out_of_bounds, rf_time, used_eles)

    ele%field_calc = custom$
   else
    call em_field_custom_ring(ele, param, s_rel, t_rel, orb, local_ref, field, calcd, err_flag, &
                                    calc_potential, use_overlap, grid_allow_s_out_of_bounds, rf_time, used_eles)
   endif
  else ! outside map boundaries
    field%B = [0.0_rp, ele%value(B_field$), 0.0_rp]
    field%E=0
 endif
 
   field_error%B=0
   field_error%E=0
   if(param%ixx == 1)call field_errors(ele, param,s_rel, orb, field_error)
   field%B = field%B + field_error%B
   field%E = field%E + field_error%E


!  write(98,'(a,2es12.4,5es12.4,3es12.4)')ele%name, ele%s_start, s_rel,ele%s_start+s_rel,ele%value(rho$)+orb%vec(1),field%B,field_error%B

 else
  call em_field_custom_inj(ele, param, s_rel,t_rel, orb, local_ref,field,calcd, err_flag, &
                                         calc_potential, use_overlap, grid_allow_s_out_of_bounds, rf_time, used_eles) 
endif

!if(index(ele%name,'QUAD2')/= 0 .and. index(ele%name,'SHORT')/=0)then
! if(first)then
!    write(197,'(9a12)')'s_rel','t_rel','Bx','By,','Bz','Ex','Ey','Ez','ele%value(field_index$)'
!    n=0
!    first=.false.
! endif
!     n=n+1
!     if((n/100)*100 == n) write(197,'(9es12.4)')s_rel,t_rel,field%B,field%E, ele%value(field_index$)
!endif

return
end

recursive SUBROUTINE em_field_custom_ring(ele, param, s_rel, t_rel, orb, local_ref, field, calcd, err_flag, &
                                         calc_potential, use_overlap, grid_allow_s_out_of_bounds, rf_time, used_eles)
  USE bmad_struct
  USE bmad_interface
  use em_field_interface, dummy=> em_field_custom_ring

  IMPLICIT NONE


  type(ele_struct) :: ele
  type(lat_param_struct) param
  type(coord_struct) orb
  type(em_field_struct) field
  type (ele_pointer_struct), allocatable, optional :: used_eles(:)
  real(rp), intent(in) :: s_rel, t_rel
  real(rp) x,y,s,t
  real(rp), optional :: rf_time
  logical  out_of_range,local_ref
  logical, optional :: calcd, err_flag, calc_potential, grid_allow_s_out_of_bounds, use_overlap

  character(32) :: r_name = 'em_field_custom'

  x = orb%vec(1)
  y = orb%vec(3)
  s = s_rel
  t = orb%t  !t_rel            
  call fields(ele,x,y,s,t,field)

  RETURN
END SUBROUTINE em_field_custom_ring