program phase_space_to_histogram
    use precision_def
    use sim_utils
    use compile_interface
    use histo_params
  implicit none
  
  real(rp), allocatable :: column(:)
  real(rp), allocatable :: twod_sum(:,:), onedsum(:)
  real(rp) xmin, xmax, ymin, ymax
  real(rp) ix_0, iy_0
  real(rp) average, rms2
  real(rp), allocatable :: xsave(:), ysave(:)
  real(rp) xcolumn, ycolumn
  character*300 file, file0
  character*600 long_string
  character*10 cxmin,cxmax,cix_bins,cymin,cymax,ciy_bins, cnumber
  character*300  dated_directories(300) 

  logical itexists
  logical special/.false./
  logical multiple_directories
  logical first_call/.true./
  logical write_column_file/.false./

  integer nargs
  integer n,ix,iy, ix_bins, iy_bins, iz
  integer lun
  integer reason
  integer iyy, ixx
  integer ix_min, ix_max, iy_min, iy_max
  integer hist_d/2/
  integer number
  integer ios
  integer at_limit(3)/0,0,0/, at_limit_0(3)
  integer nmax/20000000/  !max number of points to read
  integer order
  integer number_directories, dir_count
  integer i
  integer muon
  integer, save :: lun95

  namelist/phase_space_to_hist/file0,number,xmin,xmax,ix_bins,ymin,ymax,iy_bins, multiple_directories, write_column_file
  
  number =7
  file0 = 'all_EndOfTracking_phase_space.dat'  !as an example
  
! nargs = command_argument_count()
! if (nargs >= 1)then
!    call get_command_argument(1, file)
!    call get_command_argument(2, cnumber)
!    read(cnumber,*)number
!    call get_command_argument(3, cxmin)
!    read(cxmin,*)xmin
!    call get_command_argument(4, cxmax)
!    read(cxmax,*)xmax
!    call get_command_argument(5, cix_bins)
!    read(cix_bins,*)ix_bins
!    if(nargs > 5)then
!       call get_command_argument(6, cymin)
!       read(cymin,*)ymin
!       call get_command_argument(7, cymax)
!       read(cymax,*)ymax
!       call get_command_argument(8, ciy_bins)
!       read(ciy_bins,*)iy_bins
!    endif

    multiple_directories=.false.
    OPEN (UNIT=5, FILE='phase_space_to_hist.dat', STATUS='old', IOSTAT=ios)
    READ (5, NML=phase_space_to_hist, IOSTAT=ios)
    print *, 'ios=', ios
    rewind(unit=5)
    CLOSE(5)
    WRITE(6,NML=phase_space_to_hist)  
!      muon               Jx              Jy               x              xp               y              yp          vec(5)              pz               t               s             s_x             s_y             s_z  
  select case(number)
  case(1)  !x,xp phase space
     if(index(file0,'muon_end.dat') /=0)then
        ix=3;iy=4
     else
        ix=4;iy=5
   endif
  case(2)
     ix=6;iy=7  !y,yp phase space
     if(index(file0,'muon_end.dat') /=0)then
        ix=5;iy=6
     endif  
  case(3)
     ix=3;iy=2  !time,momentum phase space
     if(index(file0,'muon_end.dat') /=0)print *,'No time data in: ',file0
  case(4)
     if(index(file0,'particles_M4M5End_400_mod.txt')/=0)then
       ix=-1;iy=23 !x,atan(spin_x/spin_z), from  particles_M4M5End_400_mod.txt
       iz=21 
     endif

     if(index(file0,'muon_end.dat') /=0)then
        ix=3;iy=12;iz=10
     endif
     if(index(file0,'phase_space')/=0)then
        iy=14;iz=12
     endif

   case(5)  !spin vs x'
     if(index(file0,'particles_M4M5End_400_mod.txt')/=0)then
       ix=-4;iy=23 !x',spin_x, from  particles_M4M5End_400_mod.txt
       iz=21  !y_spin
     endif
     if(index(file0,'muon_end.dat') /=0)then
        ix=4;iy=12;iz=10
     endif
     if(index(file0,'phase_space')/=0)then
        ix=5;iy=14;iz=12
     endif
  case(6)   !spin angle in x-p plane vs delta
