subroutine particle_tracking (lat, con, con_value)

  use bmad
  use constraints_mod

  implicit none

  type (bmad_common_struct) com_saved
  type (lat_struct) lat
  type (constraint_struct) con
  type (coord_struct) action0, action
  type (coord_struct), allocatable, save :: orbit(:)

  integer ip, ip0, ip1, ip2
  integer it, track_state

  real(rp) con_value, dja_max, djb_max, j0, rrr

  call reallocate_coord( orbit, lat%n_ele_max )

! set parameters

  com_saved = bmad_com
  bmad_com%aperture_limit_on = .true.

  lat%ele%value(x1_limit$) = 1
  lat%ele%value(x2_limit$) = 1
  lat%ele%value(y1_limit$) = 1
  lat%ele%value(y2_limit$) = 1

! set rf voltage for correct synch tune

  call set_z_tune (lat%branch(0), -con%q_z * twopi)

! misc init

  call convert_coords ('LAB', con%track%start, lat%ele(0), 'ACTION', action0)
  ip0 = nint(sqrt(float(con%track%n_particle)))
!

  dja_max = 0
  djb_max = 0

  do ip = 0, con%track%n_particle - 1

    action = action0
    ip1 = mod(ip, ip0)
    ip2 = (ip - ip1) / ip0
    action%vec(2) = ip1 * twopi
    action%vec(4) = ip2 * pi

    call convert_coords ('ACTION', action, lat%ele(0), 'LAB', orbit(0))

    do it = 1, con%track%n_turn
      call track_all (lat, orbit, 0, track_state)
      orbit(0) = orbit(lat%n_ele_track)
      call convert_coords ('LAB', orbit(0), lat%ele(0), 'ACTION', action)
      dja_max = max(dja_max, abs(action%vec(1) - action0%vec(1)))
      djb_max = max(djb_max, abs(action%vec(3) - action0%vec(3)))
      if (track_state /= moving_forward$) then
        call convert_coords ('LAB', orbit(track_state-1), lat%ele(0), 'ACTION', action)
        rrr = float(con%track%n_turn) / it
        dja_max = max(dja_max, abs(action%vec(1) - action0%vec(1))) * rrr
        djb_max = max(djb_max, abs(action%vec(3) - action0%vec(3))) * rrr
        exit
      endif
    enddo

  enddo

  bmad_com = com_saved
  j0 = max(action0%vec(1), action0%vec(3))
  con_value = max(dja_max/j0, djb_max/j0)


end subroutine particle_tracking