LALSuite (2025)

LALSuite (the LVK Algorithm Library Suite) is a collection of routines, comprising the primary official software used by the LIGO-Virgo-KAGRA Collaboration to detect and characterize gravitational waves. As far as I can tell, the ultimate source for all spin-weighted spherical harmonics used in LALSuite is the function XLALSpinWeightedSphericalHarmonic, which cites the NINJA paper [26] as its source. Unfortunately, it cites version 1, which contained a serious error, using $\tfrac{\cos\iota}{2}$ instead of $\cos \tfrac{\iota}{2}$ and similarly for $\sin$. This error was corrected in version 2, but the citation was not updated. Nonetheless, it appears that the actual code is consistent with the corrected versions of the NINJA paper.

They also (quite separately) define Wigner's $D$ matrices in terms of the $d$ matrices, which are — in turn — defined in terms of Jacobi polynomials. For all of these, they cite Wikipedia (despite the fact that the NINJA paper defined the spin-weighted spherical harmonics in terms of the $d$ matrices). Nonetheless, the definitions in the code are consistent with the definitions in the NINJA paper, which are consistent with the definitions in the SphericalFunctions package.

Implementing formulas

Unfortunately, the lalsuite python package cannot be reliably tested on multiple platforms, so we will take a much hackier, but ultimately more robust, approach. We will simply read the C source code, and convert it to Julia code. There are three functions that we will want to compare with the definitions in this package:

COMPLEX16 XLALSpinWeightedSphericalHarmonic( REAL8 theta, REAL8 phi, int s, int l, int m );
double XLALWignerdMatrix( int l, int mp, int m, double beta );
COMPLEX16 XLALWignerDMatrix( int l, int mp, int m, double alpha, double beta, double gam );

The original source code (as of early 2025) is stored alongside this file, so we will read it in to a String and then apply a series of regular expressions to convert it to Julia code, parse it and evaluate it to turn it into runnable Julia. We encapsulate the formulas in a module so that we can test them against the SphericalFunctions package.

We begin by setting up that module, and introducing a set of basic replacements that would usually be defined in separate C headers.

module LALSuiteusing Printf: @sprintfconst I = imconst LAL_PI = πconst XLAL_EINVAL = "XLAL Error: Invalid arguments"XLAL_ERROR_VAL(val, err) = throw(ArgumentError(err))MIN(a, b) = min(a, b)gsl_sf_choose(a, b) = binomial(a, b)pow(a, b) = a^bcexp(a) = exp(a)cpolar(a, b) = a * cis(b)macro XLALPrError(msg, args...)    esc(:(@error @sprintf($msg, $(args...))))end

Next, we simply read the source file into a string.

lalsource = read(joinpath(@__DIR__, "lalsuite_SphericalHarmonics.c"), String)

Now we define a series of replacements to apply to the C code to convert it to Julia code. Note that some of these will be quite specific to this particular file, and may not be generally applicable.

