subroutine calc_jwvar(lat, co, n_turns, coord, jvar)

  use bmad

  implicit none

  type (lat_struct), target :: lat
  type (coord_struct), target :: co
  integer, target ::  n_turns
  real(rp), target :: coord(6)
  real(rp) jvar(2)

  type (coord_struct) orbit

  integer i, j
  real(rp) betax, alphax, gammax
  real(rp) betay, alphay, gammay
  real(rp) jx, jy
  real(rp) jx_avg, jx_min, jx_max, jy_avg, jy_min, jy_max

  !
  
  betax=lat%ele(0)%a%beta
  gammax=lat%ele(0)%a%gamma
  alphax=lat%ele(0)%a%alpha
  betay=lat%ele(0)%b%beta
  gammay=lat%ele(0)%b%gamma
  alphay=lat%ele(0)%b%alpha
  
  orbit = co
  orbit%vec(:) = coord(:) + co%vec(:)
  jx_avg = 0.0
  jx_min = 10.0
  jx_max = 0.0
  jy_avg = 0.0
  jy_min = 10.0
  jy_max = 0.0

  do i=1, n_turns
     jx = (gammax*(orbit%vec(1)-co%vec(1))**2 + 2*alphax*(orbit%vec(1)-co%vec(1))*(orbit%vec(2)-co%vec(2)) + betax*(orbit%vec(2)-co%vec(2))**2)/2
     jy = (gammay*(orbit%vec(3)-co%vec(3))**2 + 2*alphay*(orbit%vec(3)-co%vec(3))*(orbit%vec(4)-co%vec(4)) + betay*(orbit%vec(4)-co%vec(4))**2)/2

     if (jx .le. jx_min) then
        jx_min = jx
     end if
     if (jx .ge. jx_max) then
        jx_max = jx
     end if

     if (jy .le. jy_min) then
        jy_min = jy
     end if
     if (jy .ge. jy_max) then
        jy_max = jy
     end if

     jx_avg = jx_avg + jx
     jy_avg = jy_avg + jy

     do j=1, lat%n_ele_track       
        call track1(orbit, lat%ele(j),lat%ele(j)%branch%param, orbit )
     end do ! j element
     
     if (orbit%state .ne. alive$) then
        write (6, '(a)') 'calc_jwvar error: '
        write (6, '(a,i4,a,i4,a,i4,a)') 'Particle lost at turn: ', i, ' element ', orbit%ix_ele+1, ' plane: ', orbit%state, ' (3:-x; 4:+x; 5:-y; 6:+y; 7:z)'
        write (6, '(a)') 'Change initial coordinates in the constraint file (&TRACKING) and try again!'
        write (6, '(a,6f14.6)') 'Current initial coordinates: ', coord
        write (6, '(a)') 'To change initial coord, need read_tracking = .true. in &SINGLE_VARS'
        write (6, '(a)') 'Stop now!'
        stop
     end if
     
  end do !i turn
  jx_avg = jx_avg/n_turns
  jy_avg = jy_avg/n_turns

  jvar(1) = (jx_max-jx_min)/jx_avg
  jvar(2) = (jy_max-jy_min)/jy_avg

  
end subroutine calc_jwvar