program combine_Beam_moments_tbt
use precision_def
use muon_mod
implicit none

type (g2moment_struct) moment

 integer lun
integer number_remaining(0:10000)
integer turn, nmu,ix
integer ios
integer k, max_turn

character*300 string
character*200 file(100)/100*''/
character*16 ele_name
real(rp) moment_ave_sqr, turn_x

namelist/input/file

max_turn=0
number_remaining(:)=0


lun=lunget()

 open(unit=lun,file='dir_list.dat', STATUS='old', IOSTAT=ios)
    READ (lun, NML=input, IOSTAT=ios)
    WRITE(6,NML=input)
  print *, 'ios=', ios
   rewind(unit=lun)
    READ (lun, NML=input)
  CLOSE(lun)

do k=1,size(file)
 if(trim(file(k)) /= '')print *,' file = ',trim(file(k))
end do
! write(lun(unit),'(a7,a13,12a12,a13,a10,1x,a16)') &
!                                            'Element','Turn','x_ave','y_ave','t_ave', &
!                                            'sigma(1,1)','sigma(1,2)','sigma(2,2)', &
!                                            'sigma(3,3)','sigma(3,4)','sigma(4,4)', &
!                                            'sigma(5,5)','sigma(5,6)','sigma(6,6)', &
!                                            'dp/p sqr_ave','remaining','element name    '

 
 do k=1,size(file)
   if(trim(file(k)) == '')cycle  
   lun=lunget()
   open(unit=lun, file = file(k))
   print *,' file = ',trim(file(k))
   do while(.true.)

 read(lun,'(a)', end=99)string
  if(index(string,'Turn')/=0 .or. index(string,'harp')/=0)cycle
 print '(a)', string
 read(string, '(i7,es13.5,12es12.4,es13.4,i10,a)') &
                              ix, turn_x, moment%ave(1), moment%ave(3), moment%ave(5), &
                              moment%sigma(1,1), moment%sigma(1,2), moment%sigma(2,2), &
                              moment%sigma(3,3), moment%sigma(3,4), moment%sigma(4,4), &
                              moment%sigma(5,5), moment%sigma(5,6), moment%sigma(6,6), &
                              moment_ave_sqr, nmu, ele_name


  turn=turn_x+.5
  number_remaining(turn) = nmu + number_remaining(turn)
  max_turn = max(turn, max_turn)
print '(4i12)',turn,number_remaining(turn), max_turn, nmu
   end do

 99 continue
  close(lun)

end do ! files
print *,' max_turn = ', max_turn
print *,' write file: muons_remaining.dat'

 lun = lunget()
 open(unit=lun, file = 'muons_remaining.dat')
write(lun,'(4a19)')'turn','number remaining',' fraction remaining' , 'fraction lost after turn'
do turn = 0,max_turn
 write(lun,'(2i19,2es19.6)')turn, number_remaining(turn), float(number_remaining(turn))/number_remaining(0), &
         float(number_remaining(turn)-number_remaining(max_turn))/number_remaining(max_turn)
end do
close(lun)
end