program tuneshift_from_impedance

  use bmad
  use sim_utils
  use tune_shift_mod

 implicit none

 type (lat_struct), target :: ring
 type (coord_struct), allocatable :: bmad_orbit(:), closed_orbit(:)
 type (coord_struct) start_orb
 type(wake_ele_struct) :: wakes(100)
 type(ele_tuneshift_struct) :: ele_tuneshifts(200)
 type (normal_modes_struct) modes


 ! File handling:
 character(200) :: init_file
 character(200) :: lat_file

 ! other
 real(rp) :: total_tuneshift_x = 0., total_tuneshift_y = 0., sigma_z=0., current = 1
 integer :: ix=0, jx=0
 
!=======================================

 ! Read in parameters from file:
 namelist /parameters/ lat_file, sigma_z, current

 ! file handling: user-specified filename
 if (command_argument_count() .ne. 1) then 
    init_file = 'tuneshift_from_impedance.in' ! use default input file name 
 else ! user specified an input filename
    call get_command_argument(1, init_file)
 endif

 ! open and read contents of file
 open(1, file=init_file)
 read(1, nml=parameters)
 close(1)

 ! formality:
 bmad_com%auto_bookkeeper = .false.
 
 ! read in the lattice file:
 call bmad_parser(lat_file, ring)
 
 ! set up lattice 
 call set_on_off(rfcavity$, ring, on$)
 call closed_orbit_calc(ring, closed_orbit, 6)
 call lat_make_mat6(ring, -1, closed_orbit)
 call twiss_at_start(ring)
 call twiss_propagate_all(ring)
 call radiation_integrals(ring, closed_orbit, modes)
 
 write(*,*) "Lattice ", trim(ring%lattice), " read."
 
 if(sigma_z > 0.) then 
 		modes%sig_z = sigma_z
 end if
 ! read in wake parameters from init_file
 call parse_wakes(init_file, wakes, modes%sig_z)


 ! rescale impedances from T3P according to bunch length
 call scale_wakes(wakes, modes%sig_z)

 ! calculate the tuneshift contribution from each element
 call calc_tuneshift(init_file, wakes, ele_tuneshifts, ring, modes%sig_z, current)

! ! sum up the results:
! do ix=1,size(ele_tuneshifts)
!    if (ele_tuneshifts(ix)%name == '') cycle
!    do jx=1,size(ele_tuneshifts(ix)%tuneshift)
!       total_tuneshift_x = total_tuneshift_x + ele_tuneshifts(ix)%tuneshift(jx)%x
!       total_tuneshift_y = total_tuneshift_y + ele_tuneshifts(ix)%tuneshift(jx)%y
!    enddo
! enddo
!
!! print tuneshifts in kHz:
! write(*,*) 'I_bunch = ', current, 'mA'
! write(*,*)'tuneshift x = ', total_tuneshift_x, 'kHz'
! write(*,*)'tuneshift y = ', total_tuneshift_y, 'kHz'

end program tuneshift_from_impedance