********************************************************************************
*                                                                              *
*  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 ### fermi ###
c     
c     Analyzes Fermi charge distribution.
c
      subroutine fermi(iprint)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
      include 'commons8.inc'

      if (iprint.ne.0) then
         write(*,*) ' '
         write(*,*) '  finite nuclei with the Fermi charge',
     &        ' distribution'
      endif

c     calcualte parameters of the finite nucleus charge distribution
c     (Fremi distribution)

c     ainfcm: bohr radius in cm

      ainfcm=0.529 177 249  d-08
      fmtoau=1.0d-13/ainfcm

      if (z1.ne.0.d0) then

c        set atomic weight for centre A

         atw=z1atmass
         at3=atw**(1.d0/3.d0)
         rrmsfm = 0.836d 00*at3+0.570d 00
         tfm = 2.30d 00

c        change units from fm into bohr

         rrms = rrmsfm*fmtoau
         t = tfm*fmtoau  
         a = t/(4.0d 00*log (3.0d 00))         
         facto1 = rrms**2-(7.0d 00/5.0d 00)*(pii**2)*(a**2)
         c = sqrt (5.0d 00/3.0d 00) * sqrt (facto1) 

         ni=nni
         nc1=0
         xr2=r/2.0
         do mu=1,mxnmu
	    rr=sqrt(vxisq(mu)+vetasq(ni)-1.d0)
            xrr=xr2*rr
            if (abs(xrr-xr2).lt.c) then 
               nc1=nc1+1
               if (idbg(450).ne.0) write(*,*) 'mu,xrr',mu,xrr
            endif
         enddo

         mu=1
         nc2=0
         do ni=nni,nni/2,-1
	    rr=sqrt(vxisq(mu)+vetasq(ni)-1.d0)
            xrr=xr2*rr
            if (abs(xrr-xr2).lt.c) then
               nc2=nc2+1
               if (idbg(450).ne.0) write(*,*) 'ni,xrr',ni,xrr
            endif
         enddo

         if (iprint.ne.0) then
            write(*,*) '    center A: '
            write(*,*) '    c = ',c,' bohr', '  a = ',a,' bohr' 
            write(*,*) '    nu=-1:',nc1, ' points inside radius c'
            write(*,*) '    mu= 1:',nc2, ' points inside radius c'
            write(*,*) ' '
         endif
         elseif (z2.ne.0.d0) then

c           set atomic weight for centre B

         atw=z2atmass
         at3=atw**(1.d0/3.d0)
         rrmsfm = 0.836d 00*at3+0.570d 00
         tfm = 2.30d 00

c        change units from fm into bohr

         rrms = rrmsfm*fmtoau
         t = tfm*fmtoau  
         a = t/(4.0d 00*log (3.0d 00))         
         facto1 = rrms**2-(7.0d 00/5.0d 00)*(pii**2)*(a**2)
         c = sqrt (5.0d 00/3.0d 00) * sqrt (facto1) 

         ni=1
         nc1=0
         xr2=r/2.0
         do mu=1,mxnmu
	    rr=sqrt(vxisq(mu)+vetasq(ni)-1.d0)
            xrr=xr2*rr
            if (abs(xrr-xr2).lt.c) then 
               nc1=nc1+1
               if (idbg(450).ne.0) write(*,*) 'mu,xrr',mu,xrr
            endif
         enddo

         mu=1
         nc2=0
         do ni=nni,nni/2,-1
	    rr=sqrt(vxisq(mu)+vetasq(ni)-1.d0)
            xrr=xr2*rr
            if (abs(xrr-xr2).lt.c) then
               nc2=nc2+1
               if (idbg(450).ne.0) write(*,*) 'ni,xrr',ni,xrr
            endif
         enddo

         if (iprint.ne.0) then
            write(*,*) '     Center B: '
            write(*,*) '     c = ',c,' bohr', '  a = ',a,' bohr' 
            write(*,*) '     nu=1:',nc1, ' points inside radius c'
            write(*,*) '     mu= 1:',nc2, ' points inside radius c'
         endif
      endif

      return
      end