Subroutine  inj_coord_eff( ring, inj, synch_beam, closed_orbit, &
		coord)

  use bmad_interface
  use injlib

  implicit none

  type (lat_struct)        ring
  type (coord_struct), allocatable ::  closed_orbit(:)
  type (coord_struct), allocatable ::  inj_orbit(:)
  type (inj_struct)         inj
  type (synch_beam_struct)  synch_beam

  integer  i,j,k, index, ix_start, n, style
  integer  lun1, lun2
  integer  n_part, nturns
  integer  coord, track_state
  real  nsig
  real  lost_array(0:ring%n_ele_track), array(0:1000)
  real  totweight, weight, injweight
  real  min, max, step
  real  injval, eff
  real  injcoord(2)
  logical  lostit, random


! allocate storage
  call reallocate_coord (inj_orbit, ring%n_ele_max)

! Set these along with angmin,angmax,angstep in data file later
  nsig   = 2.5
  n_part = 500
  nturns = 50
  random = .true.
  style  = 2
  
  index  = 0

  bmad_com%aperture_limit_on = .true.

!  ***Do I need this??***
!  call twiss_and_track(ring, closed_orbit_)
  
  call inject_particle(ring, .false., 0, 0., weight, inj, &
			synch_beam,inj_orbit)
  if (coord==1) then
     print *,'Injection Position = ',inj_orbit(inj%ix_injpt)%vec(1)
     print *,'Enter Injection Position Min, Max, Step'
     read  *, min, max, step
  else if (coord==2) then
     print *,'Injection Angle = ',inj_orbit(inj%ix_injpt)%vec(2)
     print *,'Enter Injection Angle Min, Max, Step'
     read  *, min, max, step
  else
     print *,coord,'is an invalid coordinate or not implemented yet'
     return
  endif

  do injval = min, max, step

     call inject_particle(ring, .false., 0, 0., weight, inj, &
    			  synch_beam, inj_orbit)

     inj_orbit(inj%ix_injpt)%vec(coord) = injval

     !Take Beam to Injection Point
     ix_start = inj%ix_injpt
     n        = 0
     lostit   = .false.

     !***This can be done in track_injected_particle***
     !Inflate Aperture for Q34E And S34E For Zeroth Turn
     ring%ele(inj%ix_injpt)%value(x1_limit$)=0.10
     ring%ele(inj%ix_injsex)%value(x1_limit$)=0.10
     ring%ele(inj%ix_injpt)%value(x2_limit$)=0.10
     ring%ele(inj%ix_injsex)%value(x2_limit$)=0.10

     do while(n<nturns .and. .not. lostit )

 	call pulsed_bump(ring, n, inj)
  	call track_many(ring, inj_orbit, ix_start, 0, -1, track_state = track_state)

	ix_start = 0
  	
	n = n + 1

	!Return Aperture to original size
	ring%ele(inj%ix_injpt)%value(x1_limit$)=0.045
	ring%ele(inj%ix_injsex)%value(x1_limit$)=0.045
	ring%ele(inj%ix_injpt)%value(x2_limit$)=0.045
	ring%ele(inj%ix_injsex)%value(x2_limit$)=0.045

	if(track_state /= moving_forward$ ) lostit=.true.

     enddo !n<n_turns and lostit==.false.

     if( lostit ) then
        if (coord==1) then 
	    print 100,'Lost on Turn ',n,' for position = ',injval
	else if (coord==2) then
	    print 100,'Lost on Turn ',n,' for angle = ',injval
        endif
 100	format(1x,a,1x,i4,1x,a,1x,e12.5)
     else
	if (coord==1) then
	   print 110,'Beam survived for injection position = ',injval
	else if (coord==2) then
           print 110,'Beam survived for injection angle = ',injval
	endif
 110	format(1x,a,e12.5) 
	array(index)=injval
        index=index+1
     endif

    enddo ! injval

   do i=0,ring%n_ele_track
      lost_array(i)=0.
   enddo

   lun1=lunget()
   if (coord==1) then
      open(lun1,file='inj_xpos_eff.dat',status='new')
   else if (coord==2) then
      open(lun1,file='inj_xprime_eff.dat',status='new')
   endif
