program parse_fort41

use bmad
implicit none
real(rp) turn
real(rp),dimension(3):: Q                        !polarization 
real(rp) range_min/0./, range_max/0.2/
real(rp) frac_turn
integer  n             !which muon
integer j
integer num_files
integer ios
integer nevents
type (coord_struct) :: coord,ecoord ! for muon and electron
type (pauli_struct),dimension(3):: sig           !pauli matrices
character*100 string

character*72 dir/'/nfs/gm2/data1/dselezne/JobOutput/'/, name(300)/300*'0'/

namelist/files/ name, range_min, range_max

  OPEN (UNIT=5, FILE='file_name.dat', STATUS='old', IOSTAT=ios)
    READ (5, NML=files, IOSTAT=ios)
    WRITE(6,NML=files)
  print *, 'ios=', ios
   rewind(unit=5)
!    READ (5, NML=files)
  CLOSE(5)



j=1
do while(name(j) /= '0')
nevents=0
open(unit=41, file=trim(dir)//trim(name(j))//'/fort.41')
 print '(a,a)',' read file =',trim(dir)//trim(name(j))//'/fort.41'
 read(41,'(a)')string
 do while(.true.)
  read(41,'(i10,9es12.4,6es12.4,1x,3es12.4)', end=99) n,turn,coord%t,coord%path_len,ecoord%vec,coord%vec, Q
  frac_turn = turn - int(turn)
  if(range_min < frac_turn .and. frac_turn < range_max)then
   write(42, '(i10,9es12.4,6es12.4,1x,3es12.4,es14.6)') n,turn,coord%t,coord%path_len,ecoord%vec,coord%vec, Q, frac_turn
   nevents=nevents+1
  endif

 end do
99 continue
 print '(i10,a,a)',nevents,' events from file =',trim(dir)//trim(name(j))//'/fort.41'
j=j+1
end do

end