program capture_single_trajectory
 use sim_utils
 use precision_def
 implicit none

  real(rp) Bpeak, xinf/0.077/,xpinf/0./, Aperture,momentum/3.095/, index/0.11/, R0/7.112/, length_kick/3.81/
  real(rp) eta, beta, xtune, k_min, B_to_kick, Bmin, Bmax,B, dB
  real(rp) kicker_ends(6)/8.4315,9.7015,10.292,11.562,12.152,13.422/
  real(rp) kicker_start_ends(3)/8.4315,10.292,12.152/
  real(rp) kicker_length/1.27/
  real(rp) kpeak
  real(rp), allocatable :: applied_kick(:), kick(:), phi_kick(:)
  real(rp) ds, amp, delta/0.001/
  real(rp) phi, s, x
  real(rp) delta_max/0.006/
  real(rp) xmax
  real(rp) sum_kicks
  
  integer i,j,k,tot_individual_kicks,kicks_per_section/3/, n_steps/1000/
  integer nturns/10/, m_step/100/
  integer m
  integer lun, ios

  namelist/capture_single_trajectory_input/index,xinf,xpinf, Bpeak

    lun=lunget()
   open(unit=lun, file='capture_single_trajectory_input.dat', STATUS='old', IOSTAT=ios)
   READ(lun,NML=capture_single_trajectory_input, IOSTAT=ios)
   write(6, nml=capture_single_trajectory_input)
   close(lun)
   
  tot_individual_kicks = 3 * kicks_per_section
allocate(applied_kick(1:tot_individual_kicks))
allocate(kick(1:tot_individual_kicks),phi_kick(1:tot_individual_kicks))
xtune=sqrt(1.-index)
eta=R0/(1-index)
beta=R0/sqrt(1.-index)

write(18, '(a,es12.4)') ' beta  = ', beta
write(18, '(a,es12.4)')' eta  = ', eta
write(18, '(a,es12.4)')' xtune      = ', xtune
write(18, '(a,es12.4)')' Bpeak= ',Bpeak
write(18, '(a,es12.4)')' xinf= ', xinf*1000
write(18, '(a,es12.4)')' xpinf= ', xpinf*1000


B_to_kick = length_kick *0.3/momentum/1.e4*beta  ! kick[G] X  (B_to_kick) = kick[rad]*beta 


amp = 1/sqrt(beta)*((xinf-eta*delta)**2+beta*xpinf**2)**0.5

  k=0
if(kicks_per_section>0)then
  do i = 1,3  !starting end
   do k= 1,kicks_per_section
    applied_kick(k+(i-1)*kicks_per_section) = kicker_start_ends(i) + k * kicker_length/(kicks_per_section+1)
    print '(i10,2es12.4)',k+(i-1)*kicks_per_section , kicker_start_ends(i),applied_kick(k+(i-1)*kicks_per_section)
   end do
end do
    print '(a,es12.4)','cm kick at ', sum(applied_kick(:))/tot_individual_kicks/twopi/R0
else
!   applied_kick(1) =
endif
dB = 10. !G
do i= -10,10,1
B= Bpeak + i*dB   
kpeak= B/beta*B_to_kick

ds= twopi*R0/n_steps

do m=1,m_step
   kick(1:tot_individual_kicks)=0
   phi_kick(1:tot_individual_kicks)=0
   xmax=0
   delta = m*delta_max/m_step
   k=1
   do j=0,n_steps*nturns
      s= twopi*R0/n_steps * j
      phi = xtune/R0 * s
      x = eta*delta + (xinf-eta*delta)*cos(phi)+beta*xpinf*sin(phi)  + beta*sum( kick(1:tot_individual_kicks) * sin(phi-phi_kick(1:tot_individual_kicks)) )
      sum_kicks = sum(kick(1:tot_individual_kicks))
      if( minval(abs(applied_kick(k:tot_individual_kicks)-s)) <ds/2 .and. k<=tot_individual_kicks)then
         kick(k) = kpeak/tot_individual_kicks
         phi_kick(k) = phi
         k=k+1
      endif
      if(s>twopi*R0/2)xmax=max(xmax,x)
      write(11,'(es12.4,i10,2es12.4)')delta,j,phi,x
   end do
   write(11,'(/,a1,/)')'#'
   write(12,'(4es12.4)')B,delta,xmax, xmax-eta*delta
   if(delta == 0.)   write(13,'(4es12.4)')B,delta,xmax, xmax-eta*delta
end do
end do

end