subroutine write_single_particle_tbt(turn,muons, Q, nturns, ref_time)
       use bmad
       use muon_mod
       use parameters_bmad
       use naff_mod
       implicit none
       type(muon_struct), allocatable :: muons(:)
       integer i, turn
       integer j
       integer, save :: lun19, lun
       integer nmult, nfft,n, nturns
       real(rp) Q(3)
       real(rp) ref_time
       real(rp) beta
       logical first/.true./
       real(rp) freqs(10), spin_tunes(10), betatron_tunes(10)
       real(rp) cyc_freq
       real(rp) vec(6), ScrossBeta(3), betahat(3)
       real(rp) px, py, pz
       real costheta, sintheta

       real(rp), allocatable, save :: t(:)

      complex(rp), allocatable, save :: cdataxz(:), cdataxy(:), cdatayz(:), cdata_orbx(:),  cdata_orby(:) , cavex, cdatacs(:)
      complex(rp), save :: amps(10)

  i=turn
       if(first)then
        lun19 = lunget()
        open(unit=lun19, file= trim(directory)//'/'//'single_particle_tbt.dat')
         write(lun19,'(a10,1x,7a14,12x,a12,12x,2a16)')' turn ','x','xp','y','yp','z','pz','t','polarization', 'Period',' beta*c*Delta t'
!        write(lun19,'(a10,1x,7a12)')' turn ','x','xp','y','yp','z','pz','t'
        first=.false.
        n=nturns
        allocate(t(-1:n), cdataxz(n), cdataxy(n), cdatayz(n), cdatacs(n))
        allocate (cdata_orbx(n), cdata_orby(n))
        t(-1:0) = 0.
     endif
!       write(lun19,'(i10,1x,7es12.4)') i, muons(1)%coord%vec, muons(1)%coord%t
     if(i<1)return
     t(i) = muons(1)%coord%t
     beta= muons(1)%coord%beta
     write(lun19,'(i10,1x,7es14.6,3es12.4, 2es16.8)')i,muons(1)%coord%vec, muons(1)%coord%t,  Q, t(i)-t(i-1), beta*c_light*(t(i)-t(i-1))


  vec = muons(1)%coord%vec
  cdata_orbx(i) = cmplx(vec(1),vec(2))
  cdataxz(i)=cmplx(Q(3),Q(1))

  cdata_orby(i) = cmplx(vec(3),vec(4))
  cdataxy(i)=cmplx(Q(1),Q(2))

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

    vec(1:6) = muons(1)%coord%vec(1:6)
    px = vec(2)
    py = vec(4)
    pz = sqrt((vec(6)+1)**2 - vec(2)**2-vec(4)**2)
    betahat(1:3)=(/px,py,pz/)/sqrt(px**2+py**2+pz**2)
    costheta = dot_product(Q,betahat) !=cos(theta)
    ScrossBeta = cross_product(Q,betahat) 
    sintheta = dot_product(ScrossBeta,ScrossBeta)
    cdatacs(i) = cmplx(costheta,sintheta)

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

  if(i == nturns) then !this is the last turn
     nmult = int(log(float(i))/log(2.0d0))
     nfft = 2**nmult
     !  betatron tunes

     cavex=sum(cdata_orbx(1:nfft))/nfft
     cdata_orbx(1:nfft) = cdata_orbx(1:nfft)-cavex
     open(unit=lun, file = trim(directory)//'/'//'betatron_tunes.dat')

     call naff(cdata_orbx(1:nfft), betatron_tunes,amps)

     freqs = betatron_tunes * cyc_freq

     write(lun,'(a,1x,i10,a,1x,es16.8)')' fft turns = ',nfft,' cyclotron frequency = ', cyc_freq
     print '(a,16es12.4)',' cyclotron frequency = ', cyc_freq
     print '(a)', 'x-plane'

     print  '(3es16.8)', (freqs(j), real(amps(j)), aimag(amps(j)), j=1,10)
     write(lun, '(a)') 'x-plane' 
     write(lun, '(3a16)')'betatron_tune*f_cyc', 'A_real','A_img'
     write(lun,  '(3es16.8)') (freqs(j), real(amps(j)), aimag(amps(j)), j=1,10)

     call naff(cdata_orby(1:nfft), betatron_tunes,amps)

     freqs = betatron_tunes * cyc_freq

     print '(a)', 'y-plane'
     print  '(3es16.8)', (freqs(j), real(amps(j)), aimag(amps(j)), j=1,10)
     write(lun, '(a)') 'y-plane'
     write(lun, '(3a16)')'betatron_tune*f_cyc', 'A_real','A_img'
     write(lun,  '(3es16.8)') (freqs(j), real(amps(j)), aimag(amps(j)), j=1,10)
     close(unit=lun)  !betatron_tunes.dat

! spin tunes
     
     open(unit=lun, file = trim(directory)//'/'//'spin_tunes.dat')

     call naff(cdataxz(1:nfft), spin_tunes,amps)
     freqs = spin_tunes * cyc_freq

     write(lun,'(a,1x,i10,a,1x,es16.8)')' fft turns = ',nfft,' cyclotron frequency = ', cyc_freq
     print '(a,16es12.4)',' cyclotron frequency = ', cyc_freq
     print '(a)', 'z-x plane'

     print  '(3es16.8)', (freqs(j), real(amps(j)), aimag(amps(j)), j=1,10)
     write(lun, '(a)') 'z-x plane' 
     write(lun, '(3a16)')'spin_tune*f_cyc', 'A_real','A_img'
     write(lun,  '(3es16.8)') (freqs(j), real(amps(j)), aimag(amps(j)),j=1,10)

     call naff(cdataxy(1:nfft), spin_tunes,amps)
     freqs = spin_tunes * cyc_freq

     print '(a)', 'x-y plane'
     print  '(3es16.8)', (freqs(j), real(amps(j)), aimag(amps(j)),j=1,10)
     write(lun, '(a)') 'x-y plane'
     write(lun, '(3a16)')'spin_tune*f_cyc', 'A_real','A_img'
     write(lun,  '(3es16.8)') (freqs(j), real(amps(j)), aimag(amps(j)),j=1,10)

     call naff(cdatayz(1:nfft), spin_tunes,amps)
     freqs = spin_tunes * cyc_freq

     print '(a)', 'y-z plane'
     print  '(3es16.8)', (freqs(j), real(amps(j)), aimag(amps(j)),j=1,10)
     write(lun, '(a)') 'y-z plane'
     write(lun, '(3a16)')'spin_tune*f_cyc', 'A_real','A_img'
     write(lun,  '(3es16.8)') (freqs(j), real(amps(j)), aimag(amps(j)),j=1,10)

     call naff(cdatacs(1:nfft), spin_tunes,amps)
     freqs = spin_tunes * cyc_freq

     print '(a)', 'c-s plane'
     print  '(3es16.8)', (freqs(j), real(amps(j)), aimag(amps(j)),j=1,10)
     write(lun, '(a)') 'c-s plane'
     write(lun, '(3a16)')'spin_tune*f_cyc', 'A_real','A_img'
     write(lun,  '(3es16.8)') (freqs(j), real(amps(j)), aimag(amps(j)),j=1,10)

     close(unit=lun)

   endif
  
        return
       end