program virtualmachine

   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 ix_attrib, n,tot
      integer ix_lat
      integer n1, n2
      integer n_loc
      integer order
      
      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) 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*120 line, last_line
      logical err_flag/.false./, excitation_and_damping, quad_offsets/.false./
      logical, dimension(8) :: output_bool
      logical jacobian/.true./, file_exists, err, make_xfer_mat, free
      
      
      !================================
      !          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.
      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)
      print *, ' lat_file = ', lat_file
      
      
      !================================
      !         Open files
      !================================
      !inquire(file='outputs',exist=file_exists)
      !if(file_exists)then
      call system('rm -r outputs')
      call system('mkdir outputs')
      !else
      !   call system('rm -r outputs')
      !   call system('mkdir outputs')
      !endif
      ! ~~ Open files for writing ~~ !
      open(unit=13, file='outputs/bunch.dat')
      open(unit=14, file='outputs/fort.14')
      open(unit=15, file='outputs/tracking.dat')
      open(unit=16, file='outputs/jacobian.dat')
      open(unit=17, file='outputs/detectors.dat')
      !open(unit=18, file='outputs/param_loc.dat')
      ! ~~ Read file ~~ !
      open(unit=51, file=edit_file, status='old')
      
      
      !================================
      !         Parse lattice
      !================================
      call bmad_parser(lat_file, lat)
      print *,' bmad_parser done'
      call bmad_to_cesr(lat,cesr,err_flag)


      !================================
      !         Write old inputs
      !================================
      !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
      !         write(19, '(a20,'//prnt_format//')') lat%ele(ix_lat)%name, lat%ele(ix_lat)%value(ix_attrib)
      !      endif
      !      if(cesr%v_steer(i)%ix_lat /= 0)then
      !         ix_attrib = vkick$
      !         ix_lat = cesr%v_steer(i)%ix_lat
      !         write(19, '(a20,'//prnt_format//')') lat%ele(ix_lat)%name, lat%ele(ix_lat)%value(ix_attrib)
      !      endif
      !end do
      
      
      !================================
      !         Read new inputs
      !================================
      m=0
      do
         ! ~~ Check line for appropriate format ~~ !
         read(51, '(a)', iostat=ios) line
         if (ios /= 0) exit         
         ix = index(line, '=')
         if (ix == 0) cycle   
         m=m+1
         ix2 = index(line, ',')
         
         ! ~~ Parse line ~~ !
         name = adjustl(line(:ix2-1))
         attr_name = adjustl(line(ix2+1:ix-1))         
         read(line(ix+1:), *) val
         
         ! ~~ Store for Jacobian ~~ !
         param_names(m) = name
         attr_names(m) = attr_name
         param_values(m) = val
         
         ! ~~ Update lattice ~~ !
         call lat_ele_locator(name, lat, eles, n_loc, err)
         ix_lat = eles(1)%ele%ix_ele
         !print *, lat%ele(ix_lat)%name, ix_lat, lat%ele(ix_lat)%value(k1$)
         call pointer_to_attribute(eles(1)%ele, attr_name, .false., a_ptr, err, .false.)
         a_ptr%r = val
         !print *, lat%ele(ix_lat)%name, ix_lat, lat%ele(ix_lat)%value(k1$)
      end do
      call lattice_bookkeeper(lat, err)
      !stop

      
      !================================
      !        Calculate data
      !================================
      call reallocate_coord (orbit, lat%n_ele_track)
      allocate(Jac(400,100))
      call ran_seed_put(seed)
      call ran_seed_get(seed)
      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 inputs
      !================================
      !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
      !         write(19, '(a20,'//prnt_format//')') lat%ele(ix_lat)%name, lat%ele(ix_lat)%value(ix_attrib)
      !      endif
      !      if(cesr%v_steer(i)%ix_lat /= 0)then
      !         ix_attrib = vkick$
      !         ix_lat = cesr%v_steer(i)%ix_lat
      !         write(19, '(a20,'//prnt_format//')') lat%ele(ix_lat)%name, lat%ele(ix_lat)%value(ix_attrib)
      !      endif
      !end do      
      
      
      !================================
      !         Write output
      !================================
      ! ~~ Tracking data ~~ !
      write(15,'(8a12)')'s','x','y','de','betax','betay','etax','etay'
      do i=1,lat%n_ele_track
         write(15,'(8'//prnt_format//')')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 inputs ~~ !
      !do i=1,m
      !   name = param_names(i)
      !   call lat_ele_locator(name, lat, eles, n_loc, err)
      !   ix_lat = eles(1)%ele%ix_ele
      !   write(18, '(a12,'//prnt_format//')') lat%ele(ix_lat)%name, lat%ele(ix_lat)%s
      !end do
      !do i=0,100
      !   !print '(2'//prnt_format//',3a)', cesr%skew_quad(i)%ix_lat, cesr%v_steer(i)%ix_lat, ' ', &
      !   !     lat%ele(cesr%skew_quad(i)%ix_lat)%name, lat%ele(cesr%v_steer(i)%ix_lat)%name
      !   write(18, '(2'//prnt_format//',3a)') cesr%skew_quad(i)%ix_lat, cesr%v_steer(i)%ix_lat, ' ', &
      !        lat%ele(cesr%skew_quad(i)%ix_lat)%name, lat%ele(cesr%v_steer(i)%ix_lat)%name
      !end do
      
      ! ~~ Detector outputs ~~ !
      ix = index(prnt_format, '.')
      numstr = prnt_format(3:ix-1)
      write(17,'(a12,a,7a'//numstr//')') 'name', ' ', 's','x_off','y_off','beta_x','beta_y','x_disp','y_disp'
      do k=1, lat%n_ele_track
         if (index(lat%ele(k)%name, 'DET') /= 0) then
            if (index(lat%ele(k)%name,'W',back=.true.) == 7 .or. &
                 index(lat%ele(k)%name,'E',back=.true.) == 7) then
               write(17, '(a12,a,7'//prnt_format//')') lat%ele(k)%name, '=', &
                    lat%ele(k)%s, orbit(k)%vec(1), orbit(k)%vec(3), &
                    lat%ele(k)%a%beta, lat%ele(k)%b%beta, lat%ele(k)%a%eta, lat%ele(k)%b%eta 
            endif
         endif
      end do
      write(17, '(a,'//prnt_format//')') 'x_emit = ',mode%a%emittance
      write(17, '(a,'//prnt_format//')') 'y_emit = ',mode%b%emittance
      
      ! ~~ Other outputs ~~ !
      print '(a)',lat_file
      print '(3(a,'//prnt_format//'))', 'horizontal tune = ', lat%a%tune/twopi,&
           ' vertical tune = ', lat%b%tune/twopi,' synch tune = ',lat%z%tune/twopi
      print '(4(a,'//prnt_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 '(a,'//prnt_format//')', 'horizontal emittance = ',mode%a%emittance
      print '(a,'//prnt_format//')', 'vertical emittance = ',mode%b%emittance
      print '(a,'//prnt_format//')', 'bunch length = ',mode%sig_z
      print '(a,'//prnt_format//')', 'energy spread = ',mode%sigE_E
      
      ! ~~ Bunch data ~~ !
      write(13,'(a,'//prnt_format//')') 'horizontal emittance = ',mode%a%emittance
      write(13,'(a,'//prnt_format//')') 'vertical emittance = ',mode%b%emittance
      write(13,'(a,'//prnt_format//')') 'bunch length = ',mode%sig_z
      write(13,'(a,'//prnt_format//')') 'energy spread = ',mode%sigE_E
      write(13,'(a,'//prnt_format//')') 'horizontal tune = ', lat%a%tune/twopi
      write(13,'(a,'//prnt_format//')') 'vertical tune = ', lat%b%tune/twopi
      write(13,'(a,'//prnt_format//')') 'synch tune = ', lat%z%tune/twopi
      write(13,'(a,'//prnt_format//')') 'beta_x = ', lat%ele(lat%n_ele_track)%a%beta
      write(13,'(a,'//prnt_format//')') 'alpha_x = ', lat%ele(lat%n_ele_track)%a%alpha
      write(13,'(a,'//prnt_format//')') 'beta_y = ', lat%ele(lat%n_ele_track)%b%beta
      write(13,'(a,'//prnt_format//')') 'alpha_y = ', lat%ele(lat%n_ele_track)%b%alpha
      write(13,'(a,'//prnt_format//')') 'horizontal emittance = ', mode%a%emittance
      write(13,'(a,'//prnt_format//')') 'vertical emittance = ', mode%b%emittance
      write(13,'(a,'//prnt_format//')') 'bunch length = ', mode%sig_z
      write(13,'(a,'//prnt_format//')') 'energy spread = ', mode%sigE_E
      
      
      !================================
      !            Other???
      !================================
      call twiss_at_start(lat)
      call calc_z_tune(lat)
      
      !================================
      !      Calculate Jacobian
      !================================
      ! ~~ Calculate jacobian for desired outputs ~~ !
      if(jacobian)then
         print *, 'Jacobian'
         output_bool = (/.false.,.false.,.false.,.false.,.false.,.false.,.false.,.false./)
         n=0
         do i=1,m
            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
            name = param_names(i)
            attr_name = attr_names(i)
            call compute_derivative(lat,name,attr_name,do_dp,n,epsilons(i),order,output_bool)
            Jac(1:n,i)=do_dp(1:n)
         end do
         print '(i5,1x,a,a,i5,1x,a,a)', m,'(variables)','X',n,'objectives', 'jacobian'
         write(16, '(i5,1x,a,a,i5,1x,a,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
      endif
      
      !================================
      !          Close files
      !================================
      close(13)
      close(14)
      close(15)
      close(16)
      close(17)
      !close(18)
      
      close(51)

end program virtualmachine