program compute_injection_transfer_matrix
 use bmad
 implicit none
   type (lat_struct) lat
   type (ele_struct) ele_inf, ele_start
   type(coord_struct) start_orb, end_orb
   type(coord_struct), allocatable :: co(:)
   type (em_field_struct) field
   integer nargs, ix,i, ix_ele
   integer track_state
   integer lun
   integer ios
   integer ix_end
   character*60 lat_file, line
   logical err_flag
   logical local_ref/.false./, calcd/.false./
   real(rp) max_beta, min_max_beta,dist
   real(rp) s_start/0./,s_end/0./
   real(rp) energy
   real(rp) s_rel, t_rel
   namelist /inj_input/start_orb, energy

!  bmad_com%rel_tol_tracking           = 1.e-8
!  bmad_com%abs_tol_tracking           = 1.e-10
!  bmad_com%rel_tol_adaptive_tracking  = 1.e-10
  bmad_com%abs_tol_adaptive_tracking  = 1.e-8
!  bmad_com%min_ds_adaptive_tracking   = 1.e-14
!  bmad_com%fatal_ds_adaptive_tracking = 0.
     bmad_com%auto_bookkeeper=.false.
      nargs = command_argument_count()
      if(nargs == 1)then
         call get_command_argument(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.'
       print *, ' lat_file = ', lat_file
      endif

       bmad_com%min_ds_adaptive_tracking = 0.00001
       call bmad_parser (lat_file, lat, make_mats6=.true.)

     start_orb%vec(1:6)=0
     lun = lunget()
     open(unit=lun, file='inj_input.dat', STATUS ='old')
     read(lun, nml=inj_input, IOSTAT=ios)
     close(unit=lun)

       call reallocate_coord (co, lat%n_ele_max)


print *,' Just before track_all'

!       lat%ele(0)%value(E_TOT$) = energy * 1.e9
       lat%ele(1:lat%n_ele_track)%alias= 'Record'
!       call set_flags_for_changed_attribute(lat%ele(0), lat%ele(0)%value(E_TOT$))
!       call lattice_bookkeeper(lat, err_flag)
       print *,' Set energy [GeV] = ', lat%ele(0)%value(E_TOT$)/1.e9
        co(0) = start_orb
        call track_all(lat, co, 0, track_state, err_flag)
        do i = 0, lat%n_ele_track
          print '(i10,1x,a,1x,7es12.4)',i,lat%ele(i)%name,lat%ele(i)%s, co(i)%vec
        end do
       lat%ele(1:lat%n_ele_track)%alias= 'No'
       print *,' After track1 loop, track_state =', track_state, 'err_flag = ', err_flag

       print *,' Just before make_mat6'

       call lat_make_mat6 (lat, -1, co, err_flag=err_flag)

      print *,' Just after make_mat6'
    
       ele_start = lat%ele(0)
       min_max_beta = 100.

       call twiss_propagate_all(lat)
!       call element_locator('HOLE',lat,ix_ele)
!       dist = lat%ele(ix_ele)%s 

       max_beta = 0.      

        lat%ele(1:lat%n_ele_track)%alias= 'trajectory'
         print *,' Length = ', lat%ele(lat%n_ele_track)%s
         write(21,'(a12,1x,6a12,3a12)')'s', 'x','xp','y','yp','z','delta','Bx','By','Bz'
         do while(s_end < lat%ele(lat%n_ele_track)%s)

          call track_from_s_to_s(lat, s_start,s_end, start_orb, end_orb)
          ix_end = element_at_s (lat, s_end, .true.)
          s_rel = s_end - lat%ele(ix_end-1)%s
          call em_field_custom(lat%ele(ix_end), lat%param, s_rel,t_rel, end_orb, local_ref, field,calcd, err_flag)
          write(21,'(es12.4,1x,6es12.4,3es12.4)')s_end, end_orb%vec, field%B(1:3)
          call twiss_and_track_at_s(lat,s_end,ele_inf)
          write(22,'(a16,7es12.4)')ele_inf%name, ele_inf%s, ele_inf%a%beta, ele_inf%b%beta, ele_inf%a%alpha, ele_inf%b%alpha, &
                   ele_inf%x%eta,  ele_inf%x%etap
          s_end = s_end + 0.1
         end do
         close(unit=21)
         close(unit=22)
         print '(a)', ' fort.21 and fort.22 are written'

     lun = lunget()
     open(unit=lun, file='transfer_matrices_on_energy.dat')
     do i=1,lat%n_ele_track
       call type_ele (lat%ele(i), type_floor_coords=.true., nunit = lun)
     end do
     close(unit=lun)



         write(16,'(a,es12.4,2i4)')'Max beta = ',max_beta
         print '(a,es12.4)','Max beta = ',max_beta
!       min_max_beta = min(max_beta,min_max_beta)
!        endif

!      end do
!     end do
!      endif
     print *,' min_max_beta = ',min_max_beta
       end