anchor.diff:       prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3))
        
        do k3p=0,n_b(3)
  +       if (array_coefs_B(k3p,3) == 0.d0) cycle
          do k2p=0,n_b(2)
  +         if (array_coefs_B(k2p,2) == 0.d0) cycle
            do k1p=0,n_b(1)
  +           if (array_coefs_B(k1p,1) == 0.d0) cycle
              do lambdap=0,l+ntotB
                do mup=-lambdap,lambdap
                  
                  prodp=prod*array_coefs_B(k1p,1)*array_coefs_B(k2p,2)*array_coefs_B(k3p,3)*ylm(lambdap,mup,theta_BC0,phi_BC0)*array_I_B(mup,lambdap,k1p,k2p,k3p)
                  
  +               if (prodp == 0.d0) cycle
  +
                  do k=1,kmax
                    
                    ktot=ntotA+k1p+k2p+k3p+n_kl(k,l)
  
positive.code: mu=-mu
  m=-m
  factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
  factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
  sum=0.d0
  sgn = 1.d0
  do k=0,(mu-1)/2
   do i=0,lambda-mu
    sgnp = 1.d0
    do kp=0,(m-1)/2
     do ip=0,l-m
      cylm=sgn*factor1*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
      cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
      sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-(2*kp+1)+k1,(2*k+1)+(2*kp+1)+k2,i+ip+k3)
     enddo
     sgnp = -sgnp
  
negative.code:    phi_DC0=datan2(d(2)/d2,d(1)/d2)
     
     do k=1,klocmax
       if (v_k(k) == 0.d0) cycle
       do k1=0,n_a(1)
         do k2=0,n_a(2)
           do k3=0,n_a(3)
             do k1p=0,n_b(1)
               do k2p=0,n_b(2)
                 do k3p=0,n_b(3)
                   if (array_coefs(k1,k2,k3,k1p,k2p,k3p) == 0.d0) cycle
                   do l=0,ntot
                     do m=-l,l
                       coef=ylm(l,m,theta_DC0,phi_DC0)
                       if (coef == 0.d0) cycle
                       ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k)
                       if (array_R_loc(ktot,k,l) == 0.d0) cycle
                       prod=coef*array_coefs(k1,k2,k3,k1p,k2p,k3p)     &
                           *bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p)
                       accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l)
                     enddo
                   enddo
  