program analyze_fibers_ibms
use precision_def
use nr
use bmad
implicit none

 integer, parameter:: n_time_bins=700000

 type ibms_struct
  integer fiber(1:16), ntotal
  real(rp) average, rms
 end type

 type (ibms_struct) ibms(3,0:n_time_bins)
 
 type (coord_struct) start_orb
 integer ibmsNum
 integer tbin, tbin_max, lun
 integer lun1
 integer i, n, lines/0/
 integer isign/1/
 integer m
 integer k, ix
 integer fiber

 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(1:16)/-0.0225,-0.0195,-0.0165,-0.0135,-0.0105,-0.0075,-0.0045,-0.0015,0.0015,0.0045,0.0075,0.0105,0.0135,0.0165,0.0195,0.0225/
 real(rp) yfiber(1:16)/-0.0225,-0.0195,-0.0165,-0.0135,-0.0105,-0.0075,-0.0045,-0.0015,0.0015,0.0045,0.0075,0.0105,0.0135,0.0165,0.0195,0.0225/
 real(rp) thickness
 real(rp) x_avg

 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(10)/10*''/
 namelist/input/time_bin,start_turns,end_turns,file

   forall(i=1:16)ibms(:,:)%fiber(i) = 0
   
   ibms(:,:)%ntotal = 0
   ibms(:,:)%average = 0
   ibms(:,:)%rms = 0
   tbin_max = 0
ix=0
k=0
lun=lunget()

 open(unit=lun,file='ibmslist.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  
   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,'ibms')/= 0)cycle
     read(string,'(i5,1x,i5,9es18.8)',end=99)ibmsNum, fiber, start_orb%t, start_orb%s, start_orb%vec, thickness
     if(start_orb%t < 0)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)
         ibms(ibmsNum,tbin)%fiber(fiber) = ibms(ibmsNum,tbin)%fiber(fiber) + 1
         ibms(ibmsNum,tbin)%ntotal = ibms(ibmsNum,tbin)%ntotal +1
         n= ibms(ibmsNum,tbin)%ntotal

         if(ibmsNum == 3)then
             ave_last = ibms(ibmsNum,tbin)%average
             ibms(ibmsNum,tbin)%average = ((n-1)*ibms(ibmsNum,tbin)%average + xfiber(fiber))/float(n) 
             ibms(ibmsNum,tbin)%rms = ((n-1) * ibms(ibmsNum,tbin)%rms + xfiber(fiber)**2 +(n-1)*ave_last**2 - n * ibms(ibmsNum,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:3,1:m), data_ave_cos(1:3,1:m), data_rms_cos(1:3,1:m))
  allocate(data_sin(1:3,1:m), data_ave_sin(1:3,1:m), data_rms_sin(1:3,1:m))
  allocate(datac(1:3,1:m), datac_ave(1:3,1:m), datac_rms(1:3,1:m))
  do ibmsNum=1,3
   do i=1,m
    i_abs = abs(i)
    data_cos(ibmsNum,i) = float(ibms(ibmsNum,i_abs)%ntotal)
    data_ave_cos(ibmsNum,i) = ibms(ibmsNum,i_abs)%average
    data_rms_cos(ibmsNum,i) = ibms(ibmsNum,i_abs)%rms
    data_sin(ibmsNum,i) = float(ibms(ibmsNum,i_abs)%ntotal)
    data_ave_sin(ibmsNum,i) = ibms(ibmsNum,i_abs)%average
    data_rms_sin(ibmsNum,i) = ibms(ibmsNum,i_abs)%rms
    datac(ibmsNum,i) = cmplx(float(ibms(ibmsNum,i_abs)%ntotal),0.)
    datac_ave(ibmsNum,i) = cmplx(ibms(ibmsNum,i_abs)%average,0.)
    datac_rms(ibmsNum,i) = cmplx(ibms(ibmsNum,i_abs)%rms,0.)
   end do
   print *,' tbin_max = ', tbin_max,' log(two) = ', log(two), ' n= ',n,' m = ',m,' size(data_cos(ibmsNum,:) = ', size(data_cos(ibmsNum,:))

  call four1(datac(ibmsNum,1:m),isign)
  call four1(datac_ave(ibmsNum,1:m),isign)
  call four1(datac_rms(ibmsNum,1:m),isign)
 end do
   lun = lunget()
   open(unit = lun, file = 'FiberIbms.dat')
 print *,' Writing FiberIbms.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(ibmsNum,:) = ', size(data_cos(ibmsNum,:))
 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)')'#','Ibms','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)')'#','Ibms','Time bin','Average position','RMS position','Total Hits','Freq bin','FFT cmplx tot hits','FFT cmplx average','FFT complex rms'
  do ibmsNum =1,3
     if(ibmsNum /= 3)cycle
   do tbin = 1, m
     write(lun, '(2x,2i18,2es18.4,i18,1x, es18.4, 6es18.4)')ibmsNum, tbin, ibms(ibmsNum,abs(tbin))%average, ibms(ibmsNum,abs(tbin))%rms, &
                   ibms(ibmsNum,abs(tbin))%ntotal, abs(tbin)/time_bin/m/1.e6, &
                   datac(ibmsNum,tbin), datac_ave(ibmsNum, tbin), datac_rms(ibmsNum,tbin)
   end do
   do tbin =  m+1,tbin_max
     write(lun, '(2x,2i18,2es18.4,i18)')ibmsNum, tbin, ibms(ibmsNum,tbin)%average, ibms(ibmsNum,tbin)%rms, ibms(ibmsNum,tbin)%ntotal
   end do

  end do
 close(lun)
  !-----------------------------------------------------------------------
 open(unit=lun,file = 'FiberIbms_fiber.dat')
 lun1=lunget()
 open(unit=lun1,file='FiberIbms_2D-hist.dat')
  write(lun1,'(2x,4a18)')'ibms','time bin','fiber','hits'
 do k=1,ix
   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(ibmsNum,:) = ', size(data_cos(ibmsNum,:))
 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,2a18,16a10)')'#','Ibms','Time bin','fiber 1','fiber 2','fiber 3','fiber 4','fiber 5','fiber 6','fiber 7','fiber 8','fiber 9','fiber 10','fiber 11','fiber 12','fiber 13','fiber 14','fiber 15','fiber 16','x avg'
  do ibmsNum =1,3
     if(ibmsNum /= 3)cycle
     do tbin = 1, tbin_max
        x_avg = sum(xfiber(:)*ibms(ibmsNum,abs(tbin))%fiber(:))/sum(ibms(ibmsNum,abs(tbin))%fiber(:))
        write(lun, '(2x,2i18,16i10,es12.4)')ibmsNum, tbin, ibms(ibmsNum,abs(tbin))%fiber, x_avg
        do i=1,16
           write(lun1,'(2x,4i18)')ibmsNum, tbin, i, ibms(ibmsNum,abs(tbin))%fiber(i)
        end do
        
   end do
  end do
  close(lun)
  close( lun1)
  print *, 'write FiberIbms_fiber.dat'
  print *,' write FiberIbms_2D-hist.dat'
  
 
 end program