subroutine ADTS_tracking(lat, emitx, emity, npoints)
  use bmad
  use nr
  use nrutil, only: dpc

 implicit none
 type (lat_struct) lat
 type(coord_struct), allocatable :: orbit(:), co(:)
 type (coord_struct), allocatable ::  end_orb(:,:)
 real(rp) vec(6)
 real(rp) dx/0.001/
 real(rp) tune(1:3,-10:10)
 real(rp), allocatable :: p(:), x(:)
 real(rp) emitx, emity
 real(rp) tol
real(rp), allocatable :: ftx(:,:), phi_all(:,:), ft(:)
real(rp) pie
real(rp) a,b,c,d
real(rp) sum_yx2(1:3),sum_x4(1:3),sum_x2(1:3),sum_y(1:3),sum_n(1:3)
real(rp) M(2,2),V(2),W(2)
 integer i
 integer j
 integer, parameter :: nturns=1024
 integer imx(1),low,high
 integer ind
 integer lun
 integer npoints
 complex(dpc),allocatable :: cdata(:,:)
interface
 subroutine fit_cos(tune_0, x,p, tol)
  use precision_def

  implicit none
  real(rp) tune_0, tune, tol
  real(rp), allocatable :: x(:), p(:)
  end subroutine
 end interface

 logical err_flag

 call reallocate_coord(orbit,lat%n_ele_max)
 call reallocate_coord(co,lat%n_ele_max)
 allocate(end_orb(-npoints:npoints,0:nturns))
 allocate(cdata(1,0:nturns))
 allocate(ftx(-npoints:npoints,nturns),phi_all(-npoints:npoints,nturns))
  allocate(ft(nturns))
 call closed_orbit_calc(lat,orbit,i_dim=4)
 
 allocate(p(3),x(0:nturns))
 sum_yx2(1:3)=0; sum_x4(1:3)=0; sum_x2(1:3)=0; sum_y(1:3)=0; sum_n(1:3)=0
 pie = twopi/2.
 do ind =1,3,2
 dx = sqrt(emitx * lat%ele(0)%a%beta)
 if(ind == 3)dx = sqrt(emitx * lat%ele(0)%a%beta)/2.
 tol = dx/10.
 low = .5*nturns
 high = 0.8*nturns

 do i= -npoints, npoints

  ft(:)=0
  vec(:)=0
  vec(ind) = i*dx
  co(0)%vec(1:6) = orbit(0)%vec(1:6) + vec(1:6)
  end_orb(i,0)%vec(ind:ind+1) = co(0)%vec(ind:ind+1)
  x(0) = end_orb(i,0)%vec(ind) 
  cdata(1,0)=cmplx(end_orb(i,0)%vec(ind), end_orb(i,0)%vec(ind+1))*(2*(sin((pie*0)/nturns))*(sin((pie*0)/nturns)))
  do j=1, nturns
   call track_all(lat, co)
   co(0) = co(lat%n_ele_track)
   end_orb(i,j)%vec(ind:ind+1) = co(0)%vec(ind:ind+1)
   x(j) = end_orb(i,j)%vec(ind)
   cdata(1,j)=cmplx(end_orb(i,j)%vec(ind), end_orb(i,j)%vec(ind+1))*(2*(sin((pie*j)/nturns))*(sin((pie*j)/nturns)))
 end do
 
! call fit_cos(lat%a%tune,x,p, tol)
 call fourrow(cdata(:,1:nturns),1)
 forall(j=1:nturns) ft(j)=sqrt(cdata(1,j)* conjg(cdata(1,j)))
  forall(j=1:nturns) ftx(i,j)=sqrt(cdata(1,j)* conjg(cdata(1,j)))
! forall(j=1:nturns) phi_all(i.j)=atan2(aimag(cdata(1,j)),real(cdata(1,j)))
  imx = maxloc(ft,1)
  d=ft(imx(1))
  b=ft(imx(1)+1)
    if(d /= 0. .or. b /= 0.)then
      c=cos(twopi/nturns)
      A= (-(d+b*c)*(d-b)+b*sqrt(c*c*(d+b)*(d+b)-2*d*b*(2*c*c-c-1)))/(d*d+b*b+2*d*b*c)
    endif
   tune(ind,i) = imx(1)/float(nturns)+(1/twopi)*asin(A*sin(twopi/float(nturns)))
   if(i == 0)then
     if(ind == 1)tune(ind,i) = lat%a%tune/twopi
     if(ind == 3)tune(ind,i) = lat%b%tune/twopi
   endif
   print '(5i10,es14.6)',ind,i,low,high,imx(1), tune(ind,i)
