!+ 
! Program ring_aperture
!
! Program to determine the dynamic aperture of a ring
!
! Modules Needed:
!  ???
!
! Input (namelist)
!   ring_aperture.in
!
! Output
!   [particle file name].out
!-

program ring_aperture

use bmad_struct
use dynamic_aperture_mod
use ring_aperture_mod
use quick_plot
!use random_mod
!$ use omp_lib
use ptc_layout_mod
use beam_mod
 use touschek_mod

implicit none

! OpenMP shared memory
type (lat_struct), allocatable :: omp_lat(:)
type (aperture_struct), allocatable :: omp_aperture(:)
integer, allocatable :: omp_counter(:)

type (lat_struct), target :: lat
type (branch_struct), pointer :: branch
type (ele_struct), pointer :: ele0
type (bunch_struct)  :: bunch, bunch_end
type (coord_struct) :: orb0, orb_mode
type (coord_struct), allocatable :: closed_orb(:), ptc_closed_orb(:)
type (track_input_struct) :: track_input
type (aperture_struct) :: aperture
type (aperture_struct), allocatable :: aperture_list(:)
!type (ring_aperture_param_struct) :: ra_param
type (normal_modes_struct) :: normal_mode, normal_mode_ptc
type (beam_init_struct) :: beam_init
type (taylor_struct) :: t_map_1_turn(6) 
type (bunch_params_struct) :: bunch_params
type (amplitude_block_struct) :: amplitude_block

real(rp) :: Qa0, Qb0, Qa1, Qb1, Qa, Qb, sigma_mat(6,6), sigma_x, sigma_y, sigma_z
real(rp) :: emit_a, emit_b, sigma_px, sigma_py, sigma_pz, rel_gamma
real(rp) :: nJa, nJb, nJz
real(rp) :: step_a, step_b, step_z
real(rp), allocatable :: x_arr(:), y_arr(:), angle_list(:)
real(rp) :: theta_xy, theta_min, theta_max, d_theta, e_init, d_pz, pz_init
real(rp) :: ttouschek

character(100) :: in_file, lat_name, lat2_name, lat_path, base_name, particle_file_name, outfile_name

integer :: outfile, namelist_file
integer :: plot_every_n, counter
integer :: i_dim, n_angles, n_sigma_pz, i_pz, n_turns, n_particle
integer :: ix_ele, n_monitor_eles
integer :: n_char
integer :: plot_id
integer :: omp_n, omp_i, i, j, k, n_sets

logical :: save_tracks, verbose, err, plot_on, omp_on
logical :: parallel, bunch_tracking_on, ptc_calc_on, angle_scan_on, taylor_tracking_on
logical :: radiation_damping_on, radiation_fluctuations_on, rf_on, amplitude_scan_on
logical :: frequency_drift_on, touschek_calc_on
character(1), parameter :: carriage_return = achar(13) 

character(30), parameter :: r_name = 'ring_aperture'

namelist / ring_aperture_params / &
	 nJa, nJb, nJz, &
    lat_name, lat2_name, particle_file_name, &
    verbose, plot_on, parallel, plot_every_n, n_sigma_pz, n_turns, &
    n_angles, theta_min, theta_max, bunch_tracking_on, beam_init, ptc_calc_on, &
    angle_scan_on, taylor_tracking_on, radiation_damping_on, radiation_fluctuations_on, rf_on, &
    amplitude_scan_on, frequency_drift_on, touschek_calc_on, step_a, step_b, step_z
   


!------
print '(a)',  ''
print '(a)',  ' __          __           __   ___  __  ___       __   ___ '
print '(a)',  '|__) | |\ | / _`     /\  |__) |__  |__)  |  |  | |__) |__  '
print '(a)',  '|  \ | | \| \__>    /~~\ |    |___ |  \  |  \__/ |  \ |___ '                                                        
print '(a)',  ''


!------------------------------------------
!Defaults for namelist
lat_name = 'lat.bmad'
lat2_name = ''
particle_file_name = 'ReferenceParticles.dat'
verbose = .false.
parallel = .false.
plot_on = .false.
n_sigma_pz = 0
n_turns = 1000
n_angles=181
theta_min = 0
theta_max = twopi
bunch_tracking_on  = .false. 
beam_init%n_particle = 1
beam_init%n_bunch = 1
plot_every_n = 1
ptc_calc_on = .false.
touschek_calc_on = .false.
angle_scan_on = .true.
taylor_tracking_on = .false.
radiation_damping_on = .true. 
radiation_fluctuations_on = .true. 
rf_on = .true. 
amplitude_scan_on = .false.
frequency_drift_on = .false. 
nJa = 10.0_rp
nJb = 10.0_rp
nJz = 10.0_rp
step_a = 1.0_rp
step_b = 1.0_rp
step_z = 1.0_rp

