program analyze_harp_plane_hits
use precision_def
use nr
use bmad
implicit none

 integer, parameter:: n_time_bins=600000

 type (coord_struct) start_orb

 type fiber_plane_struct
   real(rp) avg(1:6), rms(1:6)
   integer ntotal
end type

type (fiber_plane_struct) harp(1:4,n_time_bins)
   
 integer fiberNum
 integer tbin, tbin_max/0/, lun, lun1
 integer ios
 integer i, n, lines/0/
 integer j

 real(rp) time_bin/1.e-8/
 real(rp) ave_last(1:6)
 integer k
 integer num_fib/0/

 character*160 string, file(100)/100*''/
 character*6 wbin
 character*100 out_name
 logical active_fib(-3:3)/7*.false./

namelist/input/time_bin,file

lun=lunget()

 open(unit=lun,file='hpd_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

do i=1,n_time_bins
 do j=1,6
  harp(1:4,i)%avg(j)=0
  harp(1:4,i)%rms(j)=0
 end do
 harp(1:4,i)%ntotal = 0
end do

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
!     write(lun1,'(a)')string
     if(index(string,'harp')/= 0)cycle
     read(string,'(i5,8es18.8)',end=99)fiberNum, start_orb%t, start_orb%s, start_orb%vec
     lines = lines + 1
     if((lines/1000000)*1000000 == lines)print *, lines

       tbin = (start_orb%t / time_bin)+.5
     
        if(tbin > n_time_bins)then
          print *,' tbin = ',tbin ,'> n_time_bins = ', n_time_bins
          stop
        endif

        tbin_max = max(tbin, tbin_max)

        if(tbin <= 0) cycle
        harp(fiberNum,tbin)%ntotal = harp(fiberNum,tbin)%ntotal +1
        n= harp(fiberNum,tbin)%ntotal

        ave_last(1:6) = harp(fiberNum,tbin)%avg(1:6)
        harp(fiberNum,tbin)%avg(1:6) = ((n-1)*harp(fiberNum,tbin)%avg(1:6) + start_orb%vec(1:6))/n
        harp(fiberNum,tbin)%rms(1:6) = ((n-1) * harp(fiberNum,tbin)%rms(1:6) + &
            start_orb%vec(1:6)**2 +(n-1)*ave_last(1:6)**2-  n* harp(fiberNum,tbin)%avg(1:6)**2   )/n
   end do

 99 continue
  close(lun)
print *, 'tbin_max = ' ,tbin_max
end do ! files

 lun = lunget()
 open(unit=lun, file = 'harp_plane_hit_anal.dat')


do fiberNum=1,4
 write(lun,'(a1)')'#'
 write(lun,'(16a12)')'harp','time bin','time','total','x avg','xp avg','y avg','yp avg','z avg','zp avg','xx','xpxp','yy','ypyp','zz','zpzp'
 write(lun,'(a1)')'#'
 do tbin=1,tbin_max
   write(lun,'(2i12,es12.4,i12,12es12.4)') fiberNum,tbin,tbin*time_bin,harp(fiberNum,tbin)%ntotal,harp(fiberNum,tbin)%avg(1:6), harp(fiberNum,tbin)%rms(1:6)
 end do
end do
close(lun)
end