anchor.diff:          else
             wb = (l_max+m_mid)*(l_max+m_mid+1)/2
           end if
  -C        write (6,*) wa,wb
           result = sqrt(wb)/wa
   *
   *   General code for l_min > 1
   *
         else
  +*
  +*     (a+b-c-1)!!(b+c-a-1)!!(c+a-b-1)!!/(a+b+c+1)!!
  +*
           l_sum = l1+l2+l3
  -        ll_max = l_sum
  +        ll_max = 2*l_max
           ll_min = 2*l_min
           ll_mid = 2*l_mid
           m = 1
  
positive.code: *   Calculate product of 3j symbols which appears in the angular       * 
  *   integral over three spherical tensors. Product is                  *
  *   (l1 l2 l3 / m1 m2 m3) (l1 l2 l3 / 0 0 0).                          *
  *   Note that because the maximum angular momentum in the projectors   *
  *   is 3, the only cases needed are l1+l2=l3 and l1+l2-2=l3 and cyclic *
  *   permutations of these.                                             *
  *                                                                      *
  *   l1,l2,l3 - l values                                                *
  *   m1,m2 - m values; m3 = -(m1+m2)                                    *
  
negative.code:         wa = m
          wb = n
          result = result*sqrt(wa/wb)
  *
  *     Remaining factor
  *
          ll = l_min+l_mid-l_max
          if (ll .gt. 0) then
            i = l_min+m_min
            j = l_min-m_min
            m = l_mid+m_mid
            n = l_mid-m_mid
            ll_min = max(0,ll-i,ll-n)
            ll_max = min(ll,j,m)
            k = 1
            do l = 1,ll_min
              k = -k*j*m/l
              j = j-1
              m = m-1
            end do
            do l = 1,ll-ll_min
              k = k*i*n/l
              i = i-1
              n = n-1
            end do
            l_sum = k
            do l = ll_min+1,ll_max
              i = i+1
              n = n+1
              k = -k*j*m*(ll-l+1)/(i*n*l)
              j = j-1
              m = m-1
              l_sum = l_sum+k
            end do
            wa = l_sum
            result = result*wa
          end if
        end if
  *
        if (mod(abs(phase+m_max+ll/2),2) .eq. 1) result = -result
  *
        return
        end
  