! OpenMP control
omp_on = .false.
!$ print *, 'OpenMP is on'
!$ omp_on = .true.

!Read namelist
in_file = 'ring_aperture.in'
if (command_argument_count() > 0) call get_command_argument(1, in_file)

namelist_file = lunget()
print *, 'Opening: ', trim(in_file)
open (namelist_file, file = in_file, status = "old")
read (namelist_file, nml = ring_aperture_params)
close (namelist_file)

!TEMP
!ra_param%verbose = verbose
!ra_param%plot_on = plot_on
!ra_param%id = 1

!Trim filename
n_char= SplitFileName(lat_name, lat_path, base_name) 

!Prepare output file name
call file_suffixer (particle_file_name, outfile_name, '.out', .true.)

!Parse Lattice
!write (*, *) 'Parsing '//trim(lat_name)
call bmad_parser (lat_name, lat)
!branch => lat%branch(0)
!print *, 'parsed: ',  lat%input_file_name

!Parse additional settings
if (lat2_name /= '') then
  print *, 'Parsing: '//trim(lat2_name)
  call bmad_parser2 (lat2_name, lat)
endif

!if (branch%param%particle /= electron$ ) then
!  call out_io (s_warn$, r_name, 'LATTICE NOT SET UP FOR ELECTRONS')
!endif


!------------------------------------------
!Track through multiple elements

!Import  BMAD-T style particles
!call import_time_distribution(particle_file_name, &
!        start_particles)


!Switch all elements to time_runge_kutta$ tracking
!do ele_id = 1, branch%n_ele_track
!   branch%ele(ele_id)%tracking_method = time_runge_kutta$
!end do

!Open output file
outfile = lunget()
open(outfile, file = outfile_name)

!Track electrons
!branch%param%particle = electron$

write(*, *) 'Lattice:', trim(lat%input_file_name)

!--------------------
! RF on/off
if (rf_on) then
  call set_on_off (rfcavity$, lat, on$)
  print *, 'RF is ON' 
  i_dim=6
else
  call set_on_off (rfcavity$, lat, off$)
  print *, 'RF is OFF'
  i_dim=4
endif

!--------------------
! Radiation 
bmad_com%radiation_damping_on = radiation_damping_on
bmad_com%radiation_fluctuations_on = radiation_fluctuations_on
print *, 'bmad_com%radiation_damping_on: ', bmad_com%radiation_damping_on
print *, 'bmad_com%radiation_fluctuations_on: ', bmad_com%radiation_fluctuations_on

!--------------------
! Apertures (CESR is roughly 9cm wide by 5 cm tall)
do i = 1, lat%n_ele_max
  if(lat%ele(i)%value(x1_limit$) == 0) lat%ele(i)%value(x1_limit$) = 0.045
  if(lat%ele(i)%value(x2_limit$) == 0) lat%ele(i)%value(x2_limit$) = 0.045
  if(lat%ele(i)%value(y1_limit$) == 0) lat%ele(i)%value(y1_limit$) = 0.025
  if(lat%ele(i)%value(y2_limit$) == 0) lat%ele(i)%value(y2_limit$) = 0.025
enddo

!--------------------
! Closed orbit calc

call reallocate_coord(closed_orb, lat)
!allocate(orb0(0:lat%n_ele_track))

call closed_orbit_calc(lat, closed_orb, i_dim)
write (*,'(a, 6es15.7)') 'Closed orbit:', closed_orb(0)%vec

! orb and element at 0 for convenience
orb0=closed_orb(0)
ele0 => lat%ele(0)



!---------------------
! Emittance calc


! Prepare for radiation integrals
call twiss_at_start(lat)
call twiss_propagate_all(lat)
!print *, 'radiation_integrals'
call radiation_integrals (lat, closed_orb, normal_mode)

