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: sum=0.d0
  mu=-mu
  factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
  factor2=dsqrt((2*l+1)/(4.d0*pi))
  sgn = 1.d0
  do k=0,(mu-1)/2
   do i=0,lambda-mu
     do ip=0,l
      cylm=sgn*factor1*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
      cylmp=factor2*coef_pm(l,ip)
      sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+k1,2*k+1+k2,i+ip+k3)
     enddo
   enddo
  
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=ylm(lambda,mu,theta_AC0,phi_AC0)*array_coefs_A(k1,1)*array_coefs_A(k2,2)*array_coefs_A(k3,3)*array_I_A(mu,lambda,k1,k2,k3)
                 if (prod == 0.d0) cycle
                 
                 do k3p=0,n_b(3)
                   do k2p=0,n_b(2)
  