! Find all  files 'domega_ddelta_vs_time.dat'  with columns 'time, dphi/de, err'
! Combine and write 'combined_domega_ddelta_vs_time.dat' with 'time, average(dphi/de), combined error'

program compile_domega_vs_time
  use sim_utils
  implicit none
  type time_phase_struct
     real(rp) time(20000),dphide(20000),err(20000)
     integer ix(20000)
  end type time_phase_struct
  type (time_phase_struct), allocatable :: time_phase(:)
  character*300 dir_file
  character*300 dir(300)
  character*20 dir_min, dir_max
  character*200 cwd
  character*200 string, new_string
  character*2 cnum
  character*100 file_string
  real(rp) ave_dphide,err
  real(rp) slice_time, dphide, junk
  
  logical itexists, first/.true./

  integer reason, i, j, n_left/0/, n_left_0/0/
  integer lun, lun2
  integer number_files
  integer ngood_file
  integer n,m, time_bins,nn
  integer k
  integer nevents/1000000/
  integer min_date,max_date, min_time, max_time, date, time
  integer nargs
  integer ix
  
  dir(1)=' '
 nargs = command_argument_count()
 if (nargs == 1)then
  call cesr_getarg(1,dir(1))    ! the data is all in a single directory input directory name on command line
  !  print *, 'len(trim(dir(1)))= ',len(trim(dir(1)))
    number_files=1
    print '(i10,1x,a19)', number_files,dir(1)
 else
    if( nargs == 2)then
       call cesr_getarg(1,dir_min)
       call cesr_getarg(2,dir_max)
       read(dir_min(1:8),*)min_date
       read(dir_min(10:15),*)min_time
       read(dir_max(1:8),*)max_date
       read(dir_max(10:15),*)max_time
    endif
    
   call getcwd(cwd)      !if all directories are to be combined, input nothing on command line
   string='ls > out.dat'
   call execute_command_line(string)
   lun=lunget()
   open(unit=lun,file='out.dat', status='old')
   n=0
  do while(.true.)
    read(lun, '(a)', IOSTAT = reason)new_string
    if(reason < 0)exit
    if(  (index(new_string,'2025')==0) .and. (index(new_string,'2019')==0)  .and. (index(new_string,'2021')==0)  .and. (index(new_string,'2022')== 0) .and. &
                                      index(new_string,'2023')==0 .and.index(new_string,'2024')==0)cycle
    n=n+1
    dir(n)=trim(new_string)
    print '(i10,1x,a19)',n, dir(n)
    number_files = n
  end do
  close(lun)
 endif
!_______________________________________________________________________

allocate(time_phase(1:number_files))
 ngood_file=0
 k=0
 do n=1, number_files
    if(nargs == 2)then
       read(dir(n)(1:8),*)date
       read(dir(n)(10:15),*)time
       if(date<min_date)cycle
       if(time<min_time)cycle
       if(date>max_date)cycle
       if(time>max_time)cycle
    endif

    string='ls '//trim(dir(n))//' > '//trim(dir(n))//'/out.dat'
    call execute_command_line(string)
    lun=lunget()
    inquire (file = trim(dir(n))//'/'//'out.dat', exist = itexists)
    if(.not. itexists)then
       print *, trim(dir(n))//'/'//'out.dat',  'does not exist'
       cycle
    endif

    open(unit=lun,file=trim(dir(n))//'/'//'out.dat', status='old')

    if(first)then
       lun2=lunget()
       open (unit=lun2,file = 'combined_domega_ddelta_vs_time.dat')
       write(lun2,'(3a16)')'time', '<dphi/de>',' err'
       first=.false.
    endif
    

    do while(.true.)
       read(lun, '(a)', IOSTAT = reason)new_string
       if(reason < 0)exit
       if(index(new_string,'domega_ddelta_vs_time.dat')/=0)then
          dir_file = trim(dir(n))//'/'//trim(new_string)
          exit
       endif
    end do
    close(unit=lun)


    print '(a,a)',' dir_file =',dir_file
    cnum='00'
    if(n<10)write(cnum(2:2),'(i1)')n
    if(n>=10)write(cnum,'(i2)')n
!      print *, 'n=',n,' cnum=',cnum
!      file_string = 'f'//cnum//' = "'//trim(dir_file)//'"'
!      print *,'file_string =',file_string
!    endif

    write(401,'(a)')'f'//cnum//' = "'//trim(dir_file)//'"'
    ngood_file=ngood_file + 1

print *,' ngood_file = ', ngood_file

    lun=lunget()
    open(unit=lun, file = dir_file)
    m=0
    time_phase(ngood_file)%ix(:) =0    
    do while(.true.)
       read(lun,'(a)',IOSTAT=reason)new_string
       if(index(new_string,'time')/=0)cycle
       if(reason<0)exit
       m=m+1
!      print '(a)',trim(new_string)
       read(new_string,*)slice_time,dphide,err,junk,ix
!       read(new_string,*)time_phase(ngood_file)%time(m),time_phase(ngood_file)%dphide(m),time_phase(ngood_file)%err(m), time_phase(ngood_file)%ix(m)
       time_phase(ngood_file)%time(ix)=slice_time
       time_phase(ngood_file)%dphide(ix)=dphide
       time_phase(ngood_file)%err(ix)=err
       time_phase(ngood_file)%ix(m) =ix

!       print '(2i10,3es12.4)',ngood_file,m, time_phase(ngood_file)%time(m),time_phase(ngood_file)%dphide(m),time_phase(ngood_file)%err(m)
    end do
    close(unit=lun)
    time_bins=m

print '(a,i10)','time_bins =', time_bins
 end do  !loop over  files (n)

 
 do m = 1, time_bins
    if(any(time_phase(1:ngood_file)%ix(m) == 0))then
       do n=1,ngood_file
          if(time_phase(n)%ix(m) == 0)then
             print '(a,i10,a,i10,a,i12)', ' file=',n, ' m=',m,' time_phase(n)%ix(m) =',time_phase(n)%ix(m)
          endif
       end do
       cycle
    endif
    do n=1, ngood_file
       ix = time_phase(n)%ix(m)

       if(time_phase(n)%time(ix) /= time_phase(1)%time(ix) .and. time_phase(n)%time(ix) /= 0 .and. time_phase(1)%time(ix) /= 0)then
          print '(a,a3,1x,i10,a3,1x,i10,a6,i10,a,es12.4)',' mismatch ',' n=',n,' m=',m,'ix=',ix,' time_dif=',time_phase(n)%time(ix)-time_phase(1)%time(ix)
          print '(a,es12.4,a,es12.4)',' time_phase(n)%time(ix) =',time_phase(n)%time(ix),' time_phase(1)%time(ix)=',time_phase(1)%time(ix)
          stop
       endif
    end do
    ave_dphide = sum(time_phase(1:ngood_file)%dphide(ix))/ngood_file
    err = sqrt(sum(time_phase(1:ngood_file)%err(ix) **2)) /ngood_file
    err = sqrt(sum((time_phase(1:ngood_file)%err(ix) /ngood_file)**2))
    write(lun2,'(3es12.4)')time_phase(1)%time(ix),ave_dphide, err
 end do
 close(unit=lun2)
 print '(a)',' write combined_domega_ddelta_vs_time.dat'
 print '(a)',' fort.401 lists directories. load to gnuplot' 
 !   call execute_command_line('cat '//trim(file_string)//' > all_Energy_vs_time_0.dat')

end program compile_domega_vs_time