program calc_jacobian

   use bmad
   use mode3_mod
   use cesr_basic_mod
   use sim_utils

      implicit none

      type (lat_struct) lat
      type (ele_struct) ele 
      type (coord_struct), allocatable :: orbit(:), co(:)
      type (normal_modes_struct) mode
      type (rad_int_all_ele_struct) rad_int
      type (cesr_struct) cesr
      type (ele_pointer_struct), allocatable :: eles(:)
      type (all_pointer_struct) a_ptr
      
      integer unit, number, m, i, nargs,j
      integer ix, ix2
      integer ix_cache/0/
      integer lun
      integer nturns, n_bunch, n_particle
      integer ios
      integer mturn/1/
      integer seed/1/
      integer bunch
      integer k,l
      integer tot
      integer ix_lat
      integer n1, n2
      integer n_loc
      integer order
      integer  n, ix_ele, ix_attr, change, ele_key
      integer ix_att
      integer ixx
      
      real(rp) initial_offset(1:6)
      real(rp) y_off, cutoff
      real(rp) set_value
      real(rp) epsilon
      real(rp) epsilons(1:85)
      real(rp) val, parm0
      real(rp) do_dp(1:400)
      real(rp) obj(1:400)
      real(rp), allocatable :: Jac(:,:)
      real(rp) value
      
      real(rp) param_values(1:100)
      character (len=64), dimension(100) :: param_names
      character (len=64), dimension(100) :: attr_names
      character (len=8), dimension(8) :: output

      character*64 name, attr_name, numstr
      character*16 prnt_format
      character*64 lat_file, edit_file
      character*16 attribute, attrib_name
      character  line*120
      character*16 ele_name
      character*1 ans



      logical err_flag/.false./, excitation_and_damping, quad_offsets/.false./
      logical, dimension(8) :: output_bool
      logical jacobian/.true./, file_exists,  make_xfer_mat, free
      logical another_variable/.true./
      logical err
      
      
      !================================
      !          Read input
      !================================
      namelist /input/lat_file,edit_file,nturns,n_particle,excitation_and_damping,initial_offset,&
           quad_offsets,y_off,cutoff,seed,prnt_format,output,jacobian,order,epsilons

      bmad_com%radiation_damping_on = .false.
      bmad_com%radiation_fluctuations_on = .false.

      lun = lunget()
      open(unit=lun, file='input.dat', STATUS ='old')
      read(lun, nml=input, IOSTAT=ios)
      rewind(lun)
      read(lun, nml=input, IOSTAT=ios)
      close(unit=lun)
      write(6, nml=input)
      print *, ' lat_file = ', lat_file
      
              
      !================================
      !         Parse lattice
      !================================
      call bmad_parser(lat_file, lat)
      print *,' bmad_parser done'
      call bmad_to_cesr(lat,cesr,err_flag)

    !================================
    !         List skew quads
    !================================

      print *,'List of all SKEW QUADS'
      print '(a10,1x,a10,1x,a16,a12)','Index','ix_lat','Element Name', 'value'
      attr_name = 'K1'
      do i=1,120
       ix_lat = cesr%skew_quad(i)%ix_lat
      if(ix_lat /= 0)then
       call pointer_to_attribute(lat%ele(ix_lat), attr_name, .false., a_ptr, err_flag, .false.)
       print '(i10,1x,i10,1x,a16,es12.4)',i,ix_lat,lat%ele(ix_lat)%name, a_ptr%r
       endif
      end do

    
    !================================
    !         List vertical steerings
    !================================
      print *
      print *,'List of all vertical steerings'
      print '(a10,1x,a10,1x,a16,a12)','Index','ix_lat','Element Name', 'value'
      attr_name = 'VKICK'
      do i=1,120
       ix_lat = cesr%v_steer(i)%ix_lat
       if(ix_lat /= 0)then
        call pointer_to_attribute(lat%ele(ix_lat), attr_name, .false., a_ptr, err_flag, .false.)
        print '(i10,1x,i10,1x,a,es12.4)',i,ix_lat,lat%ele(ix_lat)%name, a_ptr%r
       endif
      end do

    !================================
    !         Read - parameter to vary, attribute, epsilon
    !================================

    i=0
    
    do while(another_variable)
     i=i+1

      print '(a81)',' Type element <name> ,<attribute>, and <epsilon> for computing derivatives'
      read(5,'(a,$)') line
      call str_upcase(line, line)
    print *, line
      call string_trim(line, line, ix)

      ele_name = line(1:ix)
      call lat_ele_locator(line(1:ix),lat,eles,n_loc,err)     
      if(err)then
       print *,' Cannot find element ',line(1:ix), ' Try again'
       cycle
      endif

     call string_trim(line(ix+1:), line, ix)
!second word is attribute
     attr_name = line(1:ix)
     call pointer_to_attribute(eles(1)%ele, trim(attribute), .false., a_ptr, err_flag, .false.)    

     call string_trim(line(ix+1:), line, ix)
!third word is epsilon
     read(line(1:ix),*)epsilon
     print '(a12,a16,a13,a16,a16,es12.4)',' Variable = ',ele_name,' Attribute = ', attr_name, ' Epsilon = ', epsilon

      !================================
      !      Calculate Jacobian
      !================================
      ! ~~ Calculate jacobian for desired outputs ~~ !
      if(jacobian)then
         if(.not. allocated(Jac))allocate(Jac(400,100))
         print *, 'Jacobian'
         output_bool = (/.false.,.false.,.false.,.false.,.false.,.false.,.false.,.false./)
         n=0
            do k=1,8
               if(output(k)=='x_off')then
                  output_bool(1)=.true.
               else if(output(k)=='y_off')then
                  output_bool(2)=.true.
               else if(output(k)=='beta_x')then
                  output_bool(3)=.true.
               else if(output(k)=='beta_y')then
                  output_bool(4)=.true.
               else if(output(k)=='x_disp')then
                  output_bool(5)=.true.
               else if(output(k)=='y_disp')then
                  output_bool(6)=.true.
               else if(output(k)=='x_emit')then
                  output_bool(7)=.true.
               else if(output(k)=='y_emit')then
                  output_bool(8)=.true.
               endif
            end do

            call compute_derivative(lat,ele_name,attr_name,do_dp,n,epsilon,order,output_bool)
            Jac(1:n,i)=do_dp(1:n)
           endif ! compute jacobian
          print *
          print '(a,$)',' Try another variable ? (Y/N)'
          read(5, *)ans
          if(ans /= 'Y' .and. ans /= 'y')another_variable = .false.
         end do 
         m=i 
         print '(i5,1x,a,a,i5,1x,a,1x,a)', m,'(variables)','X',n,'objectives', 'jacobian'
         write(16, '(i5,1x,a,a,i5,1x,a,1x,a)') m,'(variables)','X',n,'objectives', 'jacobian'
         write(numstr, *) m
         do i=1,n
            write(16,'('//numstr//'es20.10e2)') (Jac(i,j),j=1,m)
         end do


end program calc_jacobian