subroutine compute_moments_vs_time_calo(nturns, turn, coord, calo)

! In BMAD tracking, the coordinate along the reference trajectory is the dependent variable.
!The tracking routines, propagate the coordinates of the particles in a distribution from 
!the end of one element to the end of the next.
!The phase space coordinates, and the time, are given with respect to the reference particle.
!Detectors in the g-2 ring will measure the properties of the distribution (centroid, width, etc.)
!at a particular location and within some window of time. It is therefore convenient to compute
!moments of the distribution within that window.
! In this routine the phase space coordinates at the end of an element are sorted by arrival time.
! Running moments are accumulated.

use muon_mod
use parameters_bmad
implicit none

type (coord_struct) coord
integer i
integer lun
integer, save :: n
integer j,k
integer turn
integer time_bin
integer nturns
integer ncalo
integer word_length


real(rp) num
real(rp) x_max, x_min
real(rp), allocatable :: ptop(:)

character*16 calo, number, name(1:24)/24*'NULL'/

logical first/.true./

TYPE (g2moment_struct), allocatable :: moment(:) !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), save, allocatable :: calo_running_moment(:,:)

if(turn > 0) then
   
if(bin_width == 0)return
nturns = number_turns
word_length = len(trim(calo))
number = calo(6:word_length)
if(len(trim(number)) == 1)read(calo(6:6),*)ncalo
if(len(trim(number)) == 2)read(calo(6:7),*)ncalo
name(ncalo) = calo
 
   if(first)then
      print *,'compute_moments_vs_time_calo: ncalo ',ncalo, name(ncalo)
      if(ncalo < 1 .or. ncalo > 24) then
         print *,'compute_moments_vs_time_calo: this calo number ', ncalo,' does not make sense'
         stop
      endif

   n = nturns*149.e-9/bin_width + 100
   print *,' number of time bins = ',n
   allocate(calo_running_moment(0:n,0:24))
   do i=1,n
      do j=1,24
         calo_running_moment(i,j)%product(:,:)=0
         calo_running_moment(i,j)%sum(:)=0
         calo_running_moment(i,j)%n=0
         calo_running_moment(i,j)%third(:) = 0
      end do    
   end do
   first=.false.
 endif

!Sort into time bins

     if(coord%state /= alive$)return ! .and. muons(i)%decay == .false.)cycle
     time_bin = coord%t/bin_width +0.5 ! associate integer with each bin_width time slice
     if(time_bin > n .or. time_bin <0)return
     ! running sums
     do j=1,6
        calo_running_moment(time_bin, ncalo)%product(1:6,j) = calo_running_moment(time_bin, ncalo)%product(1:6,j) + coord%vec(1:6)*coord%vec(j)
     end do
     

      !      print '(3i12,6es12.4)',i,j,time_bin,running_moment(time_bin)%product(1:6,j)
     calo_running_moment(time_bin,ncalo)%sum(1:6) = calo_running_moment(time_bin,ncalo)%sum(1:6) + coord%vec(1:6)
     calo_running_moment(time_bin,ncalo)%n = calo_running_moment(time_bin,ncalo)%n+1
     calo_running_moment(time_bin,ncalo)%third(1:6) = calo_running_moment(time_bin,ncalo)%third(1:6) +coord%vec(1:6)**3
  
  

   
elseif(turn < 0)then !compute moments and write. <turn> is the flag that tracking is complete
do ncalo = 1,24
if(name(ncalo) == 'NULL')cycle
! moments
   allocate(moment(0:n), ptop(0:n))
   ptop(:)=0.
   
  print *,' moments, n= ', n
     do i=1,n
        !     time_bin = muons(i)%coord%t/bin_width +0.5 ! associate integer with each bin_width time slice
        if(calo_running_moment(i,ncalo)%n == 0)cycle
        forall(j=1:6)moment(i)%ave(j) = calo_running_moment(i,ncalo)%sum(j)/float(calo_running_moment(i,ncalo)%n)
        num = float(calo_running_moment(i,ncalo)%n)
        do j=1,6
           do k=1,6
              moment(i)%sigma(j,k) = calo_running_moment(i,ncalo)%product(j,k)/num  -  moment(i)%ave(j)*moment(i)%ave(k)
           end do
           moment(i)%third(j) = calo_running_moment(i,ncalo)%third(j)/num -3*calo_running_moment(i,ncalo)%product(j,j)*moment(i)%ave(j)/num + 2*moment(i)%ave(j)**3
        end do
     end do
  

 lun=lunget()
 open(unit=lun,file=trim(directory)//'/'//'Time_Dep_Moments_at_'//trim(name(ncalo))//'.dat')
 print *,'Write: Time_Dep_Moments_at_'//trim(name(ncalo))//'.dat'
 write(lun,'(a10,a18,a10,22a12)')'Time bin','time','number','<x>','<px>','<y>','<py>','<z>','<pz>','<xx>','<x-px>','<px-px>','<yy>','<y-py>','<py-py>', &
                                    '<zz>','<z-pz>','<pz-pz>','<xxx>','<xpxpxp>','<yyy>','<ypypyp>','<zzz>','<zpzpzp>','<ptop>'


 do i=1,n
  if(calo_running_moment(i, ncalo)%n == 0)cycle

  if(i>15)then
   x_max = maxval(moment(i-15:i)%ave(1))
   x_min = minval(moment(i-15:i)%ave(1))
   ptop(i) = x_max-x_min
 endif

  write(lun,'(i10,es18.8,i10,21es12.4,es12.4)')i, i*bin_width,calo_running_moment(i, ncalo)%n, moment(i)%ave(:),&
                             (moment(i)%sigma(1+j,1+j),moment(i)%sigma(1+j,2+j),moment(i)%sigma(2+j,2+j), j=0,4,2), &
                             moment(i)%third(1:6), ptop(i)
 end do
 close(unit=lun)
!  call cbo_decay_time(ptop,bin_width, decay_time)


deallocate(moment,ptop)
end do !ncalo
endif

return
end