program test
  use bmad
  
  implicit none

  type (lat_struct) lat
  type (ele_struct) ele
  type (coord_struct), allocatable :: co(:)

  integer pgopen, istat1, istat2, nargs
  integer i, j,ix

  integer lun, track_state
  integer npoints
  logical ok
  character*120 line, lat_file
  real(rp) deltap/0.00005/
  real(rp) qx, qy
  real(rp) emitx,emity
  
      nargs = cesr_iargc()
      if(nargs == 1)then
         call cesr_getarg(1,lat_file)
         print *, 'Using ', trim(lat_file)
       else
        lat_file = 'bmad.'
        print '(a,$)',' Lattice file name ? (default= bmad.) '
        read(5,'(a)') line
        call string_trim(line, line, ix)
        lat_file = line
        if(ix == 0) lat_file = 'bmad.'
        print *, ' lat_file = ', lat_file
      endif
       call bmad_parser(lat_file, lat)
       call reallocate_coord (co, lat%n_ele_track)
        call closed_orbit_calc(lat,co,4)
          call lat_make_mat6(lat, ref_orb=co)
          call twiss_at_start(lat)
          print '(3es12.4)',co(0)%vec(6),lat%a%tune/twopi, lat%b%tune/twopi
          qx=lat%a%tune/twopi
          qy=lat%b%tune/twopi
          lun=lunget()
          open(unit=lun,file='tune_vs_momentum.dat')
          write(lun, '(5a12)')'momentum','h tune','v_tune','dQx','dQy'
       do i=-50,50
          co(0)%vec(6) = i*deltap
          call closed_orbit_calc(lat,co,4)
          call lat_make_mat6(lat, ref_orb=co)
          call twiss_at_start(lat)
          write(lun, '(5es12.4)')co(0)%vec(6),lat%a%tune/twopi, lat%b%tune/twopi,lat%a%tune/twopi-qx, lat%b%tune/twopi-qy
       end do
       print *,' write ', 'tune_vs_momentum.dat'
       emitx=1.e-6
       emity=1.e-6
       npoints=30
       call ADTS_tracking(lat, emitx, emity,npoints)
     end program test