subroutine qp_tune(lat, qpx,qpy, ok, element1, element2)

! Use xquneing 1 and 2 to set chromaticity to qpx and qpy
! If ok is true then chromaticity is set. 
! element1 and element2 are optional character strings
!  The default elements are raw_xquneing_1 and 2.
  
  use bmad_struct
  use bmad_interface
  use bookkeeper_mod
  
  implicit none

  type (lat_struct) lat
  type(all_pointer_struct) a_ptr1, a_ptr2
 
  real(rp) qpx, qpy, chrom_x, chrom_y, delta_e/1.e-4/
  real(rp) chrom_0x, chrom_0y, dq_dxq(2,2)
  real(rp) dxq/0.1/
  real(rp) xq1,xq2,det
  real(rp) xq1_0, xq2_0
  
  character*40 ele_name
  character*16, optional :: element1, element2
  character*16 ele_name_1, ele_name_2
  character*12 attr_name
  
  integer ix_qt1, ix_qt2, i
  integer ix
  
  logical ok, err
 
  ele_name_1 = 'RAW_XQUNEING_1'
  ele_name_2 = 'RAW_XQUNEING_2'

  if(present(element1))ele_name_1 = element1
  if(present(element2))ele_name_2 = element2
  call string_trim(ele_name_1, ele_name_1, ix)
  call str_upcase(ele_name_1, ele_name_1)
  call string_trim(ele_name_2, ele_name_2, ix)
  call str_upcase(ele_name_2, ele_name_2)
    
  ok = .true.
  ix_qt1=0
  ix_qt2=0
  i=1
  do while (ix_qt1 == 0 .or.  ix_qt2 == 0)
     ele_name = lat%ele(i)%name
     call string_trim(ele_name, ele_name, ix)
       
     if(ele_name == ele_name_1) ix_qt1=i
     if(ele_name == ele_name_2) ix_qt2=i
     if(i > lat%n_ele_max)then
        print *,' Cannot change chromaticity, no xquneing '
        return
     endif
     i=i+1
  end do 

  call chrom_calc(lat, delta_e, chrom_0x, chrom_0y)
!  print '(a29,4e12.4)',' xq1, xq2,  chrom_0x, chrom_0y ', xq1, xq2, chrom_0x, chrom_0y
  
  if(lat%ele(ix_qt1)%key == group$)then
     attr_name = 'COMMAND'
  elseif(lat%ele(ix_qt1)%key == sextupole$ .or. lat%ele(ix_qt1)%key == overlay$ )then
     attr_name = 'K2'
  endif
  
  call string_trim(attr_name, attr_name, ix)
  
  call pointer_to_attribute(lat%ele(ix_qt1), attr_name, .true., a_ptr1, err)
  call pointer_to_attribute(lat%ele(ix_qt2), attr_name, .true., a_ptr2, err)
  
  do i=1, 10

     xq1_0 = a_ptr1%r
     xq2_0 = a_ptr2%r

     a_ptr1%r = xq1_0 + dxq
     call attribute_set_bookkeeping(lat%ele(ix_qt1), attr_name, err, a_ptr1)     
     a_ptr2%r = xq2_0 
     call attribute_set_bookkeeping(lat%ele(ix_qt2), attr_name, err, a_ptr2)
     call chrom_calc(lat, delta_e, chrom_x, chrom_y)
!     print '(a29,4e12.4)',' xq1, xq2,  chrom_x, chrom_y ', xq1, xq2, chrom_x, chrom_y
     
     dq_dxq(1,1) = (chrom_x-chrom_0x)/dxq
     dq_dxq(2,1) = (chrom_y-chrom_0y)/dxq

     a_ptr1%r = xq1_0
     call attribute_set_bookkeeping(lat%ele(ix_qt1), attr_name, err, a_ptr1)     
     a_ptr2%r = xq2_0 + dxq
     call attribute_set_bookkeeping(lat%ele(ix_qt2), attr_name, err, a_ptr2)
     call chrom_calc(lat, delta_e, chrom_x, chrom_y)
!     print '(a29,4e12.4)',' xq1, xq2,  chrom_x, chrom_y ', xq1, xq2, chrom_x, chrom_y

     dq_dxq(1,2) = (chrom_x-chrom_0x)/dxq
     dq_dxq(2,2) = (chrom_y-chrom_0y)/dxq


     a_ptr2%r = xq2_0
     call attribute_set_bookkeeping(lat%ele(ix_qt2), attr_name, err, a_ptr2)

     det = dq_dxq(1,1)*dq_dxq(2,2) - dq_dxq(1,2) * dq_dxq(2,1)
     if(det == 0.)exit 

     xq1 = (dq_dxq(2,2) * (qpx-chrom_0x)  - dq_dxq(1,2)*(qpy-chrom_0y)  )/det
     xq2 = (-dq_dxq(2,1) * (qpx-chrom_0x)  + dq_dxq(1,1)*(qpy-chrom_0y)  )/det

     a_ptr1%r = xq1_0 + xq1
     call attribute_set_bookkeeping(lat%ele(ix_qt1), attr_name, err, a_ptr1)    
     a_ptr2%r = xq2_0 + xq2
     call attribute_set_bookkeeping(lat%ele(ix_qt2), attr_name, err, a_ptr2)
     call chrom_calc(lat, delta_e, chrom_0x, chrom_0y)
!     print '(a29,4e12.4)',' xq1, xq2,  chrom_0x, chrom_0y ', xq1, xq2, chrom_0x, chrom_0y
     
     if(abs(chrom_0x-qpx) < 0.1 .and. abs(chrom_0y-qpy)<0.1)then
        ok=.true.
        write(*,'(2(a16, a2, e12.4, x))') ele_name_1, ': ', xq1, ele_name_2, ': ',  xq2
        return
     endif
     
  end do
  ok=.false. 
  
  return

end subroutine qp_tune