subroutine compute_derivative(lat,name,attr_name,do_dp,n,epsilon,order,output_bool)
  use bmad
  use mode3_mod
  use sim_utils
  
    implicit none
  
    type (lat_struct) lat
    type (coord_struct), allocatable :: orbit(:)
    type (normal_modes_struct) mode
    type (rad_int_all_ele_struct) rad_int
    type(ele_pointer_struct), allocatable :: eles(:)

    type (all_pointer_struct) a_ptr
    
    integer j,m
    integer ix_cache/0/
    integer n,i
    integer ix_lat, ix_attrib
    integer tot, k, n_loc
    integer order, ix_len
    
    real(rp) epsilon
    real(rp) val, param0
    real(rp) do_dp(1:400)
    real(rp), allocatable:: objective(:,:)
    
    character*16 name, attr_name
    logical east, west, err
    logical, dimension(8) :: output_bool
    
    call reallocate_coord (orbit, lat%n_ele_track)
    allocate(objective(1:700, -order/2:order/2))
    
    !print '(a)',lat%ele(ix_lat)%name
    call lat_ele_locator(name, lat, eles, n_loc, err)
    ix_lat = eles(1)%ele%ix_ele
    call pointer_to_attribute(eles(1)%ele,attr_name,.false.,a_ptr,err,.false.)
    param0 = a_ptr%r
    do j = int(-order/2),int(order/2)
       if (j/=0) then
          n=0
          val = param0 + j* epsilon
          a_ptr%r = val
          call lat_make_mat6(lat, ix_lat)
          call closed_orbit_calc(lat,orbit,6)
          call track_all(lat,orbit)
          call lat_make_mat6(lat, -1, ref_orb = orbit)
          call twiss_at_start(lat)
          print '(/,a,es12.4)','Delta param =',j*epsilon
          print '(a14,6es12.4)','orbit(0)%vec =',orbit(0)%vec
          call twiss_propagate_all(lat)
          print '(a18,es12.4)', 'ele(0)%cmat(1,2) =',lat%ele(0)%c_mat(1,2)      !coupling
          call radiation_integrals (lat, orbit, mode, ix_cache, 0, rad_int)
          print '(a,es12.4,a,es12.4)','mode%a%emittance = ', mode%a%emittance,'mode%b%emittance =',mode%b%emittance
          do k=1,lat%n_ele_track
             west=.false.
             east=.false.
             if(index(lat%ele(k)%name, 'DET') ==0)cycle
             ix_len = len(lat%ele(k)%name)
             if(lat%ele(k)%name(ix_len:ix_len) ==  'W')west=.true.
             if(index(lat%ele(k)%name, 'E', back=.true.) == 7)east=.true.
             if(west .or. east)then
                !n=n+1
                !objective(n,j) = lat%ele(k)%c_mat(1,2)      !coupling
                !print '(2i10,1x,a,es12.4)',n,j,lat%ele(k)%name,objective(n,j)
                if(output_bool(1))then
                   n=n+1
                   objective(n,j) = orbit(k)%vec(1)             !x offset
                   !print *, 'x_off'
                endif
                if(output_bool(2))then
                   n=n+1
                   objective(n,j) = orbit(k)%vec(3)             !y offset
                   !print *, 'y_off'
                endif
                if(output_bool(3))then
                   n=n+1
                   objective(n,j) = lat%ele(k)%a%beta           !beta x
                   !print *, 'beta_x'
                endif
                if(output_bool(4))then
                   n=n+1
                   objective(n,j) = lat%ele(k)%b%beta           !beta y
                   !print *, 'beta_y'
                endif
                if(output_bool(5))then
                   n=n+1
                   objective(n,j) = lat%ele(k)%a%eta            !horizontal dispersion
                   !print *, 'x_disp'
                endif
                if(output_bool(6))then
                   n=n+1
                   objective(n,j) = lat%ele(k)%b%eta            !vertical dispersion
                   !print *, 'y_disp'
                endif
             endif
          end do
          if(output_bool(7))then
             n=n+1
             objective(n,j) = mode%a%emittance
             !print *, 'x_emit'
          endif
          if(output_bool(8))then
             n=n+1
             objective(n,j) = mode%b%emittance
             !print *, 'y_emit'
          endif
          !                    print '(2i10,es12.4)',n,j,mode%b%emittance
       endif
    end do
    if(order==2)then
       do_dp(1:n) = (objective(1:n,1)-objective(1:n,-1))/(2*epsilon)
    else if(order==4)then
       do_dp(1:n) = (objective(1:n,-2)-8*objective(1:n,-1)+8*objective(1:n,1)-objective(1:n,2))/(12*epsilon)
    else if(order==6)then
       do_dp(1:n) = (-objective(1:n,-3)+9*objective(1:n,-2)-45*objective(1:n,-1)+45*objective(1:n,1)-9*objective(1:n,2)+objective(1:n,3))/(60*epsilon)
    endif
    ! (objective(1:n,1)-objective(1:n,-1))/(2*epsilon)
    a_ptr%r=param0
    call lat_make_mat6(lat, ix_lat)
    call lattice_bookkeeper(lat,err)
    return
  end subroutine compute_derivative