!   lun2=lunget()
!   open(lun2,file='lost_array.dat',status='new')
   
   nsig=abs(nsig)

   do j=0,(index-1)

      injval=array(j)
      
      totweight = 0.0
      injweight = 0.0

      do k=1,n_part
	
	call inject_particle_at_coord(ring, style, nsig, weight, &
		inj, synch_beam, inj_orbit, injval, coord)

!	print *,'Injval - ',injval
! 	print *,'inj_orbit(inj%ix_injpt)%vec(coord) - ', &
!		inj_orbit(inj%ix_injpt)%vec(coord)

	totweight = totweight + weight

	call track_injected_particle(ring, inj, electron$, inj_orbit, &
		nturns, .false., 1, 1, i, weight, lostit, track_state)

	if( lostit ) then
	  lost_array(track_state) = lost_array(track_state) + weight
	else
	  injweight = injweight + weight
	endif

	if (mod(k,100)==0) print *,'Finished tracking particle ',k
	
      enddo ! k=1,npart

      eff = injweight/totweight

      write(lun1,200) injval, eff
 200  format(1x,f10.7,1x,f7.4)
      if (coord==1) then
         print 210,'Position - ',injval,' Efficiency - ', eff
      else if (coord==2) then
         print 210,'Angle - ',injval,' Efficiency - ', eff
      endif
 210  format(1x,a,e10.3,1x,a,f9.6)
 
!      do i = 0,ring%n_ele_track
!	lost_array(i)=lost_array(i)/totweight
!	write(lun2,220)i,ring%ele(i)%s,lost_array(i)
! 220	format(1x,i5,1x,f10.4,1x,f10.6)
!      enddo !i

    enddo !j=1,index    

   close(lun1)
!   close(lun2)

   return

   end

subroutine  inj_pinger_eff( ring, inj, synch_beam )

  use bmad_interface
  use injlib

  implicit none

  type (lat_struct)        ring
  type (coord_struct), allocatable ::  inj_orbit(:)
  type (inj_struct)         inj
  type (synch_beam_struct)  synch_beam
  type (pretzel_struct)     pretzel
  type (qtune_struct)       qtune
  type (xqune_struct)       xqune
  type (tunes_struct)       tunes

  integer  i,j,k, index, ix_start, n, style
  integer  lun1, lun2
  integer  n_part, nturns, track_state, track_state2
  integer  coord
  real  nsig
  real  lost_array(0:ring%n_ele_track)
  real  totweight, weight, injweight
  real  min, max, step
  real  pingval, kping, eff
  logical  lostit, lostit2, random
  logical  wholehog

! allocate storage
  call reallocate_coord (inj_orbit, ring%n_ele_max)
 
  wholehog=.true.
! Set these along with angmin,angmax,angstep in data file later
  nsig   = 2.5
  n_part = 500
  nturns = 50
  random = .true.
  style  = 2
  
  index  = 0

  bmad_com%aperture_limit_on = .true.

  print *,'Enter kping min, max, step:'
  read *,min,max,step
!  print *,' min,max,step=',min,max,step

!  lun1=lunget()
!  open(lun1,file='inj_pinger_eff.dat',status='unknown')

  do pingval=min,max,step

     print *,' pingval is',pingval
     
     kping=pingval/1000.

     inj%pingkick=kping

     lostit = .false.
     totweight=0.0
     injweight=0.0

     do i = 0,ring%n_ele_track
	inj_orbit(i)%vec(:) = 0.
     enddo

     call inject_particle(ring, .false., 0, 0., weight, inj, &
		synch_beam, inj_orbit)

     call track_injected_particle(ring, inj, electron$, inj_orbit, nturns, &
		.false.,1,1,1,weight,lostit,track_state)

	print *,' lostit is ',lostit

!!     if (lostit) then
	print *,'Lost at ping value: ',pingval
