!+
!
! Subroutine : CHROM_CALC_z (LAT, DELTA_E, CHROM_X, CHROM_Y, c_mat_chrom)
!
! Description: Subroutine to calculate the chromaticities by computing the 
!              tune change when then energy is changed.
!
! Arguments  : Input:
!                    LAT      -- lat_struct: Lat
!                    DELTA_E   -- Real(rp): Delta energy used for the 
!                                 calculation. If 0 then default of 1.0e-4 
!                                 is used.
!
!              Output:
!                    DELTA_E -- Real(rp): Set to 1.0e-4 if on input 
!                                DELTA_E =< 0.
!                    CHROM_X -- Real(rp): Horizontal chromaticity.
!                    CHROM_Y -- Real(rp): Vertical chromaticity.         
!-

subroutine chrom_calc_z (lat, delta_e, chrom_x, chrom_y, c_mat_chrom, err_flag)

  use bmad

  implicit none

  type (lat_struct)  lat
  type (lat_struct), save :: lat2
  type (coord_struct)  c0 
  type (coord_struct), allocatable, save :: coord(:)

  integer i, key, status

  real(rp) high_tune_x, high_tune_y, low_tune_x, low_tune_y
  real(rp) delta_e, chrom_x, chrom_y, c_mat_low(2,2), c_mat_high(2,2)
  real(rp) c_mat_chrom(2,2)

  logical err_flag, err

  !

  err_flag = .true.

  call reallocate_coord( coord, lat%n_ele_max )

  if (delta_e <= 0) delta_e = 1.0e-4
  lat2 = lat

! lower energy tune

  coord(0)%vec(:) = 0.
  coord(0)%vec(6) = -delta_e
  call closed_orbit_calc( lat2, coord, 4, err_flag = err); if (err) return
  call track_all (lat2, coord, err_flag = err); if (err) return
  call lat_make_mat6 (lat2, -1, coord, err_flag = err); if (err) return

  call twiss_at_start (lat2, status); if (status /= ok$) return
  low_tune_x = lat2%a%tune / twopi
  if (low_tune_x < 0) low_tune_x = 1 + low_tune_x
  low_tune_y = lat2%b%tune / twopi
  if (low_tune_y < 0) low_tune_y = 1 + low_tune_y
  c_mat_low = lat2%ele(0)%c_mat


! higher energy tune

  coord(0)%vec(:) = 0.
  coord(0)%vec(6) = delta_e
  call closed_orbit_calc( lat2, coord, 4, err_flag = err); if (err) return
  call track_all (lat2, coord, err_flag = err); if (err) return
  call lat_make_mat6 (lat2, -1, coord, err_flag = err); if (err) return

  call twiss_at_start (lat2, status); if (status /= ok$) return
  high_tune_x = lat2%a%tune / twopi
  if (high_tune_x < 0) high_tune_x = 1 + high_tune_x
  high_tune_y = lat2%b%tune / twopi
  if (high_tune_y < 0) high_tune_y = 1 + high_tune_y
  c_mat_high = lat2%ele(0)%c_mat

! compute the chrom

  chrom_x = (high_tune_x - low_tune_x) / (2 * delta_e)
  chrom_y = (high_tune_y - low_tune_y) / (2 * delta_e)
  c_mat_chrom = (c_mat_high - c_mat_low)/ (2 * delta_e)

  err_flag = .false.

end subroutine chrom_calc_z