program spin_energy
use bmad

real(rp) Q(3)
type (coord_struct) :: coord,ecoord ! for muon and electron
integer n, ios/0/, Nbin(0:1000)/1001*0/, Npe(0:1000)/1001*0/,Nebin(0:100)/101*0/
integer bin, ebin, bin_max/0./,tot/0./
real(rp) turn
real(rp) pe
do while(ios >= 0)
!   read(41,'(i10,9es12.4,6es12.4,1x,3es12.4)',IOSTAT=ios) n,turn,coord%t,coord%path_len,ecoord%vec,coord%vec, Q
   read(41,*,IOSTAT=ios) n,turn,coord%t,coord%path_len,ecoord%vec,coord%vec, Q
pe= (1+ecoord%vec(6))*3.095
bin = (coord%t*1e6+0.5)
tot=tot+1
Nbin(bin) = Nbin(bin)+1
ebin = pe/3.095*100 +0.5
Nebin(ebin) = Nebin(ebin)+1
print '(i10,es12.4,3i10)',n,turn,bin,Nbin(bin),tot
if(pe>0.7*3.095)Npe(bin)=Npe(bin)+1
bin_max = max(bin,bin_max)
end do

do n=0,bin_max
   write(31,'(3i10)')n,Nbin(n),Npe(n)
end do
do n=0,100
   write(32,'(2i10)')n,Nebin(n)
end do
print *,'total electrons = ', tot
end program spin_energy