!!     else if (lostit==.false. .and. wholehog==.true.) then
        print *,pingval,' survived ideal 50 loops.'
	print *,'Now calculating efficiency with ',n_part,' particles.'

	totweight=0.0
	injweight=0.0
	
	do j=1,n_part

	   call inject_particle(ring, random, style, nsig, weight, &
		inj, synch_beam, inj_orbit)
!		weight = 1
	
	   totweight=totweight+weight
	   
	   call track_injected_particle(ring, inj, electron$, inj_orbit, &
		nturns, .false., 1, 1, j, weight, lostit2, track_state2)
!		lostit2 = .false.

	   if (lostit2) then 
	      lost_array(track_state2) = lost_array(track_state2) + weight
	   else
	      injweight=injweight+weight
	   endif

	   if (mod(j,100)==0) print *,'Finished tracking particle ',j

	enddo !j=1,npart

	eff = injweight/totweight

!	write(lun1,300)pingval,eff
! 300	format(1x,f8.4,1x,f7.4)
	print 310,'Pinger Kick Value - ',pingval,' Efficiency - ',eff
 310	format(1x,a,f8.4,1x,a,f7.4)

!!     else
!	print *,pingval, 'survived 50 ideal turns'
!!     endif !lostit

  enddo !pingval
 
!  close(lun1)

  return

  end	

subroutine  inject_particle_at_coord( ring, style, nsig, weight, inj, &
		synch_beam, orbit, injval, coord)

  use bmad_interface
  use injlib

  implicit none

  type (lat_struct)  ring
  type (inj_struct)  inj
  type (synch_beam_struct)  synch_beam
  type (coord_struct), allocatable ::  orbit(:)
  type (normal_modes_struct)  mode

  integer  i, style
  real  nsig, weight
  real  vecnorm(6), vecreal(6), sigma(3)

  real injval
  integer coord

  sigma(1) =  synch_beam%a%sigma  ! H beamsize
  sigma(2) =  synch_beam%b%sigma  ! V beamsize
  sigma(3) =  synch_beam%length   ! bunch length

  call random_phase_space( style, vecnorm, sigma, nsig, weight )

!  convert from normalized coordinates to real 6-d coordinates
  vecreal(1) = vecnorm(1)   ! X
  vecreal(2) = (vecnorm(2) - synch_beam%a%alpha*vecnorm(1))/ & ! X-prime  
               synch_beam%a%beta
  vecreal(3) = vecnorm(3)   ! Y
  vecreal(4) = (vecnorm(4) - synch_beam%b%alpha*vecnorm(3))/ & ! Y-prime  
               synch_beam%b%beta
  vecreal(5) = vecnorm(5)   ! Phase Error (z position)
  vecreal(6) = synch_beam%dele/synch_beam%length * vecnorm(6)

!  Now offset phase space centroid according to energy error
  vecreal(6) = synch_beam%centroid(6) + vecreal(6)  ! Possible energy err
  vecreal(5) = synch_beam%centroid(5) + vecreal(5)  ! Possible phase err

   !Modify position of injection - default was synch_beam%centroid
  if (coord==1) then
     vecreal(1) = injval + vecreal(1) +  &
                  vecreal(6)*synch_beam%a%eta 
     vecreal(2) = synch_beam%centroid(2) + vecreal(2) +  &
                  vecreal(6)*synch_beam%a%etap 
  else if (coord==2) then
     vecreal(1) = synch_beam%centroid(1) + vecreal(1) +  &
                  vecreal(6)*synch_beam%a%eta
     vecreal(2) = injval + vecreal(2) +  &
                  vecreal(6)*synch_beam%a%etap
  end if
   
  vecreal(3) = synch_beam%centroid(3) + vecreal(3) +  &
               vecreal(6)*synch_beam%b%eta
  vecreal(4) = synch_beam%centroid(4) + vecreal(4) +  &
               vecreal(6)*synch_beam%b%etap

  do i = 1,6
    orbit(inj%ix_injpt)%vec(i) = vecreal(i)
  enddo


  end