Implement the Mourigal limit

[Lazersmoke]

This computation: k /= norm(k) + 1e-12 introduces an unphysical cutoff, and in particular it disagrees with SpinW.

However, asking for polarization_matrix([0,0,0]) is physically undefined. The correct thing to do is complicated, and not really experimentally accessible. The closest thing is to average over all k in a neighborhood of (but not actually at) zero, as described in the code comment. The other option is to simply report NaN, which makes fewer assumptions, but also makes life much more annoying for plotting purposes.

[kbarros]

Brilliant. Any opportunity to get rid of hard-coded constants is a big win. Who has the knowledge to review this?

[Lazersmoke]

Great question! Martin, Harry, and I had discussed this and decided that the value at zero is actually irrelevant, since it's not physically accessible. Further opinions on this are welcomed!

[kbarros]

Do people agree that this one looks safe to merge? It is only changes the behavior at the special k=0 point and the current behavior seems ad hoc anyway.

My only suggestion is regarding whitespace. Please use:

@test isapprox(is, is_golden; atol=1e-12)

That is, insert whitespace after comma and semicolon, whereas whitespace can be deleted around the = sign in a named argument. This is consistent with examples from BlueStyle: https://github.com/invenia/BlueStyle

[kbarros]

Here is the explanation, lifted from the PR:

        # N.B.: Where does this 2/3 come from??
        # =====================================
        # When deriving the "dipole correction" (δ - q⊗q) [with q a unit vector],
        # the following identity is used:
        #
        #     ∫ exp(i(Q+q)⋅R) d³R = (2π)³δ(Q+q)       [see 4.29 in Boothroyd]
        #
        # This integral is taken over an infinite volume (d³R over all space).
        # In reality, the sample volume is finite, so instead of a sharp delta function
        # on the right hand side, should really have a slightly blurred kernel whose shape
        # depends reciprocally on the shape of the coherent volume within the sample.
        # To zeroth order, this is just a sphere of some finite size.
        #
        # As a result, the (δ - q⊗q) should be averaged over nearby q in some way.
        # For nonzero q, this changes essentially nothing. But for q comparable to
        # 1/(sample length), the averaging is significant.
        #
        # Here, we assume a 3D spherical volume, so that the averaging is entirely
        # isotropic. Averaging the 3x3 matrix diag(0,1,1) over all rotations
        # (isotropically) gives 2/3 the identity matrix.

[kbarros]

Implemented in #288.