anchor.diff:      intold=-1.d0
       int=0.d0
       done=.false.
  -    kcp=0
  -
  -    allocate( temp_array_a(0:100))
  -    allocate( temp_array_b(0+m:200+m))
  -
  +
  +    n_1 = 2*(n)+1
  +    m_1 = 2*m+1
  +    nlm = n+m+l
  +    pi=dacos(-1.d0)
  +    a_over_b_square = (a/b)**2
  +
  +    ! Calcul first term of the sequence
  +
  +    term_a =dble_fact(nlm-1) / (dble_fact(n_1)*dble_fact(m_1))
  +    expo=0.5d0*dfloat(nlm+1)
  +    term_rap = term_a / (2.d0*gam)**expo
  +
  +    s_0_0=term_rap*a**(n)*b**(m)
  +    if(mod(nlm,2).eq.0)s_0_0=s_0_0*dsqrt(pi/2.d0)
  +
  +    ! Initialise the first recurence terme for the q loop
  +    s_q_0 = s_0_0
  +
  +    ! Loop over q for the convergence of the sequence
       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)
  -
  +      ! Init
  +      sum=0
  +      s_q_k=s_q_0
  +
  +      ! Iteration of k
         do k=0,q
  -        sum=sum+temp_array_a(k)*temp_array_b(q-2*k)
  +        sum=sum+s_q_k
  +        s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)/dble(2*(k+n)+3) ) * ( dble(q-k)/dble(k+1)) * s_q_k
         enddo
  -
  -     int=int+sum*crochet(2*q+n+m+l,gam)
  -     
  -     if(dabs(int-intold).lt.1d-15)then
  -       done=.true.
  +    
  +      int=int+sum
  +
  +      if(dabs(int-intold).lt.1d-15)then
  +        done=.true.
         else
  +
  +        !Compute the s_q+1_0
  +        s_q_0=s_q_0*(2.d0*q+nlm+1)*b**2/((2.d0*(m+q)+3)*4.d0*(q+1)*gam)
           
  +        if(mod(n+m+l,2).eq.1)s_q_0=s_q_0*dsqrt(pi/2.d0)
  +        ! Increment q
           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
  +
  +    if(q.gt.100)then
         print*,'**WARNING** bad convergence in int_prod_bessel'
       endif
   
  
positive.code: 
    double precision gam,dble_fact,fact
  
    double precision, save :: store_result(0:2,0:10)
  
    store_result(0, 0) =  1.00000000000000d0
    store_result(0, 1) =  0.166666666666667d0
    store_result(0, 2) =  8.333333333333333d-3
    store_result(0, 3) =  1.984126984126984d-4
    store_result(0, 4) =  2.755731922398589d-6
    store_result(0, 5) =  2.505210838544172d-8
    store_result(0, 6) =  1.605904383682161d-10
    store_result(0, 7) =  7.647163731819816d-13
    store_result(0, 8) =  2.811457254345521d-15
    store_result(0, 9) =  8.220635246624329d-018
    store_result(0,10) =  1.957294106339126d-020
  
    store_result(1, 0) =  0.333333333333333d0     
    store_result(1, 1) =  3.333333333333333d-002   
    store_result(1, 2) =  1.190476190476191d-003
    store_result(1, 3) =  2.204585537918871d-005   
    store_result(1, 4) =  2.505210838544172d-007   
    store_result(1, 5) =  1.927085260418594d-009   
    store_result(1, 6) =  1.070602922454774d-011   
    store_result(1, 7) =  4.498331606952833d-014   
    store_result(1, 8) =  1.479714344392379d-016
    store_result(1, 9) =  3.914588212678252d-019
    store_result(1,10) =  8.509974375387505d-022
  
    store_result(2, 0) =  6.666666666666667d-002
    store_result(2, 1) =  4.761904761904762d-003
    store_result(2, 2) =  1.322751322751323d-004
    store_result(2, 3) =  2.004168670835337d-006
    store_result(2, 4) =  1.927085260418594d-008
    store_result(2, 5) =  1.284723506945729d-010
    store_result(2, 6) =  6.297664249733966d-013
    store_result(2, 7) =  2.367542951027807d-015
    store_result(2, 8) =  7.046258782820854d-018
    store_result(2, 9) =  1.701994875077501d-020
    store_result(2,10) =  3.403989750155002d-023
  
    if (n.LE.2) then
      if (k.LE.10) then
        coef_nk = store_result(n,k)
        return 
      endif
    endif
  
    if ((n+k).GE.60) then
      coef_nk = 0.d0
      return
    endif
  
    gam=dble_fact(2*(n+k)+1)
    coef_nk=1.d0/(dble(ISHFT(1,k))*fact(k)*gam)
    
    return
  
  
negative.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
  
    integer :: n_1,m_1,nlm
    double precision :: term_A, term_B, term_rap,  expo
    double precision :: s_q_0, s_q_k, s_0_0, a_over_b_square
    double precision :: int_prod_bessel_loc
  
    logical done
  
  