program analyze_Time_dep

use precision_def
use nr
use bmad
use naff_mod
!use muon_mod
implicit none

type g2moment_struct
 real(rp) sigma(6,6), ave(6), third(6)
 real(rp) max_sigma(6,6), min_sigma(6,6), max_ave(6), min_ave(6)
end type

TYPE (g2moment_struct), allocatable :: moments(:) !moment structure for each time bin, 1000 turns => an array with 1000*149e-9/bin_width 

type running_moment_struct
  real(rp) product(1:6,1:6), sum(1:6), third(6)
  integer n
end type

type (running_moment_struct), allocatable :: running_moment(:)
 
integer, parameter:: n_time_bins=700000

 integer tbin, tbin_max/0./, lun
 integer i, n, lines/0/
 integer isign/1/
 integer m
 integer k, ix
 integer j
 integer nbins_max/149000/
 integer m_rev
 integer peak_bin

 real(rp) two/2./, zero/0./
 complex(rp), allocatable ::data_fft_n(:), data_fft_x(:), data_fft_y(:), data_fft_xx(:), data_fft_yy(:), data_fft_xxx(:), data_fft_yyy(:)
 complex(rp), allocatable:: cdata_fft_x(:),cdata_fft_y(:)
 real(rp) time_bin/1.e-8/
 real(rp), allocatable :: time(:)
 real(rp) freq 
 real(rp) freqs(10), x_betatron_tunes(10),y_betatron_tunes(10)
 real(rp) f_rev, Q_vert
 complex(rp) x_amps(10), y_amps(10)
 character*290 string
 character*40 file_name

allocate(moments(1:nbins_max), running_moment(1:nbins_max))
allocate(time(1:nbins_max))

lun=lunget()
   file_name = 'Time_Dep_Moments_at_END.dat'
   open(unit=lun, file = trim(file_name) )

   do while(.true.)
     read(lun,'(a)', end=99)string
     if(index(string,'<')/= 0)cycle
     lines = lines + 1
     if(lines > nbins_max)then
        print *,' Moment array has size smaller than number of bins in data file'
        stop
     endif
!     read(string,*,end=99)tbin, time, number, orb%vec(1:6)
! print '(a)',string
     read(string,'(i10)')i
     if(i > nbins_max)cycle
     read(string,'(i10,es18.8,i10,21(es12.4))',end=99)i, time(i),running_moment(i)%n, moments(i)%ave(:),&
                             (moments(i)%sigma(1+j,1+j),moments(i)%sigma(1+j,2+j),moments(i)%sigma(2+j,2+j), j=0,4,2), &
                             moments(i)%third(1:6)

!     if((lines/10000)*10000 == lines)print *, lines
     tbin_max = max(tbin_max, i)
  end do

 99 continue
   close(lun)

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

  allocate(data_fft_n(1:m), data_fft_x(1:m), data_fft_y(1:m), data_fft_xx(1:m), data_fft_yy(1:m), data_fft_xxx(1:m), data_fft_yyy(1:m))
  allocate(cdata_fft_x(1:m), cdata_fft_y(1:m))

   do i=1,m
    data_fft_n(i) = cmplx (running_moment(i)%n, zero)
    data_fft_x(i) = cmplx(moments(i)%ave(1),zero) * 2*sin((twopi/2*i)/m) * sin((twopi/2.*i)/m)
    data_fft_y(i) = cmplx(moments(i)%ave(3),0.) * 2*sin((twopi/2*i)/m) * sin((twopi/2.*i)/m)
    cdata_fft_x(i) = data_fft_x(i)
    cdata_fft_y(i) = data_fft_y(i)
    data_fft_xx(i) = cmplx(moments(i)%sigma(1,1),0.)
    data_fft_yy(i) = cmplx(moments(i)%sigma(3,3),0.)
    data_fft_xxx(i) = cmplx(moments(i)%third(1),0.0)
    data_fft_yyy(i) = cmplx(moments(i)%third(3),0.0)
   end do
    call four1(data_fft_n(1:m),isign)
    call four1(data_fft_x(1:m),isign)
    call four1(data_fft_y(1:m),isign)
    call four1(data_fft_xx(1:m),isign)
    call four1(data_fft_yy(1:m),isign)
    call four1(data_fft_xxx(1:m),isign)
    call four1(data_fft_yyy(1:m),isign)
    call naff(cdata_fft_x(1:m), x_betatron_tunes,x_amps)
    call naff(cdata_fft_y(1:m), y_betatron_tunes,y_amps)     
    
    print *,' tbin_max = ', tbin_max,' log(two) = ', log(two), ' n= ',n,' m = ',m,' size(data_fft_n) = ', size(data_fft_n)
    print '(a)', ' Horizontal tune '
    print  '(3es16.8)', (x_betatron_tunes(j), real(x_amps(j)), aimag(x_amps(j)), j=1,10)
    print '(a)',' Vertical Tune '
    print  '(3es16.8)', (y_betatron_tunes(j), real(y_amps(j)), aimag(y_amps(j)), j=1,10)

   lun = lunget()
   open(unit = lun, file = 'EndTime.dat')
 print *,' EndTime.dat'
 write(lun,'(a)')'#! Data from files '
 write(lun,'(a)')'#!',file_name

  write(lun, '(a9,1x,a12,7(3a14))')'time_bin','freq','Total Counts','Real fft','Imag fft',&
                                                 '<x>','Real fft','Imag fft',&
                                                 '<y>','Real fft','Imag fft',&
                                                 '<xx>','Real fft','Imag fft',&
                                                 '<yy>','Real fft','Imag fft',&
                                                 '<xxx>','Real fft','Imag fft',&
                                                 '<yyy>','Real fft','Imag fft'
     f_rev = 1./149.2e-3 !revolution frequency MHz
     m_rev = f_rev *1.e6*m*time_bin/2
     peak_bin=maxloc(abs(data_fft_y(1:m_rev)),1)
     Q_vert = float(peak_bin)/m_rev/2.
     print '(a,es12.4,a,i10,a,i10)', 'f_rev = ',f_rev,' m_rev = ',m_rev,' peak_bin = ', peak_bin 
     print *,' Q_vert = ', Q_vert
     do i = 1, m
     freq=float(i)/time_bin/m/1.e6
   write(lun,'(i10,es12.4,(i14,2es14.6),6(3es14.6))' ) i,freq,running_moment(i)%n, data_fft_n(i),&
                          moments(i)%ave(1), data_fft_x(i), &
                          moments(i)%ave(3), data_fft_y(i), &
                          moments(i)%sigma(1,1), data_fft_xx(i), &
                          moments(i)%sigma(3,3), data_fft_yy(i), &
                          moments(i)%third(1), data_fft_xxx(i), &
                          moments(i)%third(3), data_fft_yyy(i)

   end do
 
 end program