emit_a = normal_mode%a%emittance
emit_b = normal_mode%b%emittance
sigma_x = sqrt(ele0%a%beta*normal_mode%a%emittance +  (ele0%x%eta*normal_mode%sigE_E)**2)
sigma_y = sqrt(ele0%b%beta*normal_mode%b%emittance +  (ele0%y%eta*normal_mode%sigE_E)**2)
sigma_z = normal_mode%sig_z
sigma_px = sqrt(ele0%a%gamma*normal_mode%a%emittance +  (ele0%x%etap*normal_mode%sigE_E)**2)
sigma_py = sqrt(ele0%b%gamma*normal_mode%b%emittance +  (ele0%y%etap*normal_mode%sigE_E)**2)
sigma_pz = normal_mode%sigE_E
rel_gamma = ele0%value(e_tot$)/mass_of(lat%param%particle)

write (*, '(a)') '__________________________'
write (*, '(a, f15.4)') 'Relativistic gamma: ', rel_gamma
write (*, '(a, f15.7, a)') 'a emittance', emit_a*1e9, ' nm'
write (*, '(a, f15.7, a)') 'b emittance', emit_b*1e9, ' nm'
write (*, '(a, f10.7, a)') 'sigmaE/E', sigma_pz, ' '
write (*, '(a, es15.7, a)') 'beta_a', ele0%a%beta
write (*, '(a, es15.7, a)') 'beta_b', ele0%b%beta
write (*, '(a, es15.7, a)') 'alpha_a', ele0%a%alpha
write (*, '(a, es15.7, a)') 'alpha_b', ele0%b%alpha
write (*, '(a, es15.7, a)') 'eta_x',ele0%x%eta
write (*, '(a, es15.7, a)') 'etap_x',ele0%x%etap
write (*, '(a, es15.7, a)') 'sigma_x', sigma_x, ' m'
write (*, '(a, es15.7, a)') 'sigma_y', sigma_y, ' m'
write (*, '(a, es15.7, a)') 'sigma_z', sigma_z, ' m'
write (*, '(a, es15.7, a)') 'sigma_px', sigma_px, ' '
write (*, '(a, es15.7, a)') 'sigma_py', sigma_py, ' '
write (*, '(a, 3es15.7)') 'J damp     a, b, z:', normal_mode%a%j_damp, normal_mode%b%j_damp, normal_mode%z%j_damp
write (*, '(a, 3es15.7)') 'alpha damp a, b, z:', normal_mode%a%alpha_damp, normal_mode%b%alpha_damp, normal_mode%z%alpha_damp
write (*, '(a, 4es15.7)') 'I0, I1, I2, I3', normal_mode%synch_int(0:3)
write (*, '(a, es15.7, a)') 'eta_x ', ele0%x%eta, ' m'
write (*, '(a, es15.7, a)') 'etap_x ', ele0%x%etap, ''

!Qa = normal_mode%a%tune / (2.*pi)
!Qb = normal_mode%b%tune / (2.*pi) 
Qa = lat%a%tune/(2*pi)
Qb = lat%b%tune/(2*pi)
write (*, '(a, f10.5, a)') 'Qa', Qa, ' (lat%a%tune/(2*pi))'
write (*, '(a, f10.5, a)') 'Qb', Qb, ' (lat%b%tune/(2*pi))'

if (ptc_calc_on) then
  !print *, 'lat_to_ptc_layout (lat)'
  call lat_to_ptc_layout (lat)
  print *, 'ptc_emit_calc...'
  call ptc_emit_calc (lat%ele(0), normal_mode, sigma_mat, orb0)
  print *, 'ptc_emit_calc...done'
  write (*, '(a, es13.5, a)') 'ptc a emittance', normal_mode%a%emittance*1e9, ' nm'
  write (*, '(a, es13.5, a)') 'ptc b emittance', normal_mode%b%emittance*1e9, ' nm'

  print *, 'ptc_closed_orbit_calc'
  call reallocate_coord(ptc_closed_orb, lat)
  call  ptc_closed_orbit_calc (lat%branch(0), ptc_closed_orb)
  print *, 'ptc_emit_calc at ptc_closed_orb(0)'
  call ptc_emit_calc (lat%ele(0), normal_mode, sigma_mat, ptc_closed_orb(0))
  write (*, '(a, es13.5, a)') 'ptc a emittance', normal_mode%a%emittance*1e9, ' nm'
  write (*, '(a, es13.5, a)') 'ptc b emittance', normal_mode%b%emittance*1e9, ' nm'

endif



