subroutine gridsampling(sampfunc,ihowsamp,nparams,
     &  msect0,bestguess,yatguess,guessconfid0,bmax,bmin)
      implicit none
!
      integer nparams,msect0,ihowsamp
      double precision bestguess(nparams),guessconfid0,
     &  bmax(nparams),bmin(nparams),params(nparams,msect0+1),
     &  guessconfid,yatguess,beta(nparams),fatbeta
      integer i,nright,nleft,j,k,n,msect,m
      double precision accum,x1,delta,eps

      logical resetran2
      common /ran2reset/resetran2

      parameter(eps=1.0d-9)
      external sampfunc
!
      resetran2=.true.
!
      msect=msect0
!      call sampfunc(nparams,bestguess,yatguess)
!assume calculation already done for initial guess
      guessconfid=dmax1(1.0d0,guessconfid0)
      j=0
      do i=1,nparams
        if(bestguess(i).lt.bmin(i).or.bestguess(i).gt.
     &    bmax(i))then
          j=1
          if(bestguess(i).lt.bmin(i))then
            bestguess(i)=bmin(i)
          else
            bestguess(i)=bmax(i)
          endif
        endif
      enddo
      if(j.eq.1)call sampfunc(nparams,bestguess,yatguess)
      do i=1,nparams
        if(dabs(bmax(i)-bestguess(i)).lt.eps)then
          nright=0
          nleft=msect+1
        else
          if(dabs(bestguess(i)-bmin(i)).lt.eps)then
            nright=msect+1
            nleft=0
          else
            if(mod(msect,2).eq.0)then
              nright=msect/2+1
              nleft=msect/2+1
            else
              nright=msect/2+1+1
              nleft=msect/2+1
            endif
          endif
        endif
!first divide the right
        if(nright.gt.0)then
          x1=2.0d0*(bmax(i)-bestguess(i))/
     &      (dble(nright)*(guessconfid+1.0d0))
          delta=x1*(guessconfid-1.0d0)/dble(nright-1)
          accum=0.0d0
          do j=1,nright-1
            accum=accum+x1+dble(j-1)*delta
            params(i,j)=accum+bestguess(i)
          enddo
        endif
!then divide the left
        if(nleft.gt.0)then
          x1=2.0d0*(bestguess(i)-bmin(i))/
     &      (dble(nleft)*(guessconfid+1.0d0))
          delta=x1*(guessconfid-1.0d0)/dble(nleft-1)
          accum=0.0d0
          do j=1,nleft-1
            accum=accum+x1+dble(j-1)*delta
            if(nright.eq.0)then
              params(i,j)=bestguess(i)-accum
            else
              params(i,j+nright-1)=bestguess(i)-accum
            endif
          enddo      
        endif
      enddo
      do i=1,nparams
        params(i,msect+1)=bestguess(i)
      enddo
      msect=msect+1
      if(ihowsamp.eq.1)then
!Using random permutation
        call randpermut_dim_samp(msect,nparams,
     &    params(1:nparams,1:msect))
        do i=1,msect
          call sampfunc(nparams,params(1:nparams,i:i),fatbeta)
          if(fatbeta.lt.yatguess)then
            yatguess=fatbeta
            do j=1,nparams
              bestguess(j)=params(j,i)
            enddo
          endif
        enddo
      endif
      if(ihowsamp.eq.2)then
!uniform grid sampling          
        do i=1,msect**nparams
          do j=1,nparams
!the size of the larger repeated unit is msect**(nparams-j+1)
            k=i/(msect**(nparams-j+1))
            n=mod(i,(msect**(nparams-j+1)))
            if(n.eq.0)k=k-1
!k is the number of repeated units before i (not include the unit i is in)
            k=i-k*(msect**(nparams-j+1))
!now k is the position in the larger repeated unit
!
!the size of the smaller repeated unit is (msect**(nparams-j+1))/msect
            m=(msect**(nparams-j+1))/msect
            n=k/m
            if(mod(k,m).ne.0)n=n+1
            beta(j)=params(j,n)
          enddo
          call sampfunc(nparams,beta,fatbeta)
          if(fatbeta.lt.yatguess)then
            yatguess=fatbeta
            do j=1,nparams
              bestguess(j)=beta(j)
            enddo
          endif
        enddo
      endif
      return
      end subroutine gridsampling