anchor.diff:        do l=0,lmax
   
           array_R(ktot,k,l,0,lambdap)= freal & 
  -                                     *int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal)
  -      enddo
  +                                     * int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg)
  +       enddo
       enddo
     enddo
    enddo
  
positive.code:       stop 'pb in int_prod_bessel!!'
        end
  
  !!    Calculation of
  !!
  !!    I= \int dx x**l *exp(-gam*x**2) M_n(ax)
  
negative.code:       do l=0,lmax
  
        array_R(ktot,k,l,lambda,0)= freal & 
                                    * int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg)
        enddo
      enddo
    enddo
  