program test_compute_moments_vs_time
  use muon_mod
  use muon_interface
 use parameters_bmad
 use nr
 use sim_utils

 implicit none

type (muon_struct), allocatable :: muons(:)
type (kicker_params_struct) kicker_params
  integer vparam_id, vparam_col_ind, nmuons
  character*20 vparam_label 

 integer i, tot, lun
 integer time_bin
 integer nturns
 integer turn
 integer bin_min, bin_window !minimum bin and window for fitting to x vs t
 integer num
 integer ndf
 integer ios
 integer ndf_save, ndf_old

 real(rp) ptop_start_time
 real(rp) coef(4), chisq, chisq_pdf
  real(rp) coef_save(4), chisq_save, coef_old(4), chisq_old

 

 character*100 string/'TEST'/
 character*100 RFvparam_file/'RFquadparamDependence.dat'/
 character*400 read_string
 
integer test
logical itexists

nmuons = 9000

allocate(muons(nmuons))

 turn = -1
 nturns=400
 string='TEST'
 directory='test'
 bin_width=149.2e-9

 inquire (file= RFvparam_file, exist=itexists)

 if(itexists)then
    lun = lunget()
    open(unit = lun, file= RFvparam_file,STATUS='old',IOSTAT=ios)
    read_string=' '
    do while(index(read_string,'loop') == 0)
!       print *,' index(read_string,"loop") ', index(read_string,'loop')
       read(lun,'(a)')read_string
       if(index(read_string,'vparam_id =')/=0)read(read_string(12:),'(i4)')vparam_id
       if(index(read_string,'"')/=0)read(string(index(read_string,'"'):),'(a)')vparam_label
       if(index(read_string,'column1 =')/=0)read(read_string(10:),'(i5)')vparam_col_ind
    end do
    do while(ios >=0)
       read(lun, '(a)',IOSTAT = ios)read_string
       if(ios < 0)cycle
!       read(lun,'(i12,4es12.4,4es12.4,4es12.4,4es12.4,4es12.4,es12.4,i12,es12.4,2i12)')  &
       read(read_string,*)  num, rf_quad(1:4)%amp_h,rf_quad(1:4)%phi_h, rf_quad(1:4)%start_h, rf_quad(1:4)%freq_h,   coef_old(1:4), chisq_old, ndf_old, chisq_pdf,   nmuons,  tot
       rf_quad(1:4)%start_h  = rf_quad(1:4)%start_h * 1.e6
       rf_quad(1:4)%freq_h  =  rf_quad(1:4)%freq_h/1.e6
       call compute_moments_vs_time(nturns, turn,nmuons, muons, string,coef_save = coef,chisq_save=chisq,ndf_save=ndf, test=num )
       write(6,'(i12,4es12.4,4es12.4,4es12.4,4es12.4,4es12.4,es12.4,i12,es12.4,2i12)')    num, rf_quad(1:4)%amp_h,rf_quad(1:4)%phi_h, rf_quad(1:4)%start_h, rf_quad(1:4)%freq_h,   coef(1:4), chisq, ndf, chisq_pdf,   nmuons,  tot
       call write_rfquad_loop (vparam_id, vparam_label, vparam_col_ind, nmuons,kicker_params, tot, coef, chisq, ndf)
    end do
    
 else
    call compute_moments_vs_time(nturns, turn,nmuons, muons, string, test=0)
endif
end program test_compute_moments_vs_time