You mean this bit here:

            Wk[:] = np.dot(linalg.pinv(Ck, 0.1), Wk)

right?

Well, to me the only thing it does is that it includes in the final elements of the normal component of `W` (`W[2::3]`) the normal component of `Wk` scaled by variance components and both of the normal components scaled by covariance components. So, as it has to be, `W[2::3]` is mostly the scalar beamformer solution that we would get with normal orientation being applied much earlier on plus some additional components because we're now dealing with 3D solutions. Does this make any sense? I'm not satisfied with this explanation entirely, it just shows I haven't yet developed the 'feel' for the beamformer equations that I'd like to...

...but, there's another reason why I'm sure picking `W[2::3]` gives the normal component - the solution (the plot in the first example) is almost identical to what the solution was when the normal orientation picking was still implemented as a scalar beamformer (where this problem disappeared, there was no rotation, just division by a scalar), whereas `W[0::3]` and `W[1::3]` give very different results.