program analyze_harp
use precision_def
use nr
use bmad
use naff_mod
implicit none

interface
subroutine find_peak_freq(datac,nsize,freq)
 use precision_def
 implicit none
 real(rp) freq(:)
 complex(rp) datac(:)
 integer nsize
end subroutine


subroutine find_freq_near(f0,fh,fv,freq_in, freq)

 use precision_def
 implicit none
 real(rp) freq_in(:), freq(:), f0, fh, fv
 end subroutine
end interface

 integer, parameter:: n_time_bins=600000

 type harp_struct
  integer fiber(-3:3), ntotal, plane
  real(rp) average, rms, freq(100), third
  complex(rp) amp(20), fft_fiber(-3:3)
  real(rp) x
  
end type harp_struct

 type (harp_struct) harp(4,0:n_time_bins)
 type (harp_struct) tot_hits(4), FR(4), tune(4), width(4), fiber_tune(4,-3:3)

 type (coord_struct) start_orb
 integer fiberNum
 integer tbin, tbin_max, lun, lun1
 integer i, n, lines/0/
 integer isign/1/
 integer m
 integer k, ix
 integer opt_dump_spectra/79/
 integer delta_bin
 integer iy
integer tbin_m
 integer ix_timbin/1/

 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) 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(4)/1.25e-9,10.e-9,15.e-9,150.e-9/, radius/0.00025/
 real(rp) ave_last
 real(rp) closed_orbit
 real(rp) freq(3), f0/6.7/,fh/0.65/,fv/2.68/

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

 character*160 string, file(200)/200*''/
 character*6 wbin
 character*100 out_name
 logical active_fib(-3:3)/7*.false./
 logical unwrap/.false./, capture_count/.false./, harp_anal/.false./,compile_time/.false./
 character*50 dir_file_name
 namelist/input/time_bin,radius,start_turns,end_turns,  active_fib, unwrap, capture_count, harp_anal, turns_to_store, compile_time

 ! 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(:,:)%plane = 0
   harp(:,:)%x = 0
   harp(:,:)%rms = 0
   harp(:,:)%third = 0
   forall(i=1:4)tot_hits(i)%fiber(-3:3) = 0
   tbin_max = 0
ix=0
k=0
lun=lunget()

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

 if(unwrap)call unwrap_tgz

   lun1 = lunget()
   open(unit=lun1, file = 'harp_plane_hits_all.dat')

   k=0
   string = 'ls'//' */harp_plane_hits.dat > dirs.dat'
   print *,string
   call execute_command_line(string)

   ios=0
   open(unit=5, file = 'dirs.dat')
   do while(ios == 0)
    read(5,'(a)',iostat=ios)dir_file_name
   print *, dir_file_name
    if(ios== 0)then
     k=k+1
     file(k) = trim(dir_file_name)
     print *,' file = ',trim(file(k))
   end if
   end do
  close(unit=5)

!do k=1,size(file)
! if(trim(file(k)) /= '')print *,' file = ',trim(file(k))
!end do

if(capture_count)call capture(file,k, turns_to_store)

if(compile_Time)call compile_Time_Dep(file)

if(harp_anal)then
    do i=-3,3,1
     if(active_fib(i))num_fib=num_fib+1
    end do

do ix_timbin=1,4
  if(time_bin(ix_timbin) == 0.)cycle

do k=1, 3 !size(file)
   if(trim(file(k)) == '')cycle  
   lun=lunget()
   open(unit=lun, file = file(k), action = 'READ')
   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