!--------------------
! Plotting
if (plot_on) then
  call qp_open_page ('X', plot_id, 2*500.0_rp, 500.0_rp, 'POINTS')
  call qp_set_page_border (0.05_rp, 0.05_rp, 0.05_rp, 0.05_rp, '%PAGE')
  call qp_set_margin (0.05_rp, 0.05_rp, 0.05_rp, 0.05_rp, '%PAGE')
  call qp_set_line_width_basic(0)
  !call qp_calc_and_set_axis ('X', -sigma_x*12, sigma_x*12, 4, 8, 'GENERAL')
  !call qp_calc_and_set_axis ('Y', 0.0_rp,  sigma_x*12, 4, 8, 'GENERAL')
  !call plscmap1n(256)
  !call plspal1('cmap1_blue_yellow.pal',1)
  !call qp_draw_axes ("x (m)", "y (m)")
  call qp_set_line_attrib ('PLOT', color =black$)
  call qp_set_line_width_basic(0)
  !call plot_wall(lat)
endif




! Parallel memory setup
if (parallel) then
 write(*,*) '_____________'
 write(*,*) 'Parallel mode'
  omp_n = 1
  !$ if (omp_on) then
  !$   omp_n = omp_get_max_threads()
  !$   print *, 'omp_get_max_threads(): ', omp_n
  !$ else
       print *, 'OpenMP disabled, using 1 thread'
  !$ endif
  
  !Set up independent lattices
  allocate(omp_lat(omp_n))
  allocate(omp_aperture(omp_n))
  allocate(omp_counter(omp_n))
  do i=1, omp_n
    omp_lat(i) = lat
    omp_counter(i) = 0
  end do
endif



if (amplitude_scan_on) then
print *, '_____________________________________________'
print *, 'Amplitude scan'
  call amplitude_scan(lat, orb0, normal_mode, plot_on, n_turns, nJa, nJb, nJz, step_a, step_b, step_z)
endif

if (frequency_drift_on) then
print *, '_____________________________________________'
print *, 'Frequency drift'



!orb0%vec(1) = sigma_x
!orb0%vec(6) = .01_rp
!call tunes_from_tracking(lat, orb0, n_turns, Qa, Qb)
!print *, 'Qa = ', Qa
!print *, 'Qb = ', Qb

!if(.false.) then
  !call init_coord(orb_mode, ele=ele0, element_end=upstream_end$)
  !orb_mode%vec=[sqrt(ele0%a%beta*emit_a), 0.0_rp, 0.0_rp , 0.0_rp, 0.0_rp, sigma_pz]
  !call convert_coords('MODE', orb_mode, ele0, 'LAB', orb0, err)
  !call action_angle_drift(lat, orb0, Qa0, Qb0, Qa1, Qb1, n_turns)
  !write (*, '(6f14.7)') Qa0, Qa1, Qb0, Qb1, abs(Qa0-Qa1), abs(Qb0-Qb1)

  amplitude_block%a_min = sqrt(ele0%a%beta*emit_a)
  amplitude_block%b_min = 0.01*sqrt(ele0%a%beta*emit_a)
  amplitude_block%pz_min = 0.0_rp 
  amplitude_block%a_max = 20*sqrt(ele0%a%beta*emit_a)
  amplitude_block%b_max = 0.01*10*sqrt(ele0%a%beta*emit_a)
  amplitude_block%pz_max = 0.01_rp
  amplitude_block%Na = 20
  amplitude_block%Nb = 2
  amplitude_block%Npz = 20
  call amplitude_block_init(lat, amplitude_block)
  !call write_amplitude_block(6, amplitude_block)

if (plot_on) then
  call qp_set_box(1, 1, 2, 1)
  call qp_clear_box()
  call qp_calc_and_set_axis ('X', 0.0_rp, 0.01_rp, 4, 8, 'GENERAL')
  call qp_calc_and_set_axis ('Y', 0.0_rp, 0.01_rp, 4, 8, 'GENERAL')
  call qp_draw_axes ('x (m)', '\gd\d0\u')    
  call qp_set_box(2, 1, 2, 1)
  call qp_clear_box()
  call qp_calc_and_set_axis ('X', 0.0_rp, 0.5_rp, 4, 8, 'GENERAL')
  call qp_calc_and_set_axis ('Y', 0.0_rp, 0.5_rp, 4, 8, 'GENERAL')
  call qp_draw_axes ('Q_a', 'Q_b')    
endif

 call frequency_drift(lat, amplitude_block, n_turns, plot_on)
endif

  
!endif

