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:   integer n,k,m,q,l,kcp
    double precision gam,dble_fact,fact,pi,a,b
    double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
    double precision dump
  
    double precision, allocatable :: temp_array_a(:)
    double precision, allocatable :: temp_array_b(:)
  
    logical done
  
  
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
  
  