program bin_data_from_two_columns

use sim_utils
use compile_interface

implicit none

integer nargs, lun
integer lun1
integer numbins
integer reason
integer, allocatable :: count(:)
integer binwidth,i,bin
integer column1, column2
integer last_character
integer out_of_range/0/

character*300 file, long_string
character*10 ccolumn1,ccolumn2, cnumbins, cxmax, cxmin


real(rp), allocatable :: column(:)
real(rp) xmax, xmin
real(rp), allocatable :: binsum(:), p(:)

logical first_call/.true./
logical itexists

nargs = command_argument_count()

call get_command_argument(1, file)
print '(2a)','file = ',trim(file)
call get_command_argument(2, ccolumn1)
read(ccolumn1,*)column1
print '(a,i10)','column 1 = ', column1
call get_command_argument(3, ccolumn2)
read(ccolumn2,*)column2
print '(a,i10)','column 2 = ', column2
call get_command_argument(4, cnumbins)
read(cnumbins,*)numbins
print '(a,i10)','Number of bins = ', numbins
call get_command_argument(5, cxmin)
read(cxmin,*)xmin
call get_command_argument(6, cxmax)
read(cxmax,*)xmax
print '(a,es12.4,a,es12.4)','Min=',xmin,' Max=',xmax

lun=lunget()
print '(a,i10)','lun=',lun
    inquire (file= file , exist = itexists)
    if(itexists)then
       print '(a,a)',trim(file),' exists'
    else
       print '(a,a)',trim(file),' not found'
       stop
    endif

open(unit=lun, file= trim(file))

do while(.true.)
  read(lun,'(a)',IOSTAT=reason)long_string
  last_character = len_trim(long_string)
!  print '(a,i10)','last_character=',last_character
!  if(last_character == 0)cycle
  if(last_character > 0)then
    if(long_string(last_character:last_character) =='i')long_string(last_character:last_character)=' '
    if(long_string(last_character:last_character) =='u')cycle
  endif
!       print '(a12,a300)','long_string=',trim(long_string)
  if(reason < 0) exit
  if(index(long_string,'#')/=0 .or. index(long_string,'!')/=0 .or. long_string(1:10)  ==  '          ' )then
    if(first_call)then
      lun1=lunget()
      open(lun1,file=trim(file)//'_binned')
      first_call = .false.
    endif
    cycle
  endif
!       write(95,'(a)')trim(long_string(1:200))

  call get_column_data(long_string, column)

  call bin_data(column(column1),column(column2), numbins, xmin, xmax, count, binsum, p, out_of_range)

end do

write(lun1,'(a10,2a12,a10)')'bin','x','y','count' 
do i=1,numbins
  if(count(i)>0)then
    binsum(i)=binsum(i)/count(i)
    write(lun1,'(i10,2es12.4,i10)')i, p(i),binsum(i),count(i)
  endif
end do

close(unit=lun1)
close(unit=lun)

print '(a,i10)', 'Number of bins =', numbins
print '(a,es12.4,a,es12.4)',' Max =', xmax,'   Min=',xmin
print '(a,i10)',' Entries = ', sum(count) 
print '(a,a)',' Write to ', trim(file)//'_binned'
print '(i10,a)',out_of_range,' particles out of range'
end program bin_data_from_two_columns




subroutine bin_data(c1, c2, numbins, xmin, xmax, count, binsum, p, out_of_range)

 use precision_def
 implicit none

 integer numbins
 integer, allocatable :: count(:)
 integer bin
 integer out_of_range

 logical first/.true./

 real(rp) xmin, xmax
 real(rp), save :: binwidth
 real(rp), allocatable :: binsum(:), p(:)
 real(rp) c1, c2

 if(first)then
  binwidth = (xmax-xmin)/numbins
  allocate(count(0:numbins), binsum(0:numbins),p(0:numbins))
  count(:)=0
  binsum(:)=0.
  forall(bin=0:numbins) p(bin) = bin*binwidth +xmin
  first=.false.
 endif

 bin = (c1 - xmin)/binwidth +0.5
if(bin < 0 .or. bin > numbins)then
  out_of_range = out_of_range + 1
  return
endif

 count(bin) = count(bin)+1
 binsum(bin) = binsum(bin) + c2

 p(bin) = bin*binwidth + xmin

 return 
end subroutine bin_data