anchor.diff:          x=Lx(i)*Hx-(Rion(1,a)+ix*aLx)
           y=Ly(i)*Hy-(Rion(2,a)+iy*aLy)
           z=Lz(i)*Hz-(Rion(3,a)+iz*aLz)
  -        r=sqrt(x*x+y*y+z*z)
  +        r=sqrt(x**2+y**2+z**2)
           if (r<Rps(ik)) then
             j=j+1
           endif
  
positive.code:         do m=-l,l
            lm=lm+1
            uV(j,lma_tbl(lm,a))=uVr(l)*Ylm(x,y,z,l,m)
            duV(j,lma_tbl(lm,a),1)=duVr(l)*(x/r)*Ylm(x,y,z,l,m)+uVr(l)*dYlm(x,y,z,l,m,1)
            duV(j,lma_tbl(lm,a),2)=duVr(l)*(y/r)*Ylm(x,y,z,l,m)+uVr(l)*dYlm(x,y,z,l,m,2)
            duV(j,lma_tbl(lm,a),3)=duVr(l)*(z/r)*Ylm(x,y,z,l,m)+uVr(l)*dYlm(x,y,z,l,m,3)
          enddo
        enddo
      enddo
  
negative.code:         x=Lx(i)*Hx-(Rion(1,a)+ix*aLx)
          y=Ly(i)*Hy-(Rion(2,a)+iy*aLy)
          z=Lz(i)*Hz-(Rion(3,a)+iz*aLz)
          r=sqrt(x**2+y**2+z**2)
          if (r<Rps(ik)) then
            j=j+1
            if (j<=Nps) then
  