!........................................................................
!
! Subroutine : CURLY_D (LAT, CO, QUAD, synch_int2, curlyd)
!
! Description: Subroutine to compute curly D for pretzel
!
! Arguments  : Input:  LAT, QUAD, CO(0:*)
!              Output: PRETZ.curly_d    !change in curly d with pretzel
!
! Mod/Commons:    
!
! Calls      :      
!
! Author     : 
!
! Modified   :     
!             
!........................................................................
!
!
! $Log$
! Revision 1.12  2007/01/30 16:15:13  dcs
! merged with branch_bmad_1.
!
! Revision 1.8  2006/02/08 19:09:46  mjf7
! Replaced obsolete bmad subroutine calls with new versions, and found variables which needed a save attribute.
!
! Revision 1.7  2004/11/08 19:16:27  dlr
! replace n_ele_maxx with n_ele_max
!
! Revision 1.6  2003/07/22 12:12:18  mjf7
! Synchronized some orbit allocation sizes.   - mjf
!
! Revision 1.5  2003/07/17 20:46:50  mjf7
! Fixed a bug which occurred at the end of an optimization loop.
! Changed many allocatable arrays to have the save parameter for faster
! performance.   - mjf
!
! Revision 1.4  2003/07/08 19:27:50  mjf7
!
!
! Modified all subroutines to correctly use allocatable lat elements. n_ele_maxx is no longer global, but a member variable of the lat struct. -mjf
!
! Revision 1.3  2003/06/05 18:33:27  cesrulib
! synch with bmad union removal
!
! Revision 1.2  2003/04/30 17:14:48  cesrulib
! dlr's changes since last import
!
! Revision 1.1.1.1  2002/12/13 19:23:27  cesrulib
! import bmadz
!
!
!........................................................................
!

subroutine curly_d(lat, co, quad, synch_int2, curlyd)

  use bmadz_mod
  use bmad_interface
  use bmadz_interface, except_dummy => curly_d
  use zquad_lens_mod
  implicit none

  type (lat_struct) lat
  type (coord_struct), allocatable ::  co(:)
  type (coord_struct), allocatable, save :: co1(:), co2(:)
  type (zquad_struct) quad
  integer i,j
  real(rp) quad_k, quad_l, G, int_num
  real(rp) ave_eta, delta
  real(rp) synch_int2, curlyd
  real(rp) dcd

  call reallocate_coord( co1, lat%n_ele_max )
  call reallocate_coord( co2, lat%n_ele_max )

  int_num = 0.

  co1(0)%vec(6) = 0.0001
  call closed_orbit_calc( lat, co1, 4 )
  call track_all(lat, co1)

  co2(0)%vec(6) = 0.00
  call closed_orbit_calc( lat, co2, 4 )
  call track_all(lat, co2)

  delta = co1(0)%vec(6)-co2(0)%vec(6)


!  write(11,*)'  index       quad name       ave_eta         delta D'
  do j=1,quad%n
    i=quad%lens(j)%ix
    quad_k = lat%ele(i)%value(k1$)
    quad_l = lat%ele(i)%value(l$)
    if(quad_k /= 0. .and. co(i)%vec(1) /= 0.)then
     G = quad_k * (co(i)%vec(1) + co(i-1)%vec(1))/2
     ave_eta = (co1(i)%vec(1) + co1(i-1)%vec(1) - &
          co2(i)%vec(1) - co2(i-1)%vec(1))/2/delta
     Int_num = Int_num + ave_eta * G * (G*G + 2 * quad_k) * quad_l
     dcd = ave_eta * G * (G*G + 2 * quad_k) * quad_l
!     write(11,1)j,lat%ele(i)%name, ave_eta, dcd
1    format(1x,i4,5x,a12,f12.6,f12.6)
    endif
  end do

!  write(11,*)' index    name       average eta         x'
  do i=1,lat%n_ele_track
   if(lat%ele(i)%key == marker$)then
     ave_eta = (co1(i)%vec(1) - co2(i)%vec(1))/delta
!     write(11,'(1x,i4,5x,a12,2f12.6)')i,lat%ele(i)%name, ave_eta, co(i).x.pos
   endif
  end do

  curlyd = Int_num/synch_int2

  return
end subroutine curly_d