!     write(lun1,'(i5,8es18.8)')fiberNum, start_orb%t, start_orb%s, start_orb%vec
     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
    ix_hit = -99
    iy_hit = -99

    !never mind fibers. Count every plane hit

    tbin =(start_orb%t - start_turns*149.e-9)/time_bin(ix_timbin)+.5
        if(tbin > n_time_bins)then
          print *,' tbin = ',tbin ,'> n_time_bins = ', n_time_bins
          stop
        endif
         tbin_max = max(tbin, tbin_max)
            harp(fiberNum,tbin)%plane = harp(fiberNum,tbin)%plane+1
            if(fiberNum == 2 .or. fiberNum == 4)harp(fiberNum,tbin)%x = harp(fiberNum,tbin)%x+start_orb%vec(1)
            if(fiberNum == 1 .or. fiberNum == 3)harp(fiberNum,tbin)%x = harp(fiberNum,tbin)%x+start_orb%vec(3)
          

         
    do i=-3,3,1
     if(active_fib(i))then
     if((fiberNum == 1 .or. fiberNum == 3) .and. (abs(start_orb%vec(3) - yfiber(i))< radius))iy_hit=i
     if((fiberNum == 2 .or. fiberNum == 4) .and. (abs(start_orb%vec(1) - xfiber(i))< radius))ix_hit=i
        if((iy_hit > -99 .and. (fiberNum == 1 .or. fiberNum == 3)) .or. (ix_hit>-99 .and.(fiberNum == 2 .or. fiberNum== 4)))then

         harp(fiberNum,tbin)%fiber(i) = harp(fiberNum,tbin)%fiber(i) + 1
         harp(fiberNum,tbin)%ntotal = harp(fiberNum,tbin)%ntotal +1
         n= harp(fiberNum,tbin)%ntotal
         tot_hits(fiberNum)%fiber(i) = tot_hits(fiberNum)%fiber(i) + 1
         if(ix_hit > -99)then
             ave_last = harp(fiberNum,tbin)%average
             harp(fiberNum,tbin)%average = ((n-1)*harp(fiberNum,tbin)%average + xfiber(i))/float(n) 
             harp(fiberNum,tbin)%rms = ((n-1) * harp(fiberNum,tbin)%rms + xfiber(i)**2 +(n-1)*ave_last**2 - n * harp(fiberNum,tbin)%average**2 )/float(n)
         endif
         if(iy_hit>-99)then
             ave_last = harp(fiberNum,tbin)%average
             harp(fiberNum,tbin)%average = ((n-1)*harp(fiberNum,tbin)%average + yfiber(i))/float(n) 
             harp(fiberNum,tbin)%rms = ((n-1) * harp(fiberNum,tbin)%rms + yfiber(i)**2 +(n-1)*ave_last**2 - n * harp(fiberNum,tbin)%average**2)/float(n)
         endif
        exit
        endif


      endif  !active_fib
    end do
  end do


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

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

  if(.not.allocated(datac))allocate(datac(1:4,1:m), datac_ave(1:4,1:m), datac_rms(1:4,1:m), datac_third(1:4,1:m))
  
  do fiberNum=1,4
   if(fiberNum == 1 .or. fiberNum == 3) closed_orbit=sum(tot_hits(fiberNum)%fiber(:)*yfiber(:))/sum(tot_hits(fiberNum)%fiber)  
   if(fiberNum == 2 .or. fiberNum == 4) closed_orbit=sum(tot_hits(fiberNum)%fiber(:)*xfiber(:))/sum(tot_hits(fiberNum)%fiber)  
   do i=1,m
    i_abs = abs(i)
    datac(fiberNum,i) = cmplx(float(harp(fiberNum,i_abs)%ntotal),0.)
    datac_ave(fiberNum,i) = cmplx(harp(fiberNum,i_abs)%average - closed_orbit,0.)
    datac_rms(fiberNum,i) = cmplx(harp(fiberNum,i_abs)%rms,0.)
    harp(fiberNum,i)%fft_fiber(:) = cmplx(harp(fiberNum, i)%fiber(:),0.)
   end do
   print *,' tbin_max = ', tbin_max,' log(two) = ', log(two), ' n= ',n,' m = ',m,' size(data_cos(fiberNum,:) = ', size(datac(fiberNum,:))
 ! call naff(datac(fiberNum,1:m),tot_hits(fiberNum)%freq,tot_hits(fiberNum)%amp,opt_dump_spectra)
  call four1(datac(fiberNum,1:m),isign)
  call find_peak_freq(datac(fiberNum,1:m),m,FR(fiberNum)%freq)
  call four1(datac_ave(fiberNum,1:m),isign)
  call find_peak_freq(datac_ave(fiberNum,1:m),m,Tune(fiberNum)%freq)
  call four1(datac_rms(fiberNum,1:m),isign)
  call find_peak_freq(datac_rms(fiberNum,1:m),m,Width(fiberNum)%freq)
  do i=-3,3
   if(active_fib(i))then
   call four1(harp(fiberNum,1:m)%fft_fiber(i),isign)
   call find_peak_freq(harp(fiberNum,1:m)%fft_fiber(i),m,fiber_tune(fiberNum,i)%freq)
   endif !active_fib
  end do
 end do
  write(wbin,'(f6.2)')time_bin(ix_timbin)*1.e9
