program energy_time_to_2D_hist
    use precision_def
    use sim_utils
    use compile_interface
    use nr
  implicit none
  
  real(rp), allocatable :: column(:)
  real(rp), allocatable :: twod_sum(:,:), onedsum(:)
  real(rp), allocatable :: average(:),rms(:)

  real(rp) chisq, alambda
  real(rp) old_chisq
  real(rp), allocatable :: x(:), y(:), sig(:), coef(:)
  real(rp), allocatable :: covar(:,:), alpha(:,:)
  logical, allocatable :: maskcoef(:)
  integer kx,i
  
  real(rp) xmin, xmax, ymin, ymax
  real(rp) ix_0, iy_0
  real(rp) xscale,yscale
  real(rp) xconvert,yconvert, stdev
  real(rp) xoffset, yoffset

  real(rp), allocatable :: xcolumn(:), ycolumn(:)
  real(rp)  average_energy, rms_energy, stdev_energy
  real(rp) eta
  real(rp) betax, alphax, gammax

  character*300 file
  character*600 long_string
  character*10 cxmin,cxmax,cix_bins,cymin,cymax,ciy_bins
  character *25 selected_parameters(6)/'x,xp','y,yp','time-momentum','momentum-x_inj','momentum-xp_inj','momentum-time'/
  character*7 cnumber
  character*100 string
  logical itexists
  logical special/.false./
  integer nargs
  integer n,ix,iy, ix_bins, iy_bins,iz
  integer lun,lun1
  integer reason
  integer iyy, ixx
  integer ix_min, ix_max, iy_min, iy_max
  integer hist_d/2/
  integer, allocatable :: nonzero_bins(:)
  integer number
  integer at_limit(3)
  integer nmax/20000/  !max number of points to read
  integer order

  print '(a,$)',' Select parameter pair: 1-x,xp; 2-y,yp; 3-t,p; 4-p,x_inj, 5-p,xp_inj, 6-p,t '
  read(5,*) number
  
  select case (number)
  case(1)
     ix=4;iy=5  !x,xp phase space
     xmin=-0.05; xmax=0.05
     ymin = -0.01; ymax=0.01
  case(2)
     ix=6;iy=7  !y,yp phase space
     xmin=-0.05; xmax=0.05
     ymin = -0.01; ymax=0.01
  case(3)
     ix=3;iy=2  !time,momentum phase space
     xmin=-150e-9;xmax=150e-9
     ymin=-0.004; ymax= 0.004
  case(6)
     ix=2;iy=3  !momentum, time phase space
     ymin=-150e-9;ymax=100.e-9
     xmin=-0.006; xmax= 0.006
  case(4)
     ix=2;iy=4  !momentum, initial position
     xmin=-0.006; xmax=0.006
     ymin = 0.068; ymax = 0.086
  case(5)
     ix=2;iy=5  !momentum, initial angle
     xmin=-0.006; xmax=0.006
     ymin = -0.01; ymax=0.01
  end select
  
print *, 'Using: ',selected_parameters(number)
print '(2(a,i2))','ix=',ix,' iy=',iy
print '(a,4(a7,es12.4))','Default ranges: ' ,' xmin=',xmin,' xmax= ',xmax,' ymin=',ymin,' ymax = ',ymax  
print '(a,$)', ' type alternate values or carriage return to continue: '
read(5,'(a)')string
if(string(1:4) == '    ')then
   print '(a,4(a7,es12.4))','use: ' ,' xmin=',xmin,' xmax= ',xmax,' ymin=',ymin,' ymax = ',ymax
else
   read(string,*)xmin,xmax,ymin,ymax
   print '(a,4(a7,es12.4))','use: ' ,' xmin=',xmin,' xmax= ',xmax,' ymin=',ymin,' ymax = ',ymax
endif

