program simple_inj_trajectories
use bmad
implicit none
integer, parameter:: npoints=1000
 integer  lun, lun1
 integer i
 integer ios
 real(rp) xinj, xpinj,q_k, index, R0, delta, theta_k
 real(rp) beta, eta, qx, phi_k
 real(rp) s 
 real(rp) a,x,ds,phi
 character*100 outfile
 namelist/simp_inj_traj/xinj,xpinj,q_k,index,R0,delta,theta_k, outfile

lun=lunget()
 open(unit=lun,file='simp_inj_traj.dat', STATUS='old', IOSTAT=ios)


    READ (lun, NML=simp_inj_traj, IOSTAT=ios)
    WRITE(6,NML=simp_inj_traj)
  print *, 'ios=', ios
   rewind(unit=lun)
  CLOSE(lun)

   beta=R0/sqrt(1-index)
   eta = R0/(1-index)
   qx= R0/beta
   phi_k=q_k*twopi * qx

  lun1 = lunget()
  open(unit=lun1, file = 'simp_inj_traj_out.dat')

   write(lun1,'(a,es12.4)')'xinj=',xinj
   write(lun1,'(a,es12.4)')'xpinj=',xpinj
   write(lun1,'(a,es12.4)')'q_k=',q_k
   write(lun1,'(a,es12.4)')'index=',index
   write(lun1,'(a,es12.4)')'R0=',R0
   write(lun1,'(a,es12.4)')'delta=',delta
   write(lun1,'(a,es12.4)')'theta_k=',theta_k
   write(lun1,'(a,es12.4)')'phi_k=', q_k* 360 * qx
   write(lun1,'(a,a)')'outfile=', '"'//trim(outfile)//'"'
   close(lun1)

   lun1 = lunget()
   open(unit=lun1, file = 'simple_inj_trajectory.dat')



   print '(3(a,es12.4))','beta =',beta,'  eta=',eta,'  Qx=',qx

   a = sqrt((xinj-eta*delta)**2+(xpinj*beta)**2)/sqrt(beta)
   
   ds=twopi*R0/npoints
   do i=1, 2*npoints
    s = ds*i
    phi = s/beta
    x = (xinj-eta*delta)*cos(phi)+eta*delta +beta*xpinj*sin(phi)
    if(phi>phi_k)x = x + beta*theta_k *sin(phi-phi_k)

   write(lun1,'(i10,3es12.4)')i,s,phi,x
   end do
end