! fit a quadratic function
  if(i/= 0)then
   sum_yx2(ind) = sum_yx2(ind) + tune(ind,i)*end_orb(i,0)%vec(ind)**2
   sum_x4(ind) = sum_x4(ind) + end_orb(i,0)%vec(ind)**4
   sum_x2(ind) = sum_x2(ind) + end_orb(i,0)%vec(ind)**2
   sum_y(ind) = sum_y(ind) + tune(ind,i)
   sum_n(ind) = sum_n(ind) + 1
  endif
end do
M(2,2)  = sum_x4(ind)
M(1,2) = -sum_x2(ind)
M(2,1) = -sum_x2(ind)
M(1,1) = sum_n(ind)
V(1) = sum_yx2(ind)
V(2) = sum_y(ind)
W =matmul(M,V)/(M(1,1)*M(2,2) - M(1,2)*M(2,1))
if(ind == 1)print '(a,es12.4)','Quadratic dependence horizontal tune = ', W(1)
if(ind == 3)print '(a,es12.4)','Quadratic dependence vertical tune = ', W(1)
end do

lun = lunget()
open(unit=lun,file='tune_vs_amp.dat')
write(lun,'(4a14)')'x amp','x Tune ','y amp', 'y Tune'
do i=-npoints,npoints
 if(i == 0)cycle
 write(lun,'(4es14.6)')end_orb(i,0)%vec(1),tune(1,i),end_orb(i,0)%vec(3),tune(3,i)
end do
close(unit=lun)
print *,' write file tune_vs_amp.dat'
!write(11+ind-1,'(a10,1x,21es12.4)')'# Tune ', (tune(i),i=-10,10)
! write(13+ind-1,'(a10,1x,21es12.4)')'# Tune ', (tune(i),i=-10,10) 
!write(13,'(a10,1x,21i12)')'# Tune ', (imx(i),i=-10,10) 
! do j=0,nturns
! write(11+ind-1,'(i10,1x,21es12.4)') j, (end_orb(i,j)%vec(1),i=-10,10)
! write(13+ind-1,'(i10,1x,21es12.4)') j, (ftx(i,j),i=-10,10) 
!end do

 deallocate(p)
 return 
 end 

subroutine fit_cos(tune_0, x,p, tol)
use precision_def

implicit none
real(rp) tune_0, tune, tol
real(rp), allocatable :: x(:), p(:)
real(rp) c, chi
real(rp) gp(3), chi_last
integer i
integer N, iter
integer sign

interface
 function chi2(p,x)
  use precision_def
  implicit none
  real(rp), allocatable :: p(:), x(:)
  real(rp) chi2
 end function chi2

 function gradient(p,x)
  use precision_def
  implicit none

  real(rp), allocatable :: p(:), x(:)
  real(rp) gradient(size(p))
  real(rp) c2(-1:1), deltap(size(p))
 end function gradient
end interface

p(1) = x(0)
p(2) = tune_0
p(3) = 0.
N= size(x)
print *,' N = ', N
iter=0
chi=1.
sign = 1
chi_last = chi+0.1
do while(chi > tol .and. iter < 100 .and. chi < chi_last)
  chi_last = chi
   iter =iter + 1
  gp = gradient(p,x)

!  print '(i10,1x,7es12.4)',iter,p, gp, c
  c = -chi2(p,x)/dot_product(gp,gp)/iter
  p = p + c * gp
  chi = sqrt(chi2(p,x)/N)
  print '(i10,1x,6es12.4)',iter,c,p,chi2(p,x), chi
end do
if(chi > chi_last)then
    p=p-c*gp
    chi = sqrt(chi2(p,x)/N)
 endif
 print '(a10,i10,1x,6es12.4)','done',iter,c,p,chi2(p,x), chi
 return
end


function chi2(p,x)
use precision_def
implicit none
real(rp), allocatable :: p(:), x(:)
real(rp) chi2
integer i
chi2=0.
do i=0,size(x)-1
 chi2 = chi2 + (p(1)*cos(i*p(2)+p(3)) - x(i))**2
end do
return
end

function gradient(p,x)
use precision_def
implicit none

real(rp), allocatable :: p(:), x(:)
real(rp) gradient(size(p)), p0(size(p))
real(rp) c2(-1:1), deltap(size(p0))
integer i,j

interface
 function chi2(p,x)
  use precision_def
  implicit none
  real(rp), allocatable :: p(:), x(:)
  real(rp) chi2
 end function chi2
end interface


deltap(:) = 1.e-4
p0(:) = p(:)
do j=1,size(p)
 do i=-1,1,2
  p(j) = p0(j) + deltap(j) * i
  c2(i) = chi2(p,x)
 end do
 gradient(j) = (c2(1)-c2(-1))/2/deltap(j)
end do
p(:) = p0(:)
return
end