print *, 'string=', string
file = 'all_Energy_vs_time_0.dat'
! to histogram vs time and momentum

     ix_bins=100
     iy_bins = 100

     !to histograms vs xinj and momentum
     
    print '(a)',file

 ix_0 = ((xmax+xmin)/(xmax-xmin)/2)*ix_bins
 ix_0 = ((xmax+xmin)/2)
 xscale = (xmax-xmin)/ix_bins
  iy_0 = ((ymax+ymin)/(ymax-ymin)/2)*iy_bins
  iy_0 = (ymax+ymin)/2
 yscale=(ymax-ymin)/iy_bins
  allocate(twod_sum(-ix_bins/2-1:ix_bins/2+1, -iy_bins/2-1:iy_bins/2+2))
  allocate(nonzero_bins(-ix_bins/2-1:ix_bins/2+1))
  allocate(average(-ix_bins/2-1:ix_bins/2+1))
  allocate(rms(-ix_bins/2-1:ix_bins/2+1))
  allocate(x(1:ix_bins),y(1:ix_bins),sig(1:ix_bins))
  twod_sum(:,:)=0
  nonzero_bins(:)=0
  average(:)=0
  rms(:)=0
  x(:)=0
  y(:)=0
  sig(:)=0
  allocate(xcolumn(1:nmax), ycolumn(1:nmax))
  xcolumn(:)=0.
  ycolumn(:)=0.
 
  inquire (file=file, exist = itexists)
  if(.not. itexists)then
     print *,' File not found '
     stop
  endif
  
  print '(2a)', 'open ',file
  open(unit=5,file=file)
  n=0
  at_limit(1:3)=0
    do while(.true.)
       read(5,'(a)',IOSTAT=reason)long_string
       if(reason < 0) exit
       if(index(long_string,'#')/=0 .or. index(long_string,'!')/=0 .or. index(long_string,'muon')/=0)cycle
       call get_column_data(long_string, column)
       n=n+1
       if(n>nmax)then
         print '(a,i10,a)','Number of data points exceeds ',nmax,' Increase nmax'
         stop
       endif
       call create_2D_histogram_array(column,ix,iy,iz, xmin, xmax, ix_bins,ymin,ymax,iy_bins, ix_0, iy_0, twod_sum,number, at_limit)
