anchor.diff:         double precision x,coef,accu,fact,dble_fact
          accu=0.d0
          do k=0,10
  -        coef=1.d0/fact(k)/dble_fact(2*(n+k)+1)
  -        accu=accu+(x**2/2.d0)**k*coef
  +        coef=1.d0/(fact(k)*dble_fact(2*(n+k)+1))
  +        accu=accu+(0.5d0*x*x)**k*coef
          enddo
          bessel_mod_exp=x**n*accu
          end
  
positive.code:                             array_coefs_B(k1p,1)*array_coefs_B(k2p,2)*array_coefs_B(k3p,3)* &
                              array_I_B(mup,lambdap,k1p,k2p,k3p)
                           
                           do k=1,kmax
                             ktot=k1+k2+k3+k1p+k2p+k3p+n_kl(k,l)
                             accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,lambda,lambdap)
  
negative.code:     enddo
  
      do k3=0,n_a(3)
        if (array_coefs_A(k3,3) == 0.d0) cycle
        do k2=0,n_a(2)
          if (array_coefs_A(k2,2) == 0.d0) cycle
          do k1=0,n_a(1)
            if (array_coefs_A(k1,1) == 0.d0) cycle
            do lambda=0,l+ntotA
              do mu=-lambda,lambda
                
                prod=array_coefs_A(k1,1)*array_coefs_A(k2,2)*array_coefs_A(k3,3)*ylm(lambda,mu,theta_AC0,phi_AC0)*array_I_A(mu,lambda,k1,k2,k3)
                if (prod == 0.d0) cycle
                prodp=prod*bigI(0,0,l,m,n_b(1),n_b(2),n_b(3))
  
                if (prodp == 0.d0) cycle
                
                do k=1,kmax
                  ktot=k1+k2+k3+ntotB+n_kl(k,l)
  