if (angle_scan_on) then
	print *, '_____________________________________________'
	print *, 'Angle Scan'
	print *, 'n_angles: ', n_angles
	print *, 'theta min: ', theta_min
	print *, 'theta max: ', theta_max

	track_input%n_turn = n_turns
	track_input%x_init = 0.005 !sigma_x/2
	track_input%y_init = .01 !100*sigma_y/2
	track_input%accuracy = .00001 
	theta_max = pi
	theta_min = 0.0_rp
	d_theta = (theta_max-theta_min)/(n_angles-1)

	allocate(angle_list(n_angles), aperture_list(n_angles))

	! Angle list setup
	do i=1, n_angles
	  angle_list(i) = (i-1)*d_theta
	end do

	if (plot_on) then
	  call qp_set_box(1, 1, 2, 1)
	  call qp_clear_box()
	  call qp_calc_and_set_axis ('X', -20*sigma_x, 20*sigma_x, 4, 8, 'GENERAL')
	  call qp_calc_and_set_axis ('Y',      0.0_rp, 10*sigma_x, 4, 8, 'GENERAL')
	  call qp_draw_axes ('x (m)', 'y (m)')    
	endif

	pz_init = 0.0_rp
	do i_pz = 0, 10
	  pz_init = i_pz*sigma_pz
	  print *, 'pz_init, pz_init/sigma_pz :', pz_init, i_pz
	  call dynamic_aperture1_parallel(lat, orb0, angle_list, track_input, aperture_list, pz_init)

	  if (plot_on) call qp_draw_polyline(aperture_list(:)%x, aperture_list(:)%y, &
						color = qp_continuous_color(i_pz/10.0_rp) )

	  ! Info
	  write(*, '(5a15, a10, x, a )') 'angle(deg) ', 'r(mm) ', 'r/sigma_x', 'turn_lost ', 's_lost(m) ', 'ele_lost '
	  do i = 1, n_angles
		write(*, '(f15.4, es15.7, f10.3, i15, f10.4, x, a )') angle_list(i)*180/pi, &
		  1000*sqrt(aperture_list(i)%x**2 +  aperture_list(i)%y**2),  &
		  sqrt(aperture_list(i)%x**2 +  aperture_list(i)%y**2)/sigma_x,  &
		  aperture_list(i)%i_turn, lat%ele(aperture_list(i)%ix_lat)%s, trim(lat%ele(aperture_list(i)%ix_lat)%name)
	  end do     
	end do

endif