!          if(special)call create_2D_special_histogram_array(column,ix,iy, xmin, xmax, ix_bins,ymin,ymax,iy_bins, ix_0, iy_0, sum )       
       xcolumn(n) = column(ix)
       ycolumn(n) = column(iy)   
    end do
    close(unit=5)
    order =4
    call fit_2D_array(xcolumn, ycolumn, n, order,number)
    average_energy = sum(xcolumn)/float(n)
    rms_energy = sum(xcolumn(:)**2)/float(n)
    stdev_energy = sqrt(rms_energy - average_energy**2)
    print '(a,es12.4,a,es12.4)','average_energy=',average_energy, '  stdev=', stdev_energy  
    eta = 7.112/(1-0.1058)
    print '(a,es12.4,a,es12.4)','average_r=',average_energy*eta, '  stdev_r=', stdev_energy*eta  
    write(cnumber,'(i1)')number    
    lun=19
    lun1=lunget()
    open(unit=lun, file='2D_hist_'//trim(cnumber)//'.dat')
    print '(a,a)','write: ','2D_hist_'//trim(cnumber)//'.dat'

    write(lun,'(a1,a,a1,a,a1)')'#',' Using= ','"',trim(selected_parameters(number)),'"'
    write(lun,'(a1,a,i5)')'#',' type = ',number
       write(lun,'(a1,2(a,1x,es12.4))')'#', ' ix_0 = ', ix_0,'; iy_0 = ', iy_0
       write(lun,'(a1, 4(a,es12.4))')'#', ' xmax = ',xmax, '; xmin = ',xmin, '; ymax = ',ymax, ';  ymin = ',ymin
       write(lun,'(a1,a,1x,i10,1x,a,1x,i10)')'#', ' ix_bins = ', ix_bins, '; iy_bins = ', iy_bins
       write(lun,'(a1,a,es12.4,a,es12.4)')'#',' average_energy = ',average_energy, '; stdev_energy = ',stdev_energy

       open(unit=lun1, file='bin_average_rms_'//trim(cnumber)//'.dat')
       print '(a,a)','write: ','bin_average_rms_'//trim(cnumber)//'.dat'
       write(lun1,'(a1,a,a1,a,a1)')'#',' Using: ','"',trim(selected_parameters(number)),'"'
       write(lun1,'(a1,a,i5)')'#',' type = ',number       
       write(lun1,'(1a, 5a12)')'#',' bin','average','err_average','stdev','nonzero_bins'

    xconvert = (xmax-xmin)/ix_bins
    xoffset = (xmax+xmin)/2.
    yconvert = (ymax-ymin)/iy_bins
    yoffset =  (ymax+ymin)/2.
    kx=0
       do ixx=-ix_bins/2, ix_bins/2
          do iyy= -iy_bins/2, iy_bins/2
             write(lun, '(2i10,es16.8)')ixx,iyy,twod_sum(ixx,iyy)
!             write(19, '(2i10,es16.8)')ixx,iyy,twod_sum(ixx,iyy)
             average(ixx) = average(ixx) + twod_sum(ixx,iyy)*iyy
             rms(ixx) = rms(ixx) + twod_sum(ixx,iyy)*iyy**2
             if(twod_sum(ixx,iyy)/= 0)nonzero_bins(ixx) = nonzero_bins(ixx) + twod_sum(ixx,iyy) 
          end do
          if(nonzero_bins(ixx)/= 0)then
             kx = kx+1
             average(ixx) = average(ixx)/nonzero_bins(ixx)
             rms(ixx) = rms(ixx)/nonzero_bins(ixx)
             stdev=((rms(ixx)-average(ixx)**2)/nonzero_bins(ixx))**0.5
             x(kx) = ixx*xconvert
             y(kx) = average(ixx)*yconvert
             sig(kx) = 1./nonzero_bins(ixx)**.5
!             write(lun1,'(5es12.4)')ixx*xconvert,average(ixx)*yconvert, average(ixx)/sqrt(nonzero_bins(ixx))*yconvert, stdev*yconvert, nonzero_bins(ixx)
             write(lun1,'(4es12.4,i12)')ixx*xconvert+xoffset,average(ixx)*yconvert +yoffset, average(ix)/(nonzero_bins(ixx)**.5)*yconvert, stdev*yconvert, nonzero_bins(ixx)             



          endif
          
       end do

       print '(/,a)',' Fit binned array: fit to histogram '       
       allocate(coef(1:4),maskcoef(1:4))
       allocate(covar(1:4,1:4),alpha(1:4,1:4))
       coef(1:4) = 0
       coef(2) =0.18
       maskcoef(:)=.true.
       alambda = -1.
       old_chisq=100.
       chisq=0.
       print '(8a12)','iteration','chisq','old_chisq','alambda','coef1','coef2','coef3','coef4'
       do i=1,10
          call mrqmin(x(1:kx),y(1:kx), sig(1:kx),  coef,maskcoef,covar,alpha,chisq,funcs,alambda)
          print '(i10,2x,7es12.4)',i,chisq,old_chisq,alambda,coef
          if(alambda == 0)exit
          if(chisq >0.95*old_chisq)alambda=0
          old_chisq = chisq
       end do

       print '(a,4es12.4)', 'coefficients: ', coef(1:4)
       print '(a,es12.4,a,i10)',' chi2 ', chisq, ' number data points = ',kx
       print '(a)',' covariance matrix'
       do i=1,4
          print '(4es12.4)',covar(i,1:4)
       end do
       print '(a,4es12.4)',' Errors = (covar(i,i)/(number of data points-4)*chisq)**.5 =',(sqrt(covar(i,i)/(kx-4)*chisq),i=1,4)
       
       print '(a,i10)',' number of times at x limit = ',at_limit(1)
       print '(a,i10)',' number of times at y limit = ',at_limit(2)
       
       close(unit=lun)
       close(unit=lun1)
       lun=lunget()
       open(unit=lun, file='2D_hist_input_'//trim(cnumber)//'.dat')
       write(lun,'(a,i10)')'number = ',number
       write(lun,'(a9,es12.4,a11,es12.4)')'xmin =',xmin,'; xmax =',xmax
       write(lun,'(a9,es12.4,a11,es12.4)')'ymin =',ymin,'; ymax =',ymax
       write(lun,'(a9,i12,a11,i12)')'ix_bins =',ix_bins,'; iy_bins =',iy_bins
       write(lun,'(4(a,es12.4))') 'coef1=',coef(1),'; coef2=',coef(2),'; coef3=',coef(3),'; coef4=',coef(4)
       write(lun,'(a,es12.4)')'chisq=',chisq
       close(unit=lun)
       print '(a,es12.4,a,es12.4)',' xscale =', xscale,' yscale =', yscale
       print '(a)',' Use plotting_scripts/2D_hist_test.gnu or 2D_hist_test_quad.gnu to plot.  gnuplot -e "number = ?" $ps/2D_hist_test_quad.gnu'
       
  end program energy_time_to_2D_hist
  
    subroutine create_2D_histogram_array(column,ix,iy,iz,xmin, xmax, ix_bins,ymin,ymax,iy_bins,ix_0, iy_0, sum, number, at_limit, xcolumn, ycolumn)
      use precision_def
      implicit none
      real(rp), allocatable:: column(:)
      real(rp), allocatable :: sum(:,:)
      real(rp) xmin, xmax, ymin, ymax
      real(rp) ix_0, iy_0
      real(rp) size/0.5/
      real(rp) xvalue, yvalue
      real(rp), optional :: xcolumn, ycolumn
      integer ix,iy, ixx, iyy, iz
      integer ix_bins, iy_bins
      integer number
      integer at_limit(3)

      xvalue = column(ix)
      yvalue = column(iy)
      
      if(xvalue< xmin)xvalue = xmin
      if(xvalue> xmax)xvalue = xmax
      if(yvalue< ymin)yvalue = ymin
      if(yvalue> ymax)yvalue = ymax
      
      ixx = (xvalue-ix_0)/(xmax-xmin) * ix_bins + sign(size,xvalue-ix_0)
      iyy = (yvalue-iy_0)/(ymax-ymin) * iy_bins + sign(size,yvalue-iy_0)

         if(xvalue ==xmin .or. xvalue == xmax)at_limit(1) = at_limit(1)+1
         if(yvalue == ymin .or. yvalue == ymax)at_limit(2) = at_limit(2)+1

      
      sum(ixx,iyy) = sum(ixx,iyy)+1
      return
    end subroutine create_2D_histogram_array
    
! histogram betatron amplitude vs horizontal closed orbit
      subroutine create_2D_special_histogram_array(column,ix,iy,xmin, xmax, ix_bins,ymin,ymax,iy_bins,ix_0, iy_0, sum )
      use precision_def
      implicit none
      real(rp), allocatable:: column(:)
      real(rp), allocatable :: sum(:,:)
      real(rp) xmin, xmax, ymin, ymax
      real(rp) eta/8.045/, beta/7.6969/
      real(rp) fyy
      real(rp) ix_0, iy_0
      real(rp) size/0.5/
      integer ix,iy, ixx, iyy
      integer ix_bins, iy_bins
      !column 24 = x_e
      !column 9 = pz
      !column 4 = x
      !column 5 = x'

      
      ixx = column(24)/(xmax-xmin) * ix_bins  -ix_0
      ixx = ((column(24)-ix_0)/(xmax-xmin)) * ix_bins +sign(size,column(24)-ix_0) 
      fyy = ((column(4)-column(9)*eta)**2/beta+column(5)**2*beta )
      iyy = ((column(4)-column(9)*eta)**2/beta+column(5)**2*beta )/(ymax-ymin) * iy_bins + 0.5 - iy_0
      iyy = (fyy-iy_0) + sign(size,fyy-iy_0)
      sum(ixx,iyy) = sum(ixx,iyy)+1
      return
      end
      
   subroutine get_column_data(long_string, column)    
    use precision_def
    implicit none
    real(rp), allocatable :: column(:)
    character(*) long_string
    character*600 string
    integer i, ix

    if(.not. allocated(column))then
         string = long_string
         i=0
         ix=0
           do while (.true.)
              call string_trim(string(ix+1:), string, ix)
             if(ix == 0)exit
             i=i+1
          end do
          allocate(column(1:i))
       endif
       
         i=0
         ix=0
           do while (.true.)
             call string_trim(long_string(ix+1:), long_string, ix)
             if(ix == 0)exit
             i=i+1
             read(long_string(1:ix),*)column(i)        
          end do
          return
        end subroutine get_column_data
        
    subroutine create_1D_histogram_array(column,ix,xmin, xmax, ix_bins,ix_0, sum, average,rms2 )
      use precision_def
      implicit none
      real(rp), allocatable:: column(:)
      real(rp), allocatable :: sum(:)
      real(rp) xmin, xmax
      real(rp) ix_0
      real(rp) average, x2_avg,rms2
      real(rp) size/0.5/
      integer ix, ixx
      integer ix_bins
      integer n/0/
      if(n==0)then
         average =0
         x2_avg=0
      endif
      
      if(column(ix)< xmin)column(ix) = xmin
      if(column(ix)> xmax)column(ix) = xmax
      
      ixx = (column(ix)-ix_0)/(xmax-xmin) * ix_bins + sign(size,column(ix)-ix_0)
      
      sum(ixx) = sum(ixx)+1
      n=n+1
      average = (average*(n-1)+column(ix))/n
      x2_avg = (x2_avg*(n-1)+ column(ix)**2)/n
      rms2 = x2_avg-average**2
      return
    end subroutine create_1D_histogram_array