subroutine funcs(x,coef,yfit,dydcoef)
      use precision_def
      use nrtype
      implicit none
      real(rp), dimension(:), intent(in) :: x, coef
      real(rp), dimension(:), intent(out) :: yfit
      real(rp), dimension(:,:), intent(out) :: dydcoef
      integer i
      
      yfit(:)=0
      do i=1,size(coef)
        yfit(:)  = yfit(:) + coef(i)*x(:)**(i-1)
        dydcoef(:,i) = x(:)**(i-1)
!      dydcoef(:,2) = x(:)
!      dydcoef(:,3) = x(:)**2
!      dydcoef(:,4) = x(:)**3
      end do
    end subroutine funcs
    
 subroutine fit_2D_array(x, y, kx, order, number)
  use precision_def
  use sim_utils
  use compile_interface, dummy => fit_2D_array
  use nr
  implicit none
  real(rp) chisq, alambda
  real(rp) old_chisq
  real(rp), allocatable :: x(:), y(:), sig(:), coef(:)
  real(rp), allocatable :: covar(:,:), alpha(:,:)
  real(rp), allocatable :: coef_err(:)
  logical, allocatable :: maskcoef(:)
  integer kx,i
  integer lun
  integer number
  integer order
  character*7 cnumber
  character*1 cn(10)/'0','1','2','3','4','5','6','7','8','9'/

       allocate(coef(1:order),maskcoef(1:order), coef_err(1:order))
       allocate(covar(1:order,1:order),alpha(1:order,1:order))
       allocate(sig(1:kx))
       coef(1:order) = 0
       coef(2) =0.18
       sig(1:kx)=1.
       maskcoef(:)=.true.
       alambda = -1.
       old_chisq=1000000.
       chisq=0.
       print '(/,a)',' FIT_2D_ARRAY: no binning. fit to scatter plot'
       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:order)
       print '(a,es12.4,a,i10)',' chi2 ', chisq, ' number data points = ',kx
       print '(a)',' covariance matrix'
       do i=1,size(covar,1)
          print '(4es12.4)',covar(i,1:order)
       end do
       forall(i=1:order)coef_err(i)=(sqrt(covar(i,i)/(kx-order)*chisq))
       print '(a,i1,a,4es12.4)',' Errors = (covar(i,i)/(number of data points-',order,')*chisq)**.5 =',(sqrt(covar(i,i)/(kx-order)*chisq),i=1,order)
       write(cnumber,'(i1)')number
       if(number >=10)write(cnumber,'(i2)')number    
       lun=lunget()
       open(unit=lun, file='2D_hist_fit_coef_err_'//trim(cnumber)//'.dat')
       write(lun,'(a,i10)')'number = ',number
       write(lun,'(4(a6,es12.4,a1))')( 'coef'//cn(i)//'=',coef(i),';',i=1,order)
       write(lun,'(4(a10,es12.4,a1))')('coef'//cn(i)//'_err=',coef_err(i),';', i=1,order)
       write(lun,'(a,es12.4)')'chisq=',chisq
       close(unit=lun)
       print '(a,a)','write ', '2D_hist_fit_coef_err_'//trim(cnumber)//'.dat'

  return
end subroutine fit_2D_array