anchor.diff:  
     if(a.eq.0.d0.and.b.ne.0.d0)then
          
  -    if(n.ne.0)then
  -      int_prod_bessel=0.d0
  -      return
  -    endif
  -
  -    q=0
  -    intold=-1.d0
  -    int=0.d0
  -    done=.false.
  -    kcp=0
  -
  -    do while (.not.done)
  -        
  -      kcp=kcp+1
  -      int=int+coef_nk(m,q)*b**(m+2*q)*crochet(2*q+m+l,gam)
  -      
  -      if(dabs(int-intold).lt.1d-15)then
  -        done=.true.
  -      else
  -        q=q+1
  -        intold=int
  -      endif
  -    
  -    enddo
  -
  +    int = int_prod_bessel_loc(l,gam,m,b)
       int_prod_bessel=int*freal
  -    if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
       return
     endif
   
     if(a.ne.0.d0.and.b.eq.0.d0)then
  -    if(m.ne.0)then
  -      int_prod_bessel=0.d0
  -      return
  -    endif
  -
  -    q=0
  -    intold=-1.d0
  -    int=0.d0
  -    done=.false.
  -    kcp=0
  -
  -    do while (.not.done)
  -      kcp=kcp+1
  -      int=int+coef_nk(n,q)*a**(n+2*q)*crochet(2*q+n+l,gam)
  -      if(dabs(int-intold).lt.1d-15)then
  -         done=.true.
  -      else
  -         q=q+1
  -         intold=int
  -      endif
  -    enddo
  -    
  +
  +    int = int_prod_bessel_loc(l,gam,n,a)
       int_prod_bessel=int*freal
  -    if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
       return
   
     endif
  
positive.code:     intold=-1.d0
      int=0.d0
      done=.false.
      kcp=0
  
      allocate( temp_array_a(0:100))
      allocate( temp_array_b(0+m:200+m))
  
      do while (.not.done)  
  
        kcp=kcp+1
        sum=0.d0
  
        temp_array_a(q) = coef_nk(n,q)*(a**(n+2*q))
        temp_array_b(q) = coef_nk(m,q)*b**(m+2*q)
  
        do k=0,q
          sum=sum+temp_array_a(k)*temp_array_b(q-2*k)
        enddo
  
       int=int+sum*crochet(2*q+n+m+l,gam)
       
       if(dabs(int-intold).lt.1d-15)then
         done=.true.
        else
          
          q=q+1
          intold=int
        endif
  
      enddo
      
      deallocate(temp_array_a)
      deallocate(temp_array_b)
  
      int_prod_bessel=int*freal
      
      if(kcp.gt.100)then
        print*,'**WARNING** bad convergence in int_prod_bessel'
      endif
  
  
negative.code:   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
  