program analyze_fibers_harp
use precision_def
use nr
use bmad
implicit none

 integer, parameter:: n_time_bins=700000

 type harp_struct
  integer fiber(-3:3), ntotal
  real(rp) average, rms
 end type

 type (harp_struct) harp(4,0:n_time_bins)
 
 type (coord_struct) start_orb
 integer harpNum
 integer tbin, tbin_max, lun
 integer i, n, lines/0/
 integer isign/1/
 integer m
 integer k, ix
 integer fiber
 integer number_files

 real(rp) two/2./
 real(rp), allocatable :: data_cos(:,:), data_ave_cos(:,:), data_rms_cos(:,:)
 real(rp), allocatable :: data_sin(:,:), data_ave_sin(:,:), data_rms_sin(:,:)
 real(rp) time_bin/1.e-8/, radius/0.00025/
 real(rp) ave_last
 real(rp) xfiber(-3:3)/-0.039,-0.026,-0.013,0.,0.013,0.026,0.039/
 real(rp) yfiber(-3:3)/-0.039,-0.026,-0.013,0.,0.013,0.026,0.039/
 real(rp) thickness

 complex(rp), allocatable :: datac(:,:), datac_ave(:,:), datac_rms(:,:)
 integer ix_hit, iy_hit
 integer ix_timebin, ix_radius, ix_turns
 integer i_abs
 integer start_turns/0/, end_turns/0/
 integer ios

 character*120 string, file(100)/100*''/
 character*15 dir
 logical active_fibers(-3:3)/7*.false./
 namelist/input/time_bin,start_turns,end_turns,active_fibers,file

 ! initialize harp positions


!   open(unit=lun, file = 'harp_plane_hits.dat')
   forall(i=-3:3)harp(:,:)%fiber(i) = 0
   
   harp(:,:)%ntotal = 0
   harp(:,:)%average = 0
   harp(:,:)%rms = 0
   tbin_max = 0
ix=0
k=0
lun=lunget()

 open(unit=lun,file='flist.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 k=1,size(file)
   if(trim(file(k)) == '')cycle
   number_files = k   
   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,'harp')/= 0)cycle
     read(string,'(i5,1x,i5,9es18.8)',end=99)harpNum, fiber, start_orb%t, start_orb%s, start_orb%vec, thickness
!     print '(i5,1x,i5,9es18.8)', harpNum, fiber, start_orb%t, start_orb%s, start_orb%vec, thickness
     fiber = fiber - 4 !1:7 to -3:3
     if(.not. active_fibers(fiber))cycle
!     if((start_orb%t < start_turns*149.e-9 .or. start_orb%t > end_turns * 149.e-9) .and. end_turns /= 0)cycle
     if(start_orb%t < start_turns*149.e-9)cycle
     lines = lines + 1
     if((lines/1000000)*1000000 == lines)print *, lines
         tbin = (start_orb%t - start_turns*149.e-9)/time_bin+.5
         tbin_max = max(tbin, tbin_max)
         harp(harpNum,tbin)%fiber(fiber) = harp(harpNum,tbin)%fiber(fiber) + 1
         harp(harpNum,tbin)%ntotal = harp(harpNum,tbin)%ntotal +1
         n= harp(harpNum,tbin)%ntotal

         if(harpNum == 1 .or. harpNum == 3)then
             ave_last = harp(harpNum,tbin)%average
             harp(harpNum,tbin)%average = ((n-1)*harp(harpNum,tbin)%average + xfiber(fiber))/float(n) 
             harp(harpNum,tbin)%rms = ((n-1) * harp(harpNum,tbin)%rms + xfiber(fiber)**2 +(n-1)*ave_last**2 - n * harp(harpNum,tbin)%average**2 )/float(n)
          elseif(harpNum == 2 .or. harpNum == 4)then
             ave_last = harp(harpNum,tbin)%average
             harp(harpNum,tbin)%average = ((n-1)*harp(harpNum,tbin)%average + yfiber(fiber))/float(n) 
             harp(harpNum,tbin)%rms = ((n-1) * harp(harpNum,tbin)%rms + yfiber(fiber)**2 +(n-1)*ave_last**2 - n * harp(harpNum,tbin)%average**2)/float(n)
         endif

    end do
 99 continue
   close(lun)
end do ! files

