subroutine first_turn(lat,spin_tracking_on, from_orbit, circumference, make_movie, start_tracking_at_inflector_exit)
         use bmad
         use parameters_bmad
         use muon_interface, dummy => first_turn
         implicit none
         type (lat_struct) lat
         type (ele_struct) ele_at_s, ele_offset, ele
         type (ele_struct), pointer :: lord
         type (coord_struct), allocatable :: from_orbit(:), to_orbit(:)
         type (coord_struct) orb_at_s, co
         type (em_field_struct) field
         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.1/
         real(rp) s_ele
         logical spin_tracking_on, make_movie
         logical local_ref/.false./,calcd/.false./
         logical err_flag
         logical eflag
         logical local_ref_frame/.false./
         logical start_tracking_at_inflector_exit
         integer j,i, ix_end
         integer ix_ele
         integer lun
         integer nbranch/1/
         integer ie_from, track_state
  
         i=1
         lun=lunget()
         open(unit=lun, file=trim(directory)//'/'//'inj_traj_bfield.dat')
        write(lun,'(12a12)')'s', 'x','xp','y','yp','z','dp','Bx','By','Bz','state','t'

        s_end=0.00001
        step=0.1

        if(.not. start_tracking_at_inflector_exit)then
            do j=1,lat%branch(0)%n_ele_track
               if(spin_tracking_on) Q = from_orbit(j)%spin
               call write_single_particle_by_element(i,lat%branch(0)%ele(j)%name,from_orbit(j), lat%branch(0)%ele(j)%s/circumference, &
                    ((i-1)*lat%branch(0)%ele(lat%branch(0)%n_ele_track)%s +lat%branch(0)%ele(j)%s)/circumference, Q)
               call lat_make_mat6 (lat, j,from_orbit, ix_branch=0, err_flag=err_flag)
            end do

            call twiss_propagate_all(lat, ix_branch=0)
            do j=1,lat%branch(0)%n_ele_track
               s_start=0.
               if(make_movie)call step_floor_around_ring(lat,j,0,from_orbit(j-1),lat%branch(0)%ele(j-1)%floor%r, s_start)
            end do
 
            orb_at_s%vec=from_orbit(0)%vec
            ix_end=1
            
            write(lun,'(a,es12.4)')' inflector_field = ', inflector_field
            call write_injection_line_trajectory(s_end, orb_at_s%vec, from_orbit(1)%vec,lat%ele(1)%name, ix_end,lat%ele(1)%floor%r,orb_at_s%state)
            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.)          
!                   print '(es12.4,1x,i10,1x,a,1x,es12.4,1x,a)', s_end, lat%branch(0)%n_ele_track,lat%ele(lat%branch(0)%n_ele_track)%name, lat%ele(lat%branch(0)%n_ele_track)%s, ele_at_s%name
               ix_end = element_at_s (lat, s_end, .true.)
               s_rel = s_end - lat%ele(ix_end-1)%s
               call em_field_calc(lat%ele(ix_end), lat%param, s_rel, orb_at_s, local_ref, field,calcd, err_flag)
               write(lun,'(es12.4,1x,6es12.4,3es12.4,1x,i12, es12.4)')s_end, orb_at_s%vec, field%B(1:3), orb_at_s%state, orb_at_s%t 
               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
            s_end = s_end - step
         endif

! add this next bit 6/23/20 to propagate twiss through injection line and into ring
         nbranch = 1
         if(ele_at_s%a%beta /= 0 .and. ele_at_s%b%beta /=0) then
            lat%branch(nbranch)%ele(0)%a%beta = ele_at_s%a%beta
            lat% branch(nbranch)%ele(0)%a%alpha = ele_at_s%a%alpha
            lat%branch(nbranch)%ele(0)%b%beta = ele_at_s%b%beta
            lat% branch(nbranch)%ele(0)%b%alpha = ele_at_s%b%alpha
            lat%branch(nbranch)%ele(0)%x%eta = ele_at_s%x%eta
            lat%branch(nbranch)%ele(0)%x%etap = ele_at_s%x%etap
         endif
         ie_from = lat%branch(nbranch)%ix_from_ele
         
         call twiss_propagate_all(lat,ix_branch = nbranch)
         !                from_orbit(0) = orb_at_s
         call reallocate_coord(to_orbit, lat%branch(nbranch)%n_ele_max)
         to_orbit(0) = from_orbit(ie_from)
         call track_all(lat, to_orbit, nbranch, track_state, err_flag)
         s=0.
         do while(s <= lat%branch(nbranch)%ele(lat%branch(nbranch)%n_ele_track)%s)
            ix_ele= element_at_s(lat,s,choose_max=.true., ix_branch=nbranch, err_flag=eflag, s_eff=s_eff, position=co)
            call twiss_and_track_at_s(lat,s,ele_at_s, to_orbit,  orb_at_s=orb_at_s, ix_branch=nbranch, use_last=.false.)
            s_ele = lat%branch(nbranch)%ele(ix_ele)%s_start
            ele = lat%branch(nbranch)%ele(ix_ele)
            if(index(ele%name,'DFREE')/=0)then
               lord => pointer_to_lord(ele,1)
            endif
            if(index(ele%name,'KICKER')/=0)then
               print '(a)', ix_ele, ele%name
            endif
                   
                   
            call em_field_calc(lat%branch(nbranch)%ele(ix_ele), lat%param,s-s_ele,orb_at_s, local_ref_frame, field, err_flag=err_flag)
            write(lun,'(es12.4,1x,6es12.4,3es12.4,1x,i12, es12.4)')s_end+s,  orb_at_s%vec, field%B(1:3), orb_at_s%state, orb_at_s%t
            ele_offset = ele_at_s
            ele_offset%s =ele_at_s%s + s_end
            call write_injection_line_twiss(ele_offset)           
!                   write(lun,'(i10,16es12.4,1x,a16,1x,3es12.4)')ix_ele,s+(i-1)*c,s_eff, ele_at_s%a%beta, ele_at_s%a%alpha, ele_at_s%b%beta, ele_at_s%b%alpha, ele_at_s%x%eta,&
!                                          ele_at_s%x%etap, ele_at_s%y%eta, ele_at_s%y%etap, lat%branch(nbranch)%ele(ix_ele)%a%beta,&
!                                          lat%branch(nbranch)%ele(ix_ele)%b%beta,orb_at_s%vec(1),orb_at_s%vec(3), &
!                                          closed_orb(ix_ele)%vec(1), closed_orb(ix_ele)%vec(3),lat%branch(nbranch)%ele(ix_ele)%name, field%b
            s= s + step
         end do


         close(lun)
         
         return
       end subroutine first_turn