replacements = (    # Deal with annoying Windows line endings    "\r\n" => "\n",    # Deal with newlines in the middle of an assignment    r"( = .*[^;]\s*)\n" => s"\1",    # Remove a couple old, unused functions    r"(?ms)XLALScalarSphericalHarmonic.*?\n}" => "# Removed",    r"(?ms)XLALSphHarm.*?\n}" => "# Removed",    # Remove type annotations    r"COMPLEX16 ?" => "",    r"REAL8 ?" => "",    r"INT4 ?" => "",    r"int ?" => "",    r"double ?" => "",    # Translate comments    "/*" => "#=",    "*/" => "=#",    # Brackets    r" ?{" => "",    r"}.*(\n *else)" => s"\1",    r"} *else" => "else",    r"^}" => "",    "}" => "end",    # Flow control    r"( *if.*);"=>s"\1 end",  ## one-line `if` statements    "for( s=0; n-s >= 0; s++ )" => "for s=0:n",    "else if" => "elseif",    r"(?m)  break;\n *\n *case(.*?):" => s"elseif m == \1",    r"(?m)  break;\n\s*case(.*?):" => s"elseif m == \1",    r"(?m)  break;\n *\n *default:" => "else",    r"(?m)  break;\n *default:" => "else",    r"(?m)switch.*?\n *\n( *)case(.*?):" => s"\n\1if m == \2",    r"\n *break;" => "",    r"(?m)( *ans = fac;)\n" => s"\1\n  end\n",    # Deal with ugly C declarations    "f1 = (x-1)/2.0, f2 = (x+1)/2.0" => "f1 = (x-1)/2.0; f2 = (x+1)/2.0",    "sum=0, val=0" => "sum=0; val=0",    "a=0, lam=0" => "a=0; lam=0",    r"\n *fac;" => "",    r"\n *ans;" => "",    r"\n *gslStatus;" => "",    r"\n *gsl_sf_result pLm;" => "",    r"\n ?XLAL" => "\nfunction XLAL",    # Differences in Julia syntax    "++" => "+=1",    ".*" => ". *",    "./" => ". /",    ".+" => ". +",    ".-" => ". -",    # Deal with random bad syntax    "if (m)" => "if m != 0",    "case 4:" => "elseif m == 4",    "XLALPrError" => "@XLALPrError",    "__func__" => "\"\"",)

And we apply the replacements to the source code to convert it to Julia code. Note that we apply them successively, even though replace can handle multiple "simultaneous" replacements, because the order of replacements is important.

for (pattern, replacement) in replacements    global lalsource = replace(lalsource, pattern => replacement)end

Finally, we just parse the code to turn it into a runnable Julia expression, evaluate it inside this module, and we are done defining the module.

eval(Meta.parseall(lalsource))end  # module LALSuite

Tests

We can now test the LALSuite functions against the equivalent functions from the SphericalFunctions package. We will need to test approximate floating-point equality, so we set absolute and relative tolerances (respectively) in terms of the machine epsilon:

ϵₐ = 100eps()ϵᵣ = 100eps()

The spin-weighted spherical harmonics are defined explicitly, but only for

s = -2

and only up to

ℓₘₐₓ = 8

so we only test up to that point.

for (θ, ϕ) ∈ θϕrange(rng)    for (ℓ, Yˡ) ∈ SphericalFunctions.sYlm(θ, ϕ, ℓₘₐₓ, s), m ∈ -ℓ:ℓ        @test LALSuite.XLALSpinWeightedSphericalHarmonic(θ, ϕ, s, ℓ, m) ≈            Yˡ[m] atol=ϵₐ rtol=ϵᵣ    endend

Now, the Wigner $d$ matrices are defined generally, but we only need to test up to

ℓₘₐₓ = 4

because the formulas are fairly inefficient and inaccurate, and this will be sufficient to sort out any sign or normalization differences, which are the most likely sources of error.

for β ∈ βrange(rng)    for (ℓ, dˡ) ∈ SphericalFunctions.dCalculator(β, ℓₘₐₓ), m′ ∈ -ℓ:ℓ, m ∈ -ℓ:ℓ        @test LALSuite.XLALWignerdMatrix(ℓ, m′, m, β) ≈            dˡ[m′, m] atol=ϵₐ rtol=ϵᵣ    endend

We can see more-or-less by inspection that the code defines the $D$ matrix in agreement with our convention, the key line being

cexp( -(1.0I)*mp*alpha ) * XLALWignerdMatrix( l, mp, m, beta ) * cexp( -(1.0I)*m*gam );

Note that these are the signs this package uses, while its versions before 3.0 used the complex-conjugate convention.

for (α,β,γ) ∈ αβγrange(rng)    for (ℓ, 𝔇ˡ) ∈ SphericalFunctions.DCalculator(α, β, γ, ℓₘₐₓ), m′ ∈ -ℓ:ℓ, m ∈ -ℓ:ℓ        @test LALSuite.XLALWignerDMatrix(ℓ, m′, m, α, β, γ) ≈            𝔇ˡ[m′, m] atol=ϵₐ rtol=ϵᵣ    endend

These successful tests show that the spin-weighted spherical harmonics and the Wigner $d$ and $D$ matrices defined in LALSuite agree with the corresponding functions defined by the SphericalFunctions package.

This page was generated using Literate.jl.