! fourier transform of Number vs time
  n = log(float(tbin_max))/log(two)
  m = two**n 

  allocate(data_cos(1:4,1:m), data_ave_cos(1:4,1:m), data_rms_cos(1:4,1:m))
  allocate(data_sin(1:4,1:m), data_ave_sin(1:4,1:m), data_rms_sin(1:4,1:m))
  allocate(datac(1:4,1:m), datac_ave(1:4,1:m), datac_rms(1:4,1:m))
  do harpNum=1,4
   do i=1,m
    i_abs = abs(i)
    data_cos(harpNum,i) = float(harp(harpNum,i_abs)%ntotal)
    data_ave_cos(harpNum,i) = harp(harpNum,i_abs)%average
    data_rms_cos(harpNum,i) = harp(harpNum,i_abs)%rms
    data_sin(harpNum,i) = float(harp(harpNum,i_abs)%ntotal)
    data_ave_sin(harpNum,i) = harp(harpNum,i_abs)%average
    data_rms_sin(harpNum,i) = harp(harpNum,i_abs)%rms
    datac(harpNum,i) = cmplx(float(harp(harpNum,i_abs)%ntotal),0.)
    datac_ave(harpNum,i) = cmplx(harp(harpNum,i_abs)%average,0.)
    datac_rms(harpNum,i) = cmplx(harp(harpNum,i_abs)%rms,0.)
   end do
   print *,' tbin_max = ', tbin_max,' log(two) = ', log(two), ' n= ',n,' m = ',m,' size(data_cos(harpNum,:) = ', size(data_cos(harpNum,:))

!  call cosft2(data_cos(harpNum,1:m),isign)
!  call cosft2(data_ave_cos(harpNum,1:m),isign)
!  call cosft2(data_rms_cos(harpNum,1:m),isign)
!  call sinft(data_sin(harpNum,1:m))
!  call sinft(data_ave_sin(harpNum,1:m))
!  call sinft(data_rms_sin(harpNum,1:m))
  call four1(datac(harpNum,1:m),isign)
  call four1(datac_ave(harpNum,1:m),isign)
  call four1(datac_rms(harpNum,1:m),isign)
 end do
   lun = lunget()
   open(unit = lun, file = 'FiberHarp.dat')
 print *,' Writing FiberHarp.dat'
 write(lun,'(a)')'#! Data from files '

  do k=1,ix
   print *,' k = ',k,' file(k) = ', trim(file(k))
   write(lun,'(a2,1x,a)')'#!', trim(file(k))
  end do
  write(lun, *)' tbin_max = ', tbin_max,' log(two) = ', log(two), ' n= ',n,' m = ',m,' size(data_cos(harpNum,:) = ', size(data_cos(harpNum,:))
 write(lun,'(a1,1x,a,es12.4)')'#',' Time bin = ',time_bin
 write(lun,'(a1,1x,a,es12.4)')'#',' fiber_radius = ',radius
 write(lun,'(a1,1x,a,1x,i10,a,1x,i10)')'#','Start turns = ', start_turns, 'End turns = ',end_turns
!  write(lun,'(a1,1x,8a18,3a18,3a36)')'#','Harp','Time bin','Average position','RMS position','Total Hits','FFT cos tot hits','FFT cos average','FFT cos rms', &
!                                                                                  'FFT sin tot hits','FFT sin average','FFT sin rms','FFT cmplx tot hits','FFT cmplx average','FFT complex rms'
  write(lun,'(a1,1x,6a18,3a36)')'#','Harp','Time bin','Average position','RMS position','Total Hits','Freq bin','FFT cmplx tot hits','FFT cmplx average','FFT complex rms'
  do harpNum =1,4
   do tbin = 1, m
!    if(harp(harpNum,tbin)%ntotal /= 0)write(lun, '(2i12,2es12.4,i10)')harpNum, tbin, harp(harpNum,tbin)%average, harp(harpNum,tbin)%rms, harp(harpNum,tbin)%ntotal
!     write(lun, '(2x,2i18,2es18.4,i18,1x,6es18.4, 6es18.4)')harpNum, abs(tbin), harp(harpNum,abs(tbin))%average, harp(harpNum,abs(tbin))%rms, &
!                   harp(harpNum,abs(tbin))%ntotal, &
!                   data_cos(harpNum,tbin), data_ave_cos(harpNum, tbin), data_rms_cos(harpNum,tbin), &
!                   data_sin(harpNum,tbin), data_ave_sin(harpNum, tbin), data_rms_sin(harpNum,tbin), &
!                   datac(harpNum,tbin), datac_ave(harpNum, tbin), datac_rms(harpNum,tbin)
     write(lun, '(2x,2i18,2es18.4,i18,1x, es18.4, 6es18.4)')harpNum, tbin, harp(harpNum,abs(tbin))%average, harp(harpNum,abs(tbin))%rms, &
                   harp(harpNum,abs(tbin))%ntotal, abs(tbin)/time_bin/m/1.e6, &
                   datac(harpNum,tbin), datac_ave(harpNum, tbin), datac_rms(harpNum,tbin)
   end do
   do tbin =  m+1,tbin_max
     write(lun, '(2x,2i18,2es18.4,i18)')harpNum, tbin, harp(harpNum,tbin)%average, harp(harpNum,tbin)%rms, harp(harpNum,tbin)%ntotal
   end do

  end do
  close(unit=lun)
  if(number_files == 1) then
     dir=file(1)
!     dir= dir(1:index(dir,'/'))
     string = 'cp FiberHarp.dat '//trim(dir)//'/FiberHarp.dat'
     print *, string
     call execute_command_line(string)
     print *,' write FiberHarp.dat to '//trim(dir)
  endif
 
 end program