subroutine injection_channel(lat,spin_tracking_on, from_orbit, circumference, make_movie)
         use bmad
         use parameters_bmad
         use muon_interface, dummy => injection_channel
         implicit none
         type (lat_struct), target :: lat
         type (branch_struct), pointer :: branch
         type (ele_struct) ele_at_s, ele_offset
         type (coord_struct), allocatable :: from_orbit(:), to_orbit(:)
         type (coord_struct) orb_at_s, co
         type (em_field_struct) field
         type (g2twiss_struct) twiss, twiss2
         real(rp) circumference
         real(rp) Q(3)
         real(rp) s_end
         real(rp) s_start
         real(rp) s_rel
         real(rp) s, s_eff,c, step/0.01/
         real(rp) s_ele
         real(rp) T(6,6)
         real(rp) temp(6,6), m(6,6)
         logical spin_tracking_on, make_movie
         logical local_ref/.false./,calcd/.false./
         logical err_flag
         logical eflag
         logical local_ref_frame/.false./
         integer j,i, ix_end
         integer ix_ele
         integer lun
         integer nbranch/1/
         integer ie_from, track_state
         integer n
         integer k
  
         i=1
         ! use this routine when lattice twiss are given 10cm upstream of IBMS1. Use lattice file with injection channel that starts 10cm upstream from IBMS.
         !first propagate twiss from lattice file backwards to end of Q025 (-0.909 m upstream from IBMS1 which would be -0.809 m upstream from start of this special lattice
         ! the transfer matrix backwards is
         call copy_lat_to_twiss(lat,twiss,0)
         T(:,:)=0
         forall(i=1:6)T(i,i)=1.
         T(1,2) = -0.809
         T(3,4) = -0.809
         call  prop_phase_space(T, twiss, twiss2)
         ele_at_s%s=-0.809
         ele_at_s%a%beta=twiss2%betax; ele_at_s%b%beta=twiss2%betay
         ele_at_s%a%alpha=twiss2%alphax; ele_at_s%b%alpha=twiss2%alphay
         ele_at_s%x%eta=twiss2%etax; ele_at_s%y%eta=twiss2%etax
         ele_at_s%x%etap=twiss2%etapx; ele_at_s%y%etap=twiss2%etapx
         ele_at_s%a%phi=0; ele_at_s%b%phi=0
         ele_at_s%name='End_Q025'

         call write_injection_line_twiss(ele_at_s)
         
         do j=1,lat%branch(0)%n_ele_track
            call lat_make_mat6 (lat, j,from_orbit, ix_branch=0, err_flag=err_flag)
         end do

!calculation here of transport through line here is with respect to from_orbit.         
              temp = 0.
              forall(i=1:6)temp(i,i) = 1.
              branch => lat%branch(0)
              do i = 1, branch%n_ele_track
                 m = matmul(branch%ele(i)%mat6, temp)
                 temp= m
              end do
              lat%branch(0)%ele(0)%mat6=temp


                call twiss_propagate_all(lat, ix_branch=0)
                print '(a)', ' Transfer matrix of the line '
                n=lat%branch(0)%n_ele_track
                do k=0,n
                   print '(/,a,i10,a)','Element =',k,branch%ele(k)%name
                   print '(4es12.4)',(lat%branch(0)%ele(k)%mat6(i,1:4),i=1,4)
                end do
                
                print '(a)','Injection_channel.  twiss at end of injection channel '
                print '(a12,4a12)','-','betax','alphax','betay','alphay'

                print '(a12,4es12.4)','Start',lat%ele(0)%a%beta,lat%ele(0)%a%alpha,lat%ele(0)%b%beta,lat%ele(0)%b%alpha
                print '(a12,4es12.4)','End v1',lat%ele(n)%a%beta,lat%ele(n)%a%alpha,lat%ele(n)%b%beta,lat%ele(n)%b%alpha
                call  prop_phase_space(lat%branch(0)%ele(0)%mat6, twiss, twiss2)
                print '(a12,4es12.4)','Start v2',twiss%betax,twiss%alphax,twiss%betay,twiss%alphay
                print '(a12,4es12.4)','End v2',twiss2%betax,twiss2%alphax,twiss2%betay,twiss2%alphay
                s_end=0.00001
                step=0.01
                orb_at_s%vec=from_orbit(0)%vec
                ix_end=1
                lun=lunget()
                do while(s_end < lat%ele(lat%branch(0)%n_ele_track)%s + step/2.) !the step/2. is to make sure we overlap the very end
                   if(s_end > lat%ele(lat%branch(0)%n_ele_track)%s) s_end = lat%ele(lat%branch(0)%n_ele_track)%s !track to the very end

                   call twiss_and_track_at_s(lat, s_end, ele_at_s, from_orbit, orb_at_s, ix_branch=0,compute_floor_coords=.true.)          
                   ix_end = element_at_s (lat, s_end, .true.)
                   s_rel = s_end - lat%ele(ix_end-1)%s
                   call write_injection_line_trajectory(s_end, orb_at_s%vec, from_orbit(ix_end)%vec,lat%ele(ix_end)%name, ix_end, ele_at_s%floor%r,orb_at_s%state) !end_orb%vec
                   call write_injection_line_twiss(ele_at_s)

                   s_end = s_end + step
                end do


         return
         end