      subroutine fill_bas(i)
      use grid_module
      use orbitals_module
      use basins_module
      implicit none
      integer npx,npy,npz,indy,ixa,iya,iza,i,j,k,irw,icircle
      integer ifilled,level
      double precision rn,xp,yp,zp,dist,stepx,stepy,stepz
      double precision distmin
      double precision xmini,xmaxi,ymini,ymaxi,zmini,zmaxi
      double precision spreadx,spready,xn,yn,zn,b2,theta,an
      double precision xc,yc,zc,xm,ym,zm,dtheta,rmin,rmax,delt,pi
      double precision xstart,ystart,zstart,xend,yend,zend
      double precision rho,drho(3),d2rho(6),elfd,gelf(3),d2elf(6)
      pi=dacos(-1.0d0)
      dtheta=0.01d0
      npx=np(1)-1
      npy=np(2)-1
      npz=np(3)-1
      ifilled=0
      if(npx.eq.0) then
         stepx=0.0d0
      else
         stepx=(xmax-xmin)/npx
      endif
      if(npy.eq.0) then
         stepy=0.0d0
      else
         stepy=(ymax-ymin)/npy
      endif
      if(npz.eq.0) then
         stepz=0.0d0
      else
         stepz=(zmax-zmin)/npz
      endif
      xmini=xmin-0.4d0*stepx
      ymini=ymin-0.4d0*stepy
      zmini=zmin-0.4d0*stepz
      xmaxi=xmax+0.4d0*stepx
      ymaxi=ymax+0.4d0*stepy
      zmaxi=zmax+0.4d0*stepz
      icircle=iabs(basin_type(i,1))
*
*     fill around core
*
      if(icircle.eq.0) then
         rmax=0.95d0*core_radius(i)
         xstart=x(i)-1.2d0*rmax
         ystart=y(i)-1.2d0*rmax
         zstart=z(i)-1.2d0*rmax
         xend=x(i)+1.2d0*rmax
         yend=y(i)+1.2d0*rmax
         zend=z(i)+1.2d0*rmax
         if(xmin.ge.xstart) xstart=xmin
         if(ymin.ge.ystart) ystart=ymin
         if(zmin.ge.zstart) zstart=zmin
         if(xmax.le.xend) xend=xmax
         if(ymax.le.yend) yend=ymax
         if(zmax.le.zend) zend=zmax
         do xp=xstart,xend,stepx
            do yp=ystart,yend,stepy
               do zp=zstart,zend,stepz
                  rmin=dsqrt((xp-x(i))**2+(yp-y(i))**2+(zp-z(i))**2)
                  if(rmin.le.rmax) then
                        ixa=nint((xp-xmin)/stepx)
                        iya=nint((yp-ymin)/stepy)
                        iza=nint((zp-zmin)/stepz)
                        indy=1+iza+np(3)*(iya+np(2)*ixa)
                        if(indy.le.0) then
                            write(*,'(3f12.6)') xp,yp,zp
                            stop 'core '
                        endif
                        if(attract_code(indy).eq.0) then
                           attract_code(indy)=i
                           ifilled=ifilled+1
                        endif
                  endif
               enddo
            enddo
         enddo
      else 
*
*        basin (core or valence) with point, circular or spherical attractor
*
         j=basin_type(i,3)
         if(j.le.0) then
            j=1
         endif
         xn=x(i)-xat(j)
         yn=y(i)-yat(j)
         zn=z(i)-zat(j)
         b2=xn*xn+yn*yn
         rn=core_radius(i)
         if(icircle.le.2) rn=0.0d0
         if(b2.le.1.0d-2) then
            xc=0.0d0
            yc=rn
            zc=0.0d0
            xm=rn
            ym=0.0d0
            zm=0.0d0
         else
            an=rn/dsqrt(b2)
            xm=an*yn
            ym=-an*xm
            zm=0.0d0
            dist=rn/dsqrt(b2*zn*zn+b2*b2)
            xc=-dist*xn*zn
            yc=-dist*yn*zn
            zc=dist*b2
         endif
         if(natoms.eq.1) then
            icircle=4
         endif