if(wbin(1:1) == ' ')wbin(1:1) = '0'
if(wbin(2:2) == ' ')wbin(2:2) = '0'

   lun = lunget()
if(num_fib > 1) out_name = 'FiberHarp-'//wbin//'ns.dat'
if(num_fib == 1) out_name = 'FiberHarp-'//wbin//'ns.dat_central_fib'
  open(unit = lun, file = trim(out_name))
 print *, trim(out_name)
 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(datac(fiberNum,:) = ', size(datac(fiberNum,:))
 write(lun,'(a1,1x,a,es12.4)')'#',' Time bin = ',time_bin(ix_timbin)
 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
  do fiberNum = 1,4
    if(fiberNum == 1 .or. fiberNum == 3)then
       write(lun,'(a7,1x,i1,1x,a,es12.4)')'# Harp ', fiberNum, 'average_positon = ', sum(tot_hits(fiberNum)%fiber(:)*xfiber(:))/sum(tot_hits(fiberNum)%fiber)  
     else
       write(lun,'(a7,1x,i1,1x,a, es12.4)')'# Harp ', fiberNum, 'average_positon = ', sum(tot_hits(fiberNum)%fiber(:)*yfiber(:))/sum(tot_hits(fiberNum)%fiber)    
    endif
       write(lun,'(a,1x,i1,1x,a,10es12.4)')'# Harp ',fiberNum, ' FR frequencies ', FR(fiberNum)%freq(1:10)/m/time_bin(ix_timbin)/1.e6
       write(lun,'(a,1x,i1,1x,a,10es12.4)')'# Harp ',fiberNum, ' Tune frequencies ', Tune(fiberNum)%freq(1:10)/m/time_bin(ix_timbin)/1.e6
       write(lun,'(a,1x,i1,1x,a,10es12.4)')'# Harp ',fiberNum, ' Width frequencies ', Width(fiberNum)%freq(1:10)/m/time_bin(ix_timbin)/1.e6
  end do

  write(lun,'(a1,1x,2a18,4a18,3a36,a18)')'#','Harp','Time bin','Average position','RMS position','Total Hits','Freq bin', &
                                               'FFT cmplx tot hits','FFT cmplx average','FFT complex rms','plane avg position'
  do fiberNum =1,4
   do tbin = 1, tbin_max
      tbin_m = min(m, tbin)
!    if(harp(fiberNum,tbin)%ntotal /= 0)write(lun, '(2i12,2es12.4,i10)')fiberNum, tbin, harp(fiberNum,tbin)%average, harp(fiberNum,tbin)%rms, harp(fiberNum,tbin)%ntotal
     if(harp(fiberNum,tbin)%plane > 0)write(lun, '(2x,2i18,2es18.4,i18,1x,es18.4, 7es18.4)')fiberNum, tbin,  harp(fiberNum,abs(tbin))%average, harp(fiberNum,abs(tbin))%rms, &
                   harp(fiberNum,abs(tbin))%ntotal, abs(tbin)/time_bin(ix_timbin)/m/1.e6,&
                   datac(fiberNum,tbin_m), datac_ave(fiberNum, tbin_m), datac_rms(fiberNum,tbin_m), harp(fiberNum,tbin)%x/harp(fiberNum,tbin)%plane
   end do

  end do
 close(lun)
 close(lun1)

! data file to plot phase space - match of average position at harps 1&3 in tbins off by 1/4 149 and same for 2&4
delta_bin = 149.e-9/time_bin(ix_timbin)/4
lun=lunget()
open(unit=lun,file='fiberPhaseSpace.dat')
print *,' Write fiberPhaseSpace.dat'
write(lun, '(2a,1x,i10)')' Phase -space: ', 'delta_bin = ',delta_bin
write(lun, '(4a16)')' harp','time bin',' pos_ave (harp)','pos_ave( +2, delta_bin)'
do fiberNum = 1,2
do tbin = 1,tbin_max  !tbin in time_bin increments. Delta_bin = 149.e-9/time_bin/4
  write(lun,'(i10,6x,i10,6x,2(es12.4,4x))') fiberNum, tbin, harp(fiberNum,abs(tbin))%average, harp(fiberNum+2, abs(tbin + delta_bin))%average 
end do
end do
close(lun)
 
 lun = lunget()
 open(unit = lun, file = 'Fiber-'//wbin//'ns.dat')
 lun1 = lunget()
 open(unit = lun1, file = 'FiberChromaticity.dat')
 print *,' Writing Individual Fiber-'//wbin//'ns.dat'
 write(lun,'(a)')'#! Data from files '

  write(lun, '(2x,a18,2x,a18,5a18)')'Harp','fiber','time bin','freq bin','total hits','real(fft)','imag(fft)'
  do fiberNum =1,4
   do iy=-3,3
    if(active_fib(iy))then
    write(lun,'(a1,2x,2i10,10es14.6)')'#',fiberNum,iy,fiber_tune(fiberNum,iy)%freq(1:10)/m/time_bin(ix_timbin)/1.e6
    do tbin = 1, m
!     write(6, *)fiberNum, iy, abs(tbin)/time_bin(ix_timbin)/m/1.e6, &
!                   harp(fiberNum,tbin)%fiber(iy), &
!                   harp(fiberNum,tbin)%fft_fiber(iy)

     write(lun, '(2x,i18,2x,i18,2x,i18,es18.4,2x,i18,2x,2es18.4)')fiberNum, iy,tbin, abs(tbin)/time_bin(ix_timbin)/m/1.e6, &
                   harp(fiberNum,tbin)%fiber(iy), &
                   harp(fiberNum,tbin)%fft_fiber(iy)
    end do
    call find_freq_near(f0,fh,fv,fiber_tune(fiberNum,i)%freq/m/time_bin(ix_timbin)/1.e6, freq)
    if(freq(1) > 0)then
     print '(2i10,3es12.4)',fibernum,i, freq(1),freq(2),freq(2)/freq(1)
     write(lun1, '(2i10,3es12.4)')fibernum,i, freq(1),freq(2),freq(2)/freq(1)
    endif

  endif  !active_fib
  end do
  end do

end do ! time bins
endif !harp_anal 
 end program


   
subroutine find_peak_freq(datac,nsize,freq)
 use precision_def
 implicit none
 real(rp)  freq(:)
 complex(rp)  datac(:)
 real(rp), allocatable :: temp(:)
 real(rp) x,a,b,c
 integer i, n,j, nsize

! nsize = size(datac)
 allocate(temp(0:nsize+1))

!print *, ' nsize = ',nsize
n=0
do i=1,nsize 
temp(i) = sqrt(real(datac(i))**2+imag(datac(i))**2)
end do
if(all(temp(:) == 0))return
temp(1:10)=0
temp(nsize-10:nsize)=0
 temp(nsize-min(5,nsize):nsize)=0
 do while  (n<100)
  j=maxloc(temp,1)
  n=n+1
  freq(n) = j

  c = temp(j)
  a = 0.5*(temp(j+1)+temp(j-1)-2*temp(j))
  b = 0.5*(temp(j+1)-temp(j-1))
  x=0
  if(a > 0)x = -b/2/a
 print '(i10,1x,es12.4,1x,i10,1x,5es12.4)',j,temp(j),n,freq(n), (x+j)/nsize*100.,a,b,c
  freq(n) =j+ x
  temp(j-5:j+5) = 0
end do

deallocate(temp)
return
end subroutine

 subroutine find_freq_near(f0,fh,fv,freq_in, freq)

 use precision_def
 implicit none
 real(rp) freq_in(:), freq(:), f0, fh, fv

 integer i

 i =  minloc(abs(freq_in(:)-f0),1)
 freq(1) =  freq_in(i) 
 i =  minloc(abs(freq_in(:) - fh),1)
 freq(2) =  freq_in(i) 
 i =  minloc(abs(freq_in(:) - fv),1)
 freq(3) =  freq_in(i) 

 return
 end subroutine

subroutine capture(file, kfiles, turns_to_store)
use bmad
implicit none
integer, parameter :: nmax=5000
 integer  lun, lun1
 integer i, n
 integer k
 integer ios
 integer ix(nmax), number_remaining(nmax),   number_on_turn(0:nmax), incident(200)
 integer total
 integer total_inc
 integer kfiles
 integer ixx
 integer turns_to_store

 character*220  string
 character*160 file(200), local(200), dir
 character*16 element(nmax) 
real(rp)   turn(nmax), x_ave(nmax), y_ave(nmax), t_ave(nmax), vec(nmax,1:10)

   lun1 = lunget()
   open(unit=lun1, file = 'capture.dat')

do k=1,size(file)
dir = file(k)
ixx=index(dir,'/')
 local(k) = dir(1:ixx-1)
 if(trim(local(k)) /= '')print '(a10,a19)',' file = ',trim(local(k))
end do
total_inc=0
total=0
do k=1,size(file)
   if(trim(local(k)) == '')cycle  
   lun=lunget()
   local(k) = trim(local(k))//'/Beam_moments_tbt.dat'
   open(unit=lun, file = trim(local(k)))
   print '(a10,a40)',' file = ', trim(local(k))
   i=0
   do while(.true.)
     read(lun,'(a)', end=99)string

     if(index(string,'Element')/= 0)cycle
     i=i+1
     read(string,*,end=99)ix(i), turn(i), x_ave(i), y_ave(i), t_ave(i), vec(i,1:10), number_remaining(i), element(i) 
     if(element(i) == 'BEGINNING')incident(k) = number_remaining(i)
!if(i == 11)then
!  print '(i10,14es12.4,1x,i10,1x,a16)',ix(i), turn(i), x_ave(i), y_ave(i), t_ave(i), vec(i,1:10), number_remaining(i), element(i) 
!  print '(a)', string
!endif
     number_on_turn(int(turn(i))) = number_remaining(i)
   end do

 99 continue
if(i>0) write(lun1,'(a,1x,2i10,1x,es12.4)')trim(file(k)), number_on_turn(turns_to_store), incident(k), float(number_on_turn(turns_to_store))/incident(k)
   close(lun)
 total = number_on_turn(10) + total
 total_inc = incident(k)+total_inc
end do
write(lun1,'(a,1x,2i10,1x,es12.4)')'total', total, total_inc, float(total)/total_inc


return

end subroutine

subroutine compile_Time_Dep(file)

  use bmad
  use muon_mod
implicit none

TYPE (g2moment_struct),save, allocatable :: moment(:) !moment structure for each time bin, 1000 turns => an array with 1000*149e-9/bin_width 
type(g2moment_struct), save, allocatable :: sum_moment(:) 
!type running_moment_struct
!  real(rp) product(1:6,1:6), sum(1:6), third(6)
!  integer n
!end type
type (running_moment_struct),save, allocatable :: running_moment(:)
type (running_moment_struct), save, allocatable :: sum_running_moment(:)
real(rp), save, allocatable :: ptop(:),sum_ptop(:), time(:)
real(rp) y, x, sx,sxx,syx,sy,a,b, det

integer, parameter :: nmax=5000
 integer  lun, lun1
 integer i, n,j
 integer k
 integer ios
 integer kfiles
 integer ixx
 integer nturns/4200/
 integer imax/0/
 
 character*400  string
 character*160 file(200), local(200), dir
 character*16 element(nmax) 

real(rp) bin_width/149.e-9/

logical first/.true./

if(first)then
   
 n = nturns*149.e-9/bin_width + 100
 allocate(moment(0:n),running_moment(0:n),ptop(0:n))
  allocate(sum_moment(0:n),sum_running_moment(0:n), sum_ptop(0:n), time(0:n))
  first = .false.
  do i=1,n
     sum_moment(i)%ave=0
     sum_running_moment(i)%n=0
     sum_ptop(i)=0
  end do
endif
 
  lun1 = lunget()
 open(unit=lun1, file = 'time_dep_moments_all.dat')

do k=1,size(file)
dir = file(k)
ixx=index(dir,'/')
 local(k) = dir(1:ixx-1)
 if(trim(local(k)) /= '')print '(a10,a19)',' file = ',trim(local(k))
end do

do k=1,size(file)
   if(trim(local(k)) == '')cycle  
   lun=lunget()
   local(k) = trim(local(k))//'/Time_Dep_Moments_at_END.dat'
   open(unit=lun, file = trim(local(k)))
   print '(a10,a40)',' file = ', trim(local(k))
   i=0
   do while(.true.)
     read(lun,'(a)', end=99)string

     if(index(string,'Time')/= 0)cycle
  !   i=i+1

     read(string,*,end=99)i, time(i), running_moment(i)%n, moment(i)%ave(:),&
          moment(i)%sigma(1,1),moment(i)%sigma(1,2),moment(i)%sigma(2,2), &
          moment(i)%sigma(3,3),moment(i)%sigma(3,4),moment(i)%sigma(4,4), &
          moment(i)%sigma(5,5),moment(i)%sigma(5,6),moment(i)%sigma(6,6), &
                                       moment(i)%third(:), ptop(i)   
     sum_running_moment(i)%n = sum_running_moment(i)%n+running_moment(i)%n
     sum_moment(i)%ave(:) =sum_moment(i)%ave(:) +moment(i)%ave(:)*running_moment(i)%n/sum_running_moment(i)%n
     sum_ptop(i) = sum_ptop(i)+ptop(i)*running_moment(i)%n/sum_running_moment(i)%n
     imax = max(imax,i)
  end do

 99 continue
end do
   !fit to sum_ptop out to where  n < 100
  j=0
  sx=0
  sy=0
  syx=0
  sxx=0
  
  do i =1,imax
     if(sum_running_moment(i)%n < 100 .or. sum_ptop(i) == 0)cycle
     j=j+1
     y = log(sum_ptop(i))
     x = time(i)
     sx = sx + x
     sy = sy+ y
     sxx = sxx + x**2
     syx = syx+ y*x
     print '(2es12.4,10x,es12.4)',sx,sxx,syx
     print '(i12,es12.4,10x,es12.4)',j,-sx,sy
  end do
  det = -sx*sx + sxx * j
  a = (-sx * syx + sxx * sy)/det
  b = (-j*syx + sx * sy)/det
  print '(a,es12.4,a,es12.4)',' cbo decay time = ',1./b,' y0 = ',exp(a)
  
  write(lun1,'(a1,a20,es12.4,a6,es12.4)')'#',' cbo decay time = ', 1./b,'y0= ',exp(a)
  write(lun1,'(a1,a9,a12,a11,6a12,a14)')'#','bin','time','n','<x>','<xp>','<y>','<yp>','<z>','<dp>','<peak-to-peak>'
  do i=1,imax
   write(lun1,'(i10,1x,es12.4,1x,i10,1x,7es12.4)')i,time(i),sum_running_moment(i)%n,sum_moment(i)%ave(:),sum_ptop(i)
end do
print *,' time_dep_moments_all_time.dat  written'

end subroutine