anchor.diff:    do lambdap=0,lmax+ntotB
       do k=1,kmax
         do l=0,lmax
  -
           array_R(ktot,k,l,0,lambdap)=  int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg)
          enddo
       enddo
  
positive.code: !!
  !!    M_n(x) modified spherical bessel function 
  !!
        double precision function int_prod_bessel_loc(l,gam,n,a)
        implicit none
        integer n,k,l,kcp
        double precision gam,a
        double precision int,intold,coef_nk,crochet
        logical done
        k=0
        intold=-1.d0
        int=0.d0
        done=.false.
        kcp=0
        do while (.not.done)
         kcp=kcp+1
         int=int+coef_nk(n,k)*a**(n+2*k)*crochet(2*k+n+l,gam)
         if(dabs(int-intold).lt.1d-15)then
          done=.true.
         else
          k=k+1
          intold=int
         endif
        enddo
        int_prod_bessel_loc=int
        if(kcp.gt.100)print*,'**WARNING** bad convergence in int_prod_bessel'
        end
  
        double precision function int_prod_bessel_num(l,gam,n,m,a,b)
        implicit none
  
negative.code: 
    double precision gam,dble_fact,fact
  
    gam=dble_fact(2*(n+k)+1)
  
  !  coef_nk=1.d0/(dble(ISHFT(1,k))*fact(k)*gam)
    
    coef_nk=1.d0/(2.d0**k*fact(k)*gam)
    
    return
  
  