program naff_spin
  use bmad
  use naff_mod
  implicit none

  real(rp) coord_vec(6), time, Q(3)
  real(rp) freqs(10), spin_tunes(10)
  real(rp) cyc_freq
  integer i, lun
  integer reason
  integer ix
  integer, parameter :: nturn=4096
  real(rp), dimension(nturn) :: t
  character*200 string, string_out
  complex(rp), dimension(nturn) :: cdataxz, cdataxy, cdatayz
  complex(rp) amps(1)

  
  lun=lunget()
  open(unit=lun,file = 'single_particle_tbt.dat')
i=0
  do while(.true.)
  read(lun,'(a)',iostat = reason )string
  if(reason < 0 .or. i>4096) exit
  call string_trim(string, string_out, ix)
  if(index(string_out(1:ix),'turn')/=0)cycle
  i=i+1
  read(string,'(i10,1x,7es14.6,3es12.4)')i,coord_vec, time,  Q
!  print '(i10,1x,7es14.6,3es12.4)', i,coord_vec, time,  Q
  t(i) = time
  cdataxz(i)=cmplx(Q(3),Q(1))
  
  cdataxy(i)=cmplx(Q(1),Q(2))

  cdatayz(i)=cmplx(Q(2),Q(3))

  if(i > 1)cyc_freq = (i-1)/(t(i)-t(1))   

end do

close(unit=lun)
open(unit=lun, file = 'spin_tunes.dat')

call naff(cdataxz, spin_tunes,amps)

freqs = spin_tunes * cyc_freq

print '(a,16es12.4)',' cyclotron frequency = ', cyc_freq
print '(a)', 'z-x plane'
print  '(3es16.8)', freqs, real(amps), aimag(amps)
write(lun, '(a)') 'z-x plane'
write(lun,  '(3es16.8)') freqs, real(amps), aimag(amps)

call naff(cdataxy, spin_tunes,amps)

freqs = spin_tunes * cyc_freq

print '(a)', 'x-y plane'
print  '(3es16.8)', freqs, real(amps), aimag(amps)
write(lun, '(a)') 'x-y plane'
write(lun,  '(3es16.8)') freqs, real(amps), aimag(amps)
call naff(cdatayz, spin_tunes,amps)

freqs = spin_tunes * cyc_freq

print '(a)', 'y-z plane'
print  '(3es16.8)', freqs, real(amps), aimag(amps)
write(lun, '(a)') 'y-z plane'
write(lun,  '(3es16.8)') freqs, real(amps), aimag(amps)
close(unit=lun)

end program naff_spin