program virtualmachine

   use bmad
   use mode3_mod
   use cesr_basic_mod


     implicit none

      type (lat_struct) lat
      type (ele_struct) ele 
      type (ele_pointer_struct), allocatable:: eles(:)
      type (coord_struct), allocatable :: orbit(:), co(:)
      type (normal_modes_struct) mode
      type (rad_int_all_ele_struct) rad_int
      type (cesr_struct) cesr
      type (all_pointer_struct) a_ptr      
       integer unit, number, m, i, ix, nargs,j
       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 num, dec
       integer ix_attrib, n,tot
       integer ix_lat
       integer n1, n2
       integer nloc
       integer ix_ele

       real(rp) initial_offset(1:6)
       real(rp) y_off, cutoff
       real(rp) epsilon, val, parm0
       real(rp) do_dp(1:400)
       real(rp), allocatable :: Jac(:,:)
       real(rp) value
       character(len=4) :: numstr
       character(len=4) :: decstr
       character(len=8) ::  numdec
       character(len=32) :: dec_format
       character*16 name

      character*140 lat_file
      character*120 line, last_line
      logical err_flag/.false./, excitation_and_damping, quad_offsets/.false./, err, compute_jacobian/.true./
      namelist /input/nturns,n_particle, excitation_and_damping, initial_offset,quad_offsets,y_off,cutoff,seed,num,dec, compute_jacobian

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

!Initialize random number generator that is used to distribute quad offsets
      call ran_seed_put(seed)
      call ran_seed_get(seed)
      print *,' seed from ran_seed_get =',seed


      nargs = cesr_iargc()
      if(nargs == 1)then
         call cesr_getarg(1,lat_file)
         print *, 'Using ', trim(lat_file)
       else
        lat_file = 'bmad.'
        print '(a,$)',' Lattice file name ? (default= bmad.) '
        !read(5,'(a)') line
        !call string_trim(line, line, ix)
        !lat_file = line
        !if(ix == 0) lat_file = 'bmad.'
        lat_file = 'cta_2085mev_xr40m_20101215_bs.lat'
        print *, ' lat_file = ', lat_file
      endif

     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)

       open(unit=13,file='bunch.dat')

       call bmad_parser (lat_file, lat)
       print *,' bmad_parser done'
       call bmad_to_cesr(lat,cesr,err_flag)
       call reallocate_coord (orbit, lat%n_ele_track)
       allocate(Jac(400,100))
       if(quad_offsets)call implement_quad_offsets(lat,y_off,cutoff,seed)
       call closed_orbit_calc(lat,orbit,6)
       call track_all(lat,orbit)
       call lat_make_mat6(lat, -1, ref_orb = orbit)
       do i=1,lat%n_ele_track
          write(14,'(i10,1x,7es12.4)')i,lat%ele(i)%s, orbit(i)%vec
       end do
       call twiss_at_start(lat)
       call twiss_propagate_all(lat)
       call track_all(lat,orbit)
       call radiation_integrals (lat, orbit, mode, ix_cache, 0, rad_int)
       call calc_z_tune(lat)

       write(15,'(8a12)')'s','x','y','de','betax','betay','etax','etay'
       do i=1,lat%n_ele_track
        write(15,'(8es12.4)')lat%ele(i)%s,orbit(i)%vec(1),orbit(i)%vec(3),orbit(i)%vec(6),lat%ele(i)%a%beta,&
                     lat%ele(i)%b%beta, lat%ele(i)%x%eta, lat%ele(i)%y%eta
       end do

       
       
       write(numstr, '(i4)') num
       write(decstr, '(i4)') dec
       dec_format = '(a,es' // trim(adjustl(numstr)) // '.' // trim(adjustl(decstr)) // ')'
       print *, dec_format
       print '(a)',lat_file
       print '(3'//dec_format//')','horizontal tune = ', lat%a%tune/twopi,' vertical tune = ', lat%b%tune/twopi,' synch tune = ',lat%z%tune/twopi
       print '(4'//dec_format//')',' beta_x = ',lat%ele(lat%n_ele_track)%a%beta,' alpha_x = ',lat%ele(lat%n_ele_track)%a%alpha, &
               ' beta_y = ',lat%ele(lat%n_ele_track)%b%beta,' alpha_y = ',lat%ele(lat%n_ele_track)%b%alpha


       print dec_format,' horizontal emittance = ',mode%a%emittance
       print dec_format,' vertical emittance = ',mode%b%emittance
       print dec_format,' bunch length = ',mode%sig_z
       print dec_format,' energy spread = ',mode%sigE_E

       call twiss_at_start(lat)
       call calc_z_tune(lat)

       name = 'V13W'
       call lat_ele_locator(name,lat,eles,nloc,err)
       ix_ele = eles(nloc)%ele%ix_ele
       print '(a,1x,i10,1x,a)',name,eles(nloc)%ele%ix_ele, lat%ele(eles(nloc)%ele%ix_ele)%name
     print '(3i10)',lat%ele(ix_ele)%n_slave, lat%ele(ix_ele)%ix1_slave, lat%ele(ix_ele)%ix2_slave
     print *,lat%ele(lat%ele(ix_ele)%ix1_slave)%name
     ix = attribute_index(lat%ele(lat%ele(ix_ele)%ix1_slave),'VKICK')
     call pointer_to_attribute (eles(nloc)%ele, name, .false., a_ptr, err, .false.)
     if (associated (a_ptr%r)) a_ptr%r = value
     call lattice_bookkeeper(lat,err_flag)

     pause
    if(compute_jacobian)then
       m=0       
       do i=0,100
        if(cesr%skew_quad(i)%ix_lat == 0 .and. cesr%v_steer(i)%ix_lat == 0)cycle
        if(cesr%skew_quad(i)%ix_lat /= 0)then
           ix_attrib = k1$
           ix_lat = cesr%skew_quad(i)%ix_lat
           m=m+1
           call compute_derivative(lat, ix_lat, ix_attrib, do_dp,n)
           Jac(1:n,m)= do_dp(1:n) 
           print '(i10,1x,a,1x,20es12.4)',m,lat%ele(ix_lat)%name,do_dp(1:20)
        endif
        if(cesr%v_steer(i)%ix_lat /= 0)then
           ix_attrib = vkick$
           ix_lat = cesr%v_steer(i)%ix_lat
           m=m+1
           call compute_derivative(lat, ix_lat, ix_attrib, do_dp,n)
           Jac(1:n,m)= do_dp(1:n) 
           print '(i10,1x,a,1x,20es12.4)',m,lat%ele(ix_lat)%name,do_dp(1:20)
        endif
      end do
      write(16, '(i5,1x,a,a,i5,1x,a,a)') m,'(variables)','X',n,'objectives', 'jacobian'
      do i=1,n
       write(16,'(85es12.4)')(Jac(i,j),j=1,m)
      end do
    endif
        close(13)
end