anchor.diff:        stop 'pb in int_prod_bessel!!'
         end
   
  +
  +      double precision function int_prod_bessel_num_soph_p(l,gam,n,m,a,b,arg)
  +      implicit none
  +      integer n,m,l
  +      double precision gam,a,b,arg,arg_new
  +      double precision bessel_mod,factor
  +logical not_done
  +double precision bigA,xold,x,dx,accu,intnew,intold,intold2,u,v,freal
  +integer iter
  +double precision n0,nI,i,eps
  +      u=(a+b)/(2.d0*dsqrt(gam))
  +      arg_new=u**2-arg
  +      freal=dexp(arg_new)
  +      v=u/dsqrt(gam)
  +
  +bigA=v+dsqrt(-dlog(1.d-15)/gam)
  +n0=5
  +accu=0.d0
  +dx=bigA/(n0-1.d0)
  +iter=0
  +do i=1.d0,n0
  + x=(i-1.d0)*dx
  + accu=accu+x**l*dexp(-gam*(x-v)**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x)
  +enddo
  +accu=accu*freal
  +intold=accu*dx
  +
  +eps=1.d-08
  +nI=n0-1.d0
  +dx=dx/2.d0
  +not_done=.true.
  +
  +do while(not_done)
  + iter=iter+1
  + accu=0.d0
  + do i=1.d0,nI
  +  x=dx+(i-1.d0)*2.d0*dx
  +  accu=accu+dx*x**l*dexp(-gam*(x-v)**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x)
  + enddo
  + accu=accu*freal
  + intnew=intold/2.d0+accu
  + if(iter.gt.1.and.dabs(intnew-intold).lt.eps.and.dabs(intnew-intold2).lt.eps)then
  +  not_done=.false.
  + else
  +  intold2=intold
  +  intold=intnew
  +  dx=dx/2.d0
  +  nI=2.d0*nI
  + endif
  +enddo
  +int_prod_bessel_num_soph_p=intold
  +end
  +
   !!    Calculation of
   !!
   !!    I= \int dx x**l *exp(-gam*x**2) M_n(ax)
  
positive.code: ! __   __  _ __  ___  ___ _   _  __| | ___  
  ! \ \ / / | '_ \/ __|/ _ \ | | |/ _` |/ _ \ 
  !  \ V /  | |_) \__ \  __/ |_| | (_| | (_) |
  !   \_/   | .__/|___/\___|\__,_|\__,_|\___/ 
  !         | |                               
  !         |_|                               
  
  
negative.code:       do l=0,lmax
  
          array_R(ktot,k,l,0,lambdap)= freal & 
                                       * 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
  