********************************************************************************
*                                                                              *
*  2DHF version 1-2003                                                         *
*  Copyright (C) 1996  Jacek Kobus, Leif Laaksonen, Dage Sundholm              *
*                                                                              *
*  This software may be used and distributed according to the terms            *
*  of the GNU General Public License, see README and COPYING.                  *
*                                                                              *
********************************************************************************
c ### nmucalc ###
c
c     Calculates nmu for given nni and practical infinity values 
c     (a single grid case is assumed) 
c
      function nmucalc()
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      character*8 sigma,pi,delta,phi
      character*8 spac,ger,uger
      include 'commons8.inc'

c     mxnmu --  maximum no. of grid points in mu variable
c               mxmnu should be of the form 6N+1
c     hni - step in ni variable
c     hmu - step in mu variable
      
c     determine step size in nmu variable so as hmu is approximately
c     equal to hni

c     adjust, if necessary, the number of grid points in nu variable
      k1=5
      k2=6
 12   nord=(nni-1)/k2
      if (k2*nord.ne.nni-1) then
         nni=nni-1
         goto 12
      endif

c     estimate nmu assuming that hmu=hni

      hni=pii/dble(nni-1)

      xmi0=2.d0*rgrid(1)/r
      xmi0=log(xmi0+sqrt(xmi0*xmi0-1.d0))
c      hmu(1)=xmi0/dble(mxnmu-1)
      mxnmu=xmi0/hni+1
      nmu(1)=mxnmu

c     adjust, if necessary, the number of grid points in mu variable
      ig=1
 20   continue
      ichk1=0
      ichk2=0
 11   nord=(nmu(ig)-1)/k2
      if (k2*nord.ne.nmu(ig)-1) then
         nmu(ig)=nmu(ig)-1
         ichk2=ichk2+1
         goto 11
      endif
      if (ichk1.ne.0.or.ichk2.ne.0) goto 20

      mxnmu=nmu(1)

c     calculate hmu(1) for the adjusted value of mxnmu
      xmi0=2.d0*rgrid(1)/r
      xmi0=log(xmi0+sqrt(xmi0*xmi0-1.d0))
      hmu(1)=xmi0/dble(mxnmu-1)

      nmucalc=mxnmu
      
      return
      end