*
*       point and spherical attractors
*
         if(icircle.ne.3) then
            irw=0
            k=0
            distmin=1.0d6
            do k=1,natoms
               dist=(x(i)-xat(k))**2+(y(i)-yat(k))**2+(z(i)-zat(k))**2
               if((dist.lt.distmin).and.(dist.gt.core_radius(i))) then
                  irw=k
                  distmin=dist
               endif
            enddo
            rmax=rn
            yp=y(i)
            zp=z(i)
            xp=x(i)
            elfd=val_elf(i)
            delt=0.95d0*val_elf(i)
            do while(elfd.ge.delt)
               rmax=rmax+0.04d0
               xp=x(i)+rmax*(xat(irw)-x(i))/distmin
               yp=y(i)+rmax*(yat(irw)-y(i))/distmin
               zp=z(i)+rmax*(zat(irw)-z(i))/distmin
               level=0
               call elf_eval(rho,drho,d2rho,elfd,gelf,d2elf,xp,yp,zp,
     .                       level)
               if(rmax.gt.6.0d0) then
                    elfd=0.0d0
               endif
            enddo
            if (icircle.le.2) then
               rmin=0.0d0
               xp=x(i)
               if((xp.ge.xmini).and.(xp.le.xmaxi).and.
     .            (yp.ge.ymini).and.(yp.le.ymaxi).and.
     .            (zp.ge.zmini).and.(zp.le.zmaxi)) then
                  ixa=nint((xp-xmin)/stepx)
                  iya=nint((yp-ymin)/stepy)
                  iza=nint((zp-zmin)/stepz)
                  indy=1+iza+np(3)*(iya+np(2)*ixa)
                  if(indy.le.0) then 
                         write(*,'(3f12.6)') xp,yp,zp
                         stop 'icircle.le.2 '
                     endif
                  attract_code(indy)=i
                  ifilled=ifilled+1
               endif
            else
               rmin=1.1d0*core_radius(irw)
            endif
            ifilled=0
            xstart=x(i)-1.2d0*rmax
            ystart=y(i)-1.2d0*rmax
            zstart=z(i)-1.2d0*rmax
            xend=x(i)+1.2d0*rmax
            yend=y(i)+1.2d0*rmax
            zend=z(i)+1.2d0*rmax
            if(xmin.ge.xstart) xstart=xmin
            if(ymin.ge.ystart) ystart=ymin
            if(zmin.ge.zstart) zstart=zmin
            if(xmax.le.xend) xend=xmax
            if(ymax.le.yend) yend=ymax
            if(zmax.le.zend) zend=zmax
            do xp=xstart,xend,stepx
               do yp=ystart,yend,stepy
                  do zp=zstart,zend,stepz
                     xc=dsqrt((xp-x(i))**2+(yp-y(i))**2+(zp-z(i))**2)
                     if((xc.gt.rmin).and.(xc.lt.rmax)) then
                        ixa=nint((xp-xmin)/stepx)
                        iya=nint((yp-ymin)/stepy)
                        iza=nint((zp-zmin)/stepz)
                        indy=1+iza+np(3)*(iya+np(2)*ixa)
                        if(indy.le.0) then 
                            write(*,'(3f12.6)') xp,yp,zp
                            stop 'sphere '
                        endif
                        if((attract_code(indy).eq.0).and.
     .                       (dble(elf(indy)).ge.delt)) then
                           attract_code(indy)=i
                           ifilled=ifilled+1
                        endif
                     endif
                  enddo
               enddo
            enddo
         else
*
*      circular attractor
*
            b2=dsqrt(b2+zn*zn)
            xn=xn/b2
            yn=yn/b2
            zn=zn/b2
            do spreadx=-0.2d0,0.2d0,0.05d0
               do spready=0.7d0,1.3d0,0.05d0
                  do theta=0.0d0,2.0d0*pi,dtheta
                     xp=x(i)+(xm*dcos(theta)+xc*dsin(theta))*spready+
     .                       spreadx*xn
                     yp=y(i)+(ym*dcos(theta)+yc*dsin(theta))*spready+
     .                       spreadx*yn
                     zp=z(i)+(zm*dcos(theta)+zc*dsin(theta))*spready+
     .                       spreadx*zn
                     dist=dsqrt((xp-x(i))**2+(yp-y(i))**2+(zp-z(i))**2)
                     if((xp.ge.xmini).and.(xp.le.xmaxi).and.
     .                  (yp.ge.ymini).and.(yp.le.ymaxi).and.
     .                  (zp.ge.zmini).and.(zp.le.zmaxi)) then
                        ixa=nint((xp-xmin)/stepx)
                        iya=nint((yp-ymin)/stepy)
                        iza=nint((zp-zmin)/stepz)
                        indy=1+iza+np(3)*(iya+np(2)*ixa)
                        if(indy.le.0) then 
                         write(*,'(3f12.6)') xp,yp,zp
                         stop 'cercle '
                     endif
                        if(attract_code(indy).eq.0) then
                           attract_code(indy)=i
                           ifilled=ifilled+1
                        endif
                     endif
                  enddo
               enddo
            enddo
         endif
      endif
      return
      end
