program fieldmap_test
  use bmad
  use runge_kutta_mod
  use quad_scrape_parameters

  implicit none

  type (lat_struct), target:: lat
  type (coord_struct), allocatable :: co(:)
  type (em_field_struct) field

  integer pgopen, istat1, istat2
  integer i, j, k,l

  integer lun
  integer nargs, iargc, ios
  integer end, nturns
  integer n
  integer track_state
  integer nmuons
  integer ix
  integer ngrids
  integer nbranch

  real(rp) efield, bfield, scale_efield_bfield/1./
  real(rp) Q(3), deltae, deltax
  real(rp) gamma_magic,beta_magic, p0c_magic, amu_measured/0.00116592091/, amu
  character*120 line, lat_file, lat_file_name/''/
  logical scraping_on_save, err
 
  namelist/input/ lat_file_name, scraping_on, init_quad_focus, init_quad_steer, quad_ramp_start_time, quad_ramp_end_time        !turn on quad scraping, scraping focus fraction, scraping steering fraction, start time and end time

  OPEN (UNIT=5, FILE='fieldmap_test_input.dat', STATUS='old', IOSTAT=ios)
    READ (5, NML=input, IOSTAT=ios)
    WRITE(6,NML=input)
  print *, 'ios=', ios
   rewind(unit=5)
    READ (5, NML=input)
  CLOSE(5)

! bmad_com%spin_tracking_on = .true.

!nargs = command_argument_count()
!  if (nargs == 1) then
!    call get_command_argument(1, lat_file)
!    !print *, 'Using ', trim(lat_file)
!  else if(lat_file_name /= '')then
!     lat_file = lat_file_name
!  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.'
!  endif

lat_file = lat_file_name
    print *, ' lat_file = ', trim(lat_file)

 ! bmad_com%spin_tracking_on = .true.
  bmad_com%auto_bookkeeper=.false.
  bmad_com%min_ds_adaptive_tracking = 0.0
  bmad_com%rel_tol_adaptive_tracking = 1d-10        ! Runge-Kutta tracking relative tolerance.
  bmad_com%abs_tol_adaptive_tracking = 1d-12       ! Runge-Kutta tracking absolute tolerance.

  scraping_on_save = scraping_on
!  scraping_on = .false.
  call bmad_parser (lat_file, lat)
  call twiss_at_start(lat)
  print '(a19,2f14.4)',' Fractional Tune   ',lat%a%tune/twopi,lat%b%tune/twopi
  call reallocate_coord (co, lat%n_ele_max)
  call closed_orbit_calc(lat, co)
  open(unit=22, file = 'closed_orbit_fieldmap_test.dat')
  do i=0,lat%n_ele_track
   write(22, '(i10,1x,a16,1x,es12.4,1x,6es12.4)')i,lat%ele(i)%name,lat%ele(i)%s, co(i)%vec
  end do
  close(unit=22)

 call lat_make_mat6(lat,  -1, co)
 call twiss_at_start(lat)  
  print '(a19,2f14.4)',' Fractional Tune   ',lat%a%tune/twopi,lat%b%tune/twopi
!  co(0)%vec(3) = co(0)%vec(3) - 0.0025 

print*,' closed orbit written'
 scraping_on = scraping_on_save
  do i=1,1000
    call track_all(lat,co)
    do j=1,lat%n_ele_track
     write(21, '(2i10,es14.6,6es14.6)')i,j,(i-1)+lat%ele(j)%s/lat%ele(lat%n_ele_track)%s,co(j)%vec
    end do
     write(23, '(i10,1x,es12.4, es14.6,6es14.6)')i,i+lat%ele(j)%s/lat%ele(lat%n_ele_track)%s,co(lat%n_ele_track)%t,co(lat%n_ele_track)%vec
    co(0) = co(lat%n_ele_track)
 end do
print *,'1000 turns complete'

 nbranch=0
 call twiss_propagate_cm_general(lat,nbranch)

 do i=1,lat%n_ele_track
  if(associated(lat%ele(i)%grid_field))then
    ngrids = size(lat%ele(i)%grid_field)
    print '(a16,1x,i10,1x,6es12.4)',lat%ele(i)%name, ngrids, lat%ele(i)%grid_field(1:ngrids)%field_scale
  endif
  if(index(lat%ele(i)%name,'QUAD1_LONG')/= 0)print '(a,1x,es12.4)',lat%ele(i)%name, value_of_attribute(lat%ele(i), 'QUAD_VOLTAGE', err)
 end do

end