program count_zeros

use bmad
implicit none
character*250 string
character*60 file_name
character*20 word(15)
character*20 n_max_word, n_group_word
integer nargs,iargc,i,j,n
integer n_group/20/,nzero, n_max
real(rp), allocatable :: time(:), pol(:)
real(rp), allocatable::time_zero(:)
real(rp) freq,freq_0/2.29078556E+05/

  nargs = command_argument_count()
  if(nargs >= 1)then
  call get_command_argument(1,file_name)
  call get_command_argument(2,n_max_word)
  call get_command_argument(3,n_group_word)
  read(n_max_word,*)n_max
  read(n_group_word,*)n_group
  print *, 'Using ', trim(file_name), ' n_max = ',n_max,' n_group = ',n_group
  else
   print '(a37,$)',' single particle by element file ? '
   read(5,'(a)')file_name
   print *, 'Using ', trim(file_name)
  endif

 open(unit=11, file=file_name)
 read(11,'(a)')string
! print *,string
 n=0
 do while(.true.)
  read(11,'(a)', end=99)string
  n=n+1
 end do
99 continue
close(11)
print *,' number of lines = ', n
allocate(time(1:n))
allocate(pol(1:n))
allocate(time_zero(1:n/(n_group-1)))

 open(unit=11, file=file_name)
 read(11,'(a)')string
 n=0
 do while(.true.)
  read(11,'(a)', end=199)string
  i=0
  do j=1,15
   call string_trim(string(i+1:),string,i)
   word(j) = string(1:i)
  end do
  n=n+1
  read(word(9),*)time(n)
  read(word(15),*)pol(n)
 end do
199 continue
close(11)

nzero = 0
if(n_max == 0) n_max = n
do i=n_group/2+1,n_max- n_group/2   !for each group of n_group points, if half of points are > 0 and half are < 0 then count 
!  print '(11es12.4)',pol(i-n_group/2:i)
!  print '(11es12.4)',pol(i+1:n_group/2)
  if((all(pol(i-n_group/2:i) > 0) .and. all(pol(i+1:i+n_group/2) < 0)) .or. &
       (all(pol(i-n_group/2:i) < 0) .and. all(pol(i+1:i+n_group/2) > 0)) )then
    nzero = nzero + 1
    time_zero(nzero) = time(i)
    write(13,'(i10,1x,es16.8)')nzero,time_zero(nzero)
  endif
end do
freq = (nzero - 1)/(time_zero(nzero)-time_zero(1))/2.
write (13,'(a,es16.8,a,es16.8)')' f = ',freq ,' omega = ', twopi*freq
write (6,'(a,es16.8,a,es16.8)')' f = ',freq ,' omega = ', twopi*freq, ' delta_omega/omega = ', (freq - freq_0)/freq_0
do i=1,n_max
 write(12,'(i10,1x,2es16.8)')i,time(i),pol(i)
end do
end