! particles_M4M5End_400_mod.txt
!     #x y z Px Py Pz t PDGid EventID TrackID ParentID Weight Bx By Bz Ex Ey Ez ProperTime PathLength PolX PolY PolZ InitX InitY InitZ InitT InitKE
     if(index(file0,'particles_M4M5End_400_mod.txt')/=0)then
       ix=-1;iy=21 !x',spin_x, from  particles_M4M5End_400_mod.txt
       iz=23  !y_spin
     endif
     if(index(file0,'muon_end.dat') /=0)then
        ix=8;iy=12;iz=10
     endif
     if(index(file0,'phase_space')/=0)then
        ix=9;iy=12;iz=13
     endif
     if(index(file0,'EndofM5_Valetov_withInit')/=0)then
        ix=9;iy=12;iz=13
     endif
     
  case(7)   !acos((spin'p)/|s||p| angle vs x - spin angle in x-p plane
     if(index(file0,'particles_M4M5End_400_mod.txt')/=0)then
        ix=-1;iy=21 !x, spin_z, from
        iz=22  !x_spin
     endif
     if(index(file0,'EndofM5_Valetov_withInit')/=0)then
        ix=4;iy=12;iz=13
     endif

     if(index(file0,'muon_end.dat') /=0)then
        ix=8;iy=10;iz=11
     endif
     if(index(file0,'phase_space')/=0)then
        ix=4;iy=12;iz=13
     endif
  case(8)   !acos((spin'p)/|s||p| angle in x-p plane vs x'
     if(index(file0,'particles_M4M5End_400_mod.txt')/=0)then
        ix=-4;iy=21 !x, spin_z, from
        iz=22  !x_spin
     endif
     if(index(file0,'muon_end.dat') /=0)then
        ix=9;iy=10;iz=11
     endif
     if(index(file0,'phase_space')/=0)then
        ix=5;iy=12;iz=13
     endif
     if(index(file0,'EndofM5_Valetov_withInit')/=0)then
        ix=5;iy=12;iz=13
     endif
  case(9)   !acos((spin'p)/|s||p| angle vs delta xz plane
     if(index(file0,'particles_M4M5End_400_mod.txt')/=0)then
        ix=-4;iy=21 !x, spin_z, from
        iz=22  !x_spin
     endif
     if(index(file0,'muon_end.dat') /=0)then
        ix=9;iy=10;iz=11
     endif
     if(index(file0,'phase_space')/=0)then
        ix=5;iy=12;iz=13
     endif
     if(index(file0,'EndofM5_Valetov_withInit')/=0)then
        ix=5;iy=12;iz=13
     endif
  case(10)   !acos((spin'p)/|s||p| angle vs invariant amplitude
     if(index(file0,'EndofM5_Valetov_withInit')/=0)then
          betax=5.2394950385E+01     !twiss parameters of raw distribution    
          alphax=1.4788683882E+01    
          gammax=4.1932508641E+00
          etax = -0.0481589
          etapx= 1.5424866889E-02
          betay=4.56
          alphay=3.10
          gammay=2.33
          etay = 2.8104057071E-02  
          etapy =  -7.0962834848E-03
        ix=5;iy=12;iz=13
        ytotal=0.
        entry=0.
     endif
   case(11)  !delta vs spin angle in x-p plane
     if(index(file0,'EndofM5_Valetov_withInit')/=0 .or. index(file0,'phase_space')/=0)then
        ix=9;iy=12;iz=13
     endif
   case(12)  !spin angle in x-p plane vs y'
     if(index(file0,'EndofM5_Valetov_withInit')/=0 .or. index(file0,'phase_space')/=0)then
        ix=7;iy=12;iz=13
     endif
   case(13)  !spin angle in x-p plane vs y
     if(index(file0,'EndofM5_Valetov_withInit')/=0 .or. index(file0,'phase_space')/=0)then
        ix=6;iy=12;iz=13
     endif
   case(14)  !sy vs yp
     if(index(file0,'EndofM5_Valetov_withInit')/=0 .or. index(file0,'phase_space')/=0)then
        ix=7;iy=13;iz=13
     endif
   case(15)  !sy vs p
     if(index(file0,'EndofM5_Valetov_withInit')/=0 .or. index(file0,'phase_space')/=0)then
        ix=6;iy=13;iz=13
     endif
   case(16)  !sx vs xp
     if(index(file0,'EndofM5_Valetov_withInit')/=0 .or. index(file0,'phase_space')/=0)then
        ix=5;iy=12;iz=13
     endif
   case(17)  !sx vs x
     if(index(file0,'EndofM5_Valetov_withInit')/=0 .or. index(file0,'phase_space')/=0)then
        ix=4;iy=12;iz=13
     endif
   case(18)  !phi vs delta (phi read from column 15)
     if(index(file0,'EndofM5_Valetov_withInit')/=0 .or. index(file0,'phase_space')/=0)then
        ix=9;iy=15
     endif

   end select


 ix_0 = ((xmax+xmin)/(xmax-xmin)/2)*ix_bins
 ix_0 = ((xmax+xmin)/2)
 if(ymax-ymin >0)then
  iy_0 = ((ymax+ymin)/(ymax-ymin)/2)*iy_bins
  iy_0 = (ymax+ymin)/2
  allocate(twod_sum(-ix_bins/2-1:ix_bins/2+1, -iy_bins/2-1:iy_bins/2+2))
  twod_sum(:,:)=0
else
   allocate(onedsum(-ix_bins/2-1:ix_bins/2+1))
   onedsum(:)=0
endif

  allocate(xsave(1:nmax), ysave(1:nmax))  

! these next few lines allow looping over all dated directories

number_directories=1
if(multiple_directories)call get_all_dated_directories(number_directories, dated_directories)
dir_count=1
n=0
do while(dir_count <= number_directories)  
if(number_directories > 1)then
   file = trim(dated_directories(dir_count))//'/'//trim(file0)
  else
   file = file0
endif
!

  inquire (file=file, exist = itexists)
  if(.not. itexists)then
     print *,' File not found '
     stop
  endif
  
  print '(2a)', 'open ',file
  open(unit=5,file=file)
  dir_count=dir_count+1
  

  
    do while(.true.)
       read(5,'(a)',IOSTAT=reason)long_string
!       print '(a12,a100)','long_string=',trim(long_string)
       if(reason < 0) exit
       if(index(long_string,'#')/=0 .or. index(long_string,'!')/=0 .or. long_string(1:10)  ==  '          ' .or. index(long_string,'muon')/=0)then
          if(first_call .and. write_column_file)then
             lun95=lunget()
             open(lun95,file=trim(file0)//'_trimmed')
             write(lun95,'(a2,a)')'# ',file0
             write(lun95,'(a10,14a16)')'# muon','jx,','jy','x','px/p','y','py/p','vec(5)','dp/p','t','s','sx','sy','sz','phi'
             first_call = .false.
          endif
         cycle
       endif
!       write(95,'(a)')trim(long_string(1:200))
       call get_column_data(long_string, column)
       muon = column(1)

       if(hist_d == 2)then
          if(n>nmax)then
            print '(a,i10,a)','Number of data points exceeds ',nmax,' Increase nmax'
            stop
          endif
          at_limit_0(1:3) = at_limit(1:3) 
          if(.not. special)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, xcolumn, ycolumn)

          if(write_column_file .and. ycolumn < 98.9 .and. xcolumn > -98.9 )write(lun95,'(i11,14es16.8)')muon,column(2:3),column(4:8),column(9:11), column(12:14),ycolumn

          if(xmin < xcolumn .and. xcolumn < xmax .and. ymin < ycolumn .and. ycolumn < ymax)then
             
             if(any(at_limit(1:3)-at_limit_0(1:3)>0))cycle
            n=n+1
            xsave(n) =xcolumn
            ysave(n)=ycolumn
            write(97,'(i10,2es12.4)')n, xcolumn, ycolumn
           endif
          if(special)call create_2D_special_histogram_array(column,ix,iy, xmin, xmax, ix_bins,ymin,ymax,iy_bins, ix_0, iy_0, twod_sum )
       else
          call create_1D_histogram_array(column,ix, xmin, xmax, ix_bins, ix_0, onedsum, average, rms2 )
       endif
       
          
    end do !read through file to end
  end do !read through list of dated directories


    order = 2
    do i=1,n
      write(96,'(i10,2es12.4)')i,xsave(i),ysave(i)
    end do
    call fit_2D_array(xsave, ysave, n,order, number)
!    call fit_2D_array(ysave, xsave, n,order, number)
       lun=lunget()
       write(cnumber,'(i1)')number
       open(unit=lun, file='2D_hist_fit_coef_err_'//trim(cnumber)//'.dat',access="append")
       write(lun,'(2(a,es12.4))')'xmax=',xmax,'; xmin=',xmin
       write(lun,'(2(a,es12.4))')'ymax=',ymax,'; ymin=',ymin
!       write(lun,'(a,es12.4,a,es12.4)')' average x =',sum(xsave(:))/n,' average y = ',sum(ysave(:))/n
       close(lun)


    lun=lunget()
    if(ymax-ymin>0)then
       open(unit=lun, file='2D_hist.dat')
       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
       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)
          end do
       end do
       close(unit=lun)
       print '(a,a,a)','write', ' 2D_hist.dat.',' Plot with $ps/2D_histogram.gnu'
       print '(a,es12.4)',' sum(twod_sum) ',sum(twod_sum)
       call compress_2D(ix_bins,iy_bins,twod_sum, xmin, xmax, ymin, ymax,number,file)
       print '(a,i10)',' number of times at x limit = ',at_limit(1)
       print '(a,i10)',' number of times at y limit = ',at_limit(2)
       print '(a,i10)',' number of times |spin|=0  = ',at_limit(3)       
       if(entry >0)print '(a,es12.4)',' Average phi =', ytotal/entry
       print '(a)',' Plot with   /plotting_scripts/phase_space_to_hist.gnu'
    else
       open(unit=lun, file='1D_hist.dat')
       write(lun,'(a1,(a,1x,es12.4))')'#', ' ix_0 = ', ix_0
       write(lun,'(a1, 2(a,es12.4))')'#', ' xmax = ',xmax, ' xmin = ',xmin
       write(lun,'(a1,a,1x,i10,1x)')'#', ' ix_bins = ', ix_bins
       write(lun,'(a1,a,es12.4,a,es12.4)')'#',' average = ', average,' rms2 = ', rms2
       do ixx=-ix_bins/2, ix_bins/2
             write(lun, '(i10,es16.8)')ixx,onedsum(ixx)
       end do
       close(unit=lun)
    endif
    
       
  end program phase_space_to_histogram
  
    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 sim_utils
      use precision_def
      use histo_params
      implicit none
      real(rp), allocatable:: column(:)
      real(rp), allocatable :: sum(:,:)
      real(rp) xmin, xmax, ymin, ymax
      real(rp) ix_0, iy_0
      real(rp) yvalue
      real(rp) xvalue
      real(rp) px,py,pz,ptot,xp
      real(rp) deltap
      real(rp) spinx,spiny,spinz,stot,sdotp_norm
      real(rp) half/0.5/
      real(rp), optional :: xcolumn, ycolumn
      real(rp) ampx2, ampy2
      real(rp) gamma0
      real(rp) xb, xpb, yb, ypb
      real(rp) khat(3), svec(3), pvec(3), svec_p(3), sxp(3)
      real(rp) sdotp
      integer ix,iy, ixx, iyy,iz
      integer ix_bins, iy_bins
      integer number
      integer at_limit(3)

      if(present(xcolumn))then
        xcolumn=-99.
        ycolumn=99.
      endif
! ix<0 means diktys style      
      
      if(ix <0)then
         px=column(4) 
         py=column(5)
         pz=column(6)
      else
         
         px = column(5)
         py = column(7)
         pz = sqrt((1+column(9))**2-px**2-py**2)
      endif
      
      xvalue = column(iabs(ix))
      
      if(iy < 15)then
         spinx = column(iy)
         spiny = column(iy+1)
         spinz = column(Iy+2)

         ptot = sqrt(px**2+py**2+pz**2)
         stot = sqrt(spinx**2+spiny**2+spinz**2)
         if(stot < 0.00001)then
            at_limit(3) = at_limit(3)+1
            return
         endif
         
            
         sdotp_norm= (spinx*px + spiny*py + spinz * pz)/ptot/stot   !standard def

 
         khat = (/0.0_rp,-pz,py/)         
         khat = khat/sqrt(dot_product(khat,khat))  !unit normal to x-p plane
         svec = (/spinx,spiny,spinz/) 
         svec_p = (svec-dot_product(svec,khat)*khat)
         svec_p = svec_p/sqrt(dot_product(svec_p,svec_p))   ! polarization in the x-p plane
         pvec=(/px,py,pz/)
         sdotp = dot_product(svec_p,pvec)
         sdotp_norm = dot_product(svec_p,pvec)/sqrt(dot_product(svec_p,svec_p))/ptot   ! phase of component of spin in x-p plane


         if((number == 7 .or. number == 4)  .and. ix == -1)xvalue=xvalue/1000. !mm to m
         if((number == 8 .or. number == 5) .and. ix == -4)xvalue = px/pz
         if(number == 10)then
            xb =column(4) - etax*column(9)
            xpb =column(5) - etapx*column(9)
            ampx2 = xb**2*gammax + 2*alphax*xb*xpb + xpb**2*betax
            yb =column(6) - etay*column(9)
            ypb =column(7) - etapy*column(9)
            ampx2 = yb**2*gammay + 2*alphay*yb*ypb + ypb**2*betay
            xvalue = sqrt(ampx2+ampy2)
!            xvalue = column(9)
            yvalue = column(iy)
         endif

         if(number == 7 .or. number == 8 .or. number == 6 .or. number == 10 .or. number == 12 .or. number == 13)then
            yvalue = acos(sdotp_norm) ! * sign(1.0_rp, svec_p(1))
            !             call signed_phi(svec, pvec,yvalue)
            call RightOrLeft(pvec,svec,sxp)
            yvalue = atan2(sxp(2),sdotp)
            if(yvalue <0)yvalue = yvalue + twopi
           write(80,'(8es12.4)')xvalue,svec,pvec,yvalue
             deltap=column(9)
         endif

         if(number == 11) then
            xvalue = acos(sdotp_norm)
            yvalue = column(9)
         endif

         
         
         if(number == 10)then
           deltap=column(9)
           gamma0=29.3
           yvalue =  acos(sdotp_norm)-twopi/2. + 4*twopi*(gamma0*(1+deltap)-gamma0) * 1./(gamma0**2-1) + 9.1780E-01 -3.5225E-06          
           ytotal = ytotal + yvalue
           yvalue = abs(yvalue)
           entry = entry + 1
!           yvalue =   twopi/2.- 4*twopi*(gamma0*(1+deltap)) * 1./(gamma0**2-1)
         endif


         if(number == 9)then
            sdotp_norm = (spinx*px+spinz*pz)/((spinx*spinx+spinz*spinz)*(px*px+pz*pz))**.5
            yvalue = acos(sdotp_norm)
            deltap=column(9)
         endif
            
!         if(number == 5 .and. ix == 1) xvalue = px/ptot
         if(number == 6  .and. ix == -1)  xvalue = (ptot-3094.35)/3094.35         
         if(number == 4 .or. number == 5)then
            yvalue=twopi/4.
            if(abs(column(iy))>0.)yvalue = atan(column(iz)/column(iy))   !arctan(spin_x/spin_z)
         endif
      endif  !iy<15
      

         if(number >= 14 .and. number <= 18)yvalue = column(iy)
         
         if(xvalue< xmin)xvalue = xmin
         if(xvalue> xmax)xvalue = xmax
         ixx = (xvalue-ix_0)/(xmax-xmin) * ix_bins + sign(half,xvalue-ix_0)

         if(yvalue< ymin)then
!              print '(a,2es12.4)',' ymin,yvalue =', ymin, yvalue
              yvalue = ymin
         endif
         if(yvalue> ymax)then
!              print '(a,2es12.4)',' ymax,yvalue =', ymax, yvalue
              yvalue = ymax
         endif
         iyy = (yvalue-iy_0)/(ymax-ymin) * iy_bins + sign(half,yvalue-iy_0)

!         if(xvalue < xmin .or. xvalue > xmax)return
!         if(yvalue < ymin .or. yvalue > ymax)return
         
!      print '(2es12.4)',xvalue,yvalue      
         if(xvalue <=xmin .or. xvalue >= xmax)then
           at_limit(1) = at_limit(1)+1
           return
         endif
         if(yvalue <= ymin .or. yvalue >= ymax)then
           at_limit(2) = at_limit(2)+1
           return
         endif

      if(present(xcolumn))then
        xcolumn = xvalue
        ycolumn = yvalue
      endif
      sum(ixx,iyy) = sum(ixx,iyy)+1
      return
      end

! 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) half/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(half,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(half,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
    integer error

    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
!          print '(a,i10)',' number of columns=',i
          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),*, iostat=error)column(i)
 !            print '(2i10,3x,a,3x,es12.4)',i,ix, long_string(1:ix), column(i)
             if(error /= 0 .and. error /= 59 .and. error /= 5010)print *,'read column data read error =',error
          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) half/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(half,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
    
    subroutine compress_2D(ix_bins,iy_bins,twod_sum, xmin,xmax,ymin,ymax, number, file)
    use precision_def
    use sim_utils
    use compile_interface, except_dummy=>compress_2D
    implicit none

  real(rp), allocatable :: twod_sum(:,:)
  real(rp), allocatable :: average(:),rms(:),nonzero_bins(:)
  real(rp) xmin, xmax, ymin, ymax
  real(rp) xconvert, yconvert
  real(rp) xoffset, yoffset
  real(rp) m(2,2),v(2), fit_params(2), minverse(2,2)
  real(rp) stdev
  integer ix_bins, iy_bins, ixx,iyy
  integer lun
  integer number
  integer ind
  character*7 cnumber
  character*300 file
  
    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))

    nonzero_bins(:)=0
    average(:)=0
    rms(:)=0
    m(:,:)=0
    v(:)=0
    lun=lunget()
    write(cnumber,'(i1)')number
    if(number>=10)write(cnumber,'(i2)')number
!    print '(a,i10,a)',' number = ',number,cnumber
    open(unit=lun, file='bin_average_rms_'//trim(cnumber)//'.dat')
    ind = max(index(file,'.dat'), index(file,'.txt'))
    write(lun,'(a8,a)')'file = ','"'//file(1:ind-1)//'"'
    write(lun,'(a4)')'exit'
    write(lun,'(1a, a10,5a12)')'#',' bin','average','err_average','rms ','nonzero_bins','sum'

    xconvert = (xmax-xmin)/ix_bins
    xoffset =  (xmax+xmin)/2.
    yconvert = (ymax-ymin)/iy_bins
    yoffset = (ymax+ymin)/2.

    
           do ixx=-ix_bins/2, ix_bins/2
          do iyy= -iy_bins/2, iy_bins/2
             average(ixx) = average(ixx) + twod_sum(ixx,iyy)*(iyy*yconvert+yoffset)
             rms(ixx) = rms(ixx) + twod_sum(ixx,iyy)*(iyy*yconvert+yoffset)**2
             nonzero_bins(ixx) = nonzero_bins(ixx) + twod_sum(ixx,iyy) 
          end do
          if(nonzero_bins(ixx)/= 0)then
             average(ixx) = average(ixx)/nonzero_bins(ixx)
             rms(ixx) = rms(ixx)/nonzero_bins(ixx)
             stdev = ((rms(ixx) - average(ixx)**2)/nonzero_bins(ixx))**.5

             write(lun,'(6es12.4)')ixx*xconvert+xoffset,average(ixx),  (average(ixx))/sqrt(nonzero_bins(ixx)),  stdev, nonzero_bins(ixx), average(ixx)
! chi2 = sum (a*x_i+b- y_i)^2/N_i
             m(1,1) = m(1,1) + (ixx*xconvert+xoffset)**2/nonzero_bins(ixx)
             m(1,2) = m(1,2) + (ixx*xconvert+xoffset)/nonzero_bins(ixx)
             m(2,1) = m(1,2)
             m(2,2) = m(2,2) + 1./nonzero_bins(ixx)
             v(1) = ixx*xconvert*average(ixx)/nonzero_bins(ixx)
             v(2) = average(ixx)/nonzero_bins(ixx)
          endif
          
       end do
       minverse = -m
       minverse(1,1)=m(2,2) 
       minverse(2,2)=m(1,1)
       fit_params = matmul(minverse,v)
!       chi2 = sum((v(1)*ixx*xconvert+v(2) -iyy*yconvert)**2/nonzero_bins(ixx))
       print '(2(a,es12.4,a,es12.4))','xconvert = ', xconvert,' yconvert = ', yconvert,' xoffset=', xoffset,' yoffset=',yoffset
       print  '(a,2es12.4)', 'fit_params =',fit_params
       
       close(unit=lun)
       print '(2a)',' write ','bin_average_rms_'//trim(cnumber)//'.dat'
end subroutine compress_2D