program virtualmachine

   use bmad
   use mode3_mod


     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
      
       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

       real(rp) initial_offset(1:6)
       real(rp) y_off, cutoff

       character(len=4) :: numstr
       character(len=4) :: decstr
       character(len=8) ::  numdec
       character(len=32) :: dec_format

      character*140 lat_file
      character*120 line, last_line
      logical err_flag/.false./, excitation_and_damping, quad_offsets/.false./

      namelist /input/nturns,n_particle, excitation_and_damping, initial_offset,quad_offsets,y_off,cutoff,seed,num,dec

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

      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 reallocate_coord (orbit, lat%n_ele_track)
       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)

        forall(i=1:6)orbit(0)%vec(i) = orbit(0)%vec(i) + initial_offset(i)

        if(excitation_and_damping)then
           bmad_com%radiation_damping_on = .true.
           bmad_com%radiation_fluctuations_on = .true.
           print '(a,l)',' radiation_damping_on = ',bmad_com%radiation_damping_on
           print '(a,l)',' radiation_fluctuations_on = ',bmad_com%radiation_fluctuations_on 
        endif

        ix=0
        do i=1,nturns
           call track_all(lat, orbit, err_flag = err_flag)
!        write(21,'(i10,6es12.4)')orbit(lat%n_ele_track)%vec(1:6)
        end do

        close(13)
end