if (bunch_tracking_on) then
  beam_init%a_emit = emit_a
  beam_init%b_emit = emit_b
  beam_init%sig_z = sigma_z
  beam_init%sig_pz = sigma_pz
  beam_init%center = closed_orb(0)%vec
  call init_bunch_distribution (ele0, lat%param, beam_init, 0, bunch)
  call init_bunch_distribution (ele0, lat%param, beam_init, 0, bunch_end)
  !bunch%particle(:)%vec(1) =bunch%particle(:)%vec(1) + 0.002*ele0%a%eta
  !bunch%particle(:)%vec(2) =bunch%particle(:)%vec(2) + 0.002*ele0%a%etap
  !bunch%particle(:)%vec(6) = 0.002
  
  
  n_particle =  size(bunch%particle)
  print *, 'n_particle: ', n_particle
  if (plot_on) then
    call qp_clear_page ()
    counter = 0
  endif
  
  print *, '_____________________________________________'
  if (touschek_calc_on) then 
	  print *, 'Touschek calculations'
	  print *, 'normal_mode%pz_aperture: ', normal_mode%pz_aperture
	  normal_mode%pz_aperture = 0.02_rp
	  print *, 'manually setting 2% normal_mode%pz_aperture: ', normal_mode%pz_aperture
	  print *, 'bunch charge (C) :',  bunch%charge_tot
	  lat%param%n_part= abs(bunch%charge_tot/e_charge)
	  print *, 'setting lat%param%n_part: ', lat%param%n_part

	  call touschek_lifetime(normal_mode, ttouschek, lat)
	  print *, 'Touschek lifetime (hours): ', ttouschek/3600
  endif

  if (taylor_tracking_on) then
    print *, 'Calculating one-turn map...' 
    call transfer_map_calc (lat, t_map_1_turn, err, 0, 0, integrate = .false., one_turn = .true.)
    print *, 'Calculating one-turn map...Done' 
  endif

  write (*,*) '_____________________________________________'
  write (*,*) 'Tracking'
  write (*, '(a8, x, a8, 6a15)') 'turn', 'n_live ', 'emit_a(nm)', 'emit_b(nm)', &
        'sigma_x(m)', 'sigma_y(m)', 'sigma_z(m)', 'sigma_pz(1)'

 call calc_bunch_params (bunch, bunch_params, err)
      write (*, '(i8, x, i8,  2f15.7, 4es15.5)') 0,  bunch_params%n_particle_live, &
         bunch_params%a%norm_emit/rel_gamma*1e9, bunch_params%b%norm_emit/rel_gamma*1e9, &
         sqrt(bunch_params%sigma(1,1)), sqrt(bunch_params%sigma(3,3)), &
         sqrt(bunch_params%sigma(5,5)), sqrt(bunch_params%sigma(6,6))
      !print *, 'plot particle bounds', minval(x_arr), maxval(x_arr ), minval(y_arr), maxval(y_arr )

  do j = 1, n_turns
    !print *, 'Tracking turn: ', j
    !do i = 1, size(bunch%particle)
    !  call track_taylor (bunch%particle(i)%vec, t_map_1_turn, bunch_end%particle(i)%vec)
    !  bunch%particle(i)%vec = bunch_end%particle(i)%vec
    !end do
    
    if ((.not. parallel) .and. (.not. taylor_tracking_on) ) then
      ! Track the ring using normal tracking
      do i=1, lat%n_ele_track
        !print *, 'tracking ele: ', i
        call track1_bunch(bunch, lat, lat%ele(i), bunch_end, err)
        bunch = bunch_end
      end do
      
    else if (taylor_tracking_on) then
      !print *, 'Taylor tracking, turn: ', j
      !$OMP parallel &
      !$OMP default(private), &
      !$OMP shared(n_particle, t_map_1_turn, bunch, bunch_end)
      !$OMP do schedule(dynamic)
      do i = 1, n_particle
        call track_taylor (bunch%particle(i)%vec, t_map_1_turn, bunch_end%particle(i)%vec)
        bunch%particle(i)%vec = bunch_end%particle(i)%vec
      end do   
      !$OMP end do
      !$OMP end parallel
    else  
      !$OMP parallel &
      !$OMP default(private), &
      !$OMP shared(n_particle, omp_lat, bunch, bunch_end)
      !$OMP do schedule(dynamic)
      do i=1, n_particle
        omp_i = 1
        !$ omp_i = omp_get_thread_num()+1
        !print *, omp_i, i
        do k=1, omp_lat(omp_i)%n_ele_track
          call track1(bunch%particle(i),  omp_lat(omp_i)%ele(k), omp_lat(omp_i)%param, bunch_end%particle(i))
          bunch%particle(i) = bunch_end%particle(i)
        end do
      end do
      !$OMP end do
      !$OMP end parallel
      bunch = bunch_end
    endif
    
    counter = counter + 1
    !if (plot_on .and. (counter==plot_every_n)) then
      call calc_bunch_params (bunch, bunch_params, err)
      if (err) then
        print *, 'err: calc_bunch_params'
        call calc_bunch_params (bunch, bunch_params, err)
      endif 
      !print *, bunch_params

      write (*, '(i8, x, i8,  2f15.7, 4es15.5)') j,  bunch_params%n_particle_live, &
         bunch_params%a%norm_emit/rel_gamma*1e9, bunch_params%b%norm_emit/rel_gamma*1e9, &
         sqrt(bunch_params%sigma(1,1)), sqrt(bunch_params%sigma(3,3)), &
         sqrt(bunch_params%sigma(5,5)), sqrt(bunch_params%sigma(6,6))
      !print *, 'plot particle bounds', minval(x_arr), maxval(x_arr ), minval(y_arr), maxval(y_arr )
    
    if (plot_on .and. (counter==plot_every_n)) then
      counter = 0
      
      call qp_set_box(1, 1, 2, 1)
      call qp_clear_box()
      call plot_particle_phase_space(bunch%particle, 1, 2, 12*sigma_x, 12*sigma_px, 'x (m)', 'p\dx\u (rad)')

      call qp_set_box(2, 1, 2, 1)
      call qp_clear_box()
      call plot_particle_phase_space(bunch%particle, 5, 6, 12*sigma_z, 12*sigma_pz, 'z (m)', 'p\dz\u (1)')

    endif

  end do

endif



!Close output file
close(outfile)

print *, "Written: ", outfile_name

!--------------------
! Plotting control
if (plot_on) then
   write (*, '(a)', advance = 'NO') ' Hit any key to end program: '
   read (*, '(a)') lat_name
  call qp_close_page 
endif



end program