202 lines
7.1 KiB
Go
202 lines
7.1 KiB
Go
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
package core
|
||
|
||
import (
|
||
"math"
|
||
"testing"
|
||
)
|
||
|
||
// TestOrthogonalPolynomials checks H, L, T and U against their closed
|
||
// forms for the first degrees.
|
||
func TestOrthogonalPolynomials(t *testing.T) {
|
||
x := mustFloats(t, []float64{-0.7, 0, 0.4, 1.1})
|
||
cases := []struct {
|
||
name string
|
||
n int
|
||
want func(v float64) float64
|
||
}{
|
||
{"H0", 0, func(v float64) float64 { return 1 }},
|
||
{"H1", 1, func(v float64) float64 { return 2 * v }},
|
||
{"H2", 2, func(v float64) float64 { return 4*v*v - 2 }},
|
||
{"H3", 3, func(v float64) float64 { return 8*v*v*v - 12*v }},
|
||
{"T0", 0, func(v float64) float64 { return 1 }},
|
||
{"T2", 2, func(v float64) float64 { return 2*v*v - 1 }},
|
||
{"T3", 3, func(v float64) float64 { return 4*v*v*v - 3*v }},
|
||
{"U0", 0, func(v float64) float64 { return 1 }},
|
||
{"U1", 1, func(v float64) float64 { return 2 * v }},
|
||
{"U2", 2, func(v float64) float64 { return 4*v*v - 1 }},
|
||
}
|
||
for _, tc := range cases {
|
||
var got *Array
|
||
var err error
|
||
switch {
|
||
case tc.name[0] == 'H' && len(tc.name) == 2:
|
||
got, err = Hermite(tc.n, x)
|
||
case tc.name[0] == 'T':
|
||
got, err = ChebyshevT(tc.n, x)
|
||
default:
|
||
got, err = ChebyshevU(tc.n, x)
|
||
}
|
||
if err != nil {
|
||
t.Fatalf("%s: %v", tc.name, err)
|
||
}
|
||
for i := range x.Len() {
|
||
want := tc.want(x.FloatAt(i))
|
||
if math.Abs(got.FloatAt(i)-want) > 1e-12*(1+math.Abs(want)) {
|
||
t.Fatalf("%s(%v) = %v, want %v", tc.name, x.FloatAt(i), got.FloatAt(i), want)
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestLaguerrePolynomials checks L against closed forms and the
|
||
// orthogonality weight sanity at a sample point.
|
||
func TestLaguerrePolynomials(t *testing.T) {
|
||
x := mustFloats(t, []float64{-0.3, 0, 0.6, 2.1})
|
||
closed := []struct {
|
||
n int
|
||
alpha float64
|
||
want func(v float64) float64
|
||
}{
|
||
{0, 0, func(v float64) float64 { return 1 }},
|
||
{1, 0, func(v float64) float64 { return 1 - v }},
|
||
{2, 0, func(v float64) float64 { return 1 - 2*v + v*v/2 }},
|
||
{2, 1, func(v float64) float64 { return 3 - 3*v + v*v/2 }},
|
||
}
|
||
for _, tc := range closed {
|
||
l, err := Laguerre(tc.n, tc.alpha, x)
|
||
if err != nil {
|
||
t.Fatalf("Laguerre(%d, %v): %v", tc.n, tc.alpha, err)
|
||
}
|
||
for i := range x.Len() {
|
||
want := tc.want(x.FloatAt(i))
|
||
if math.Abs(l.FloatAt(i)-want) > 1e-12*(1+math.Abs(want)) {
|
||
t.Fatalf("L_%d^%v(%v) = %v, want %v", tc.n, tc.alpha, x.FloatAt(i), l.FloatAt(i), want)
|
||
}
|
||
}
|
||
}
|
||
if _, err := Laguerre(-1, 0, x); err == nil {
|
||
t.Fatal("expected an error for a negative degree")
|
||
}
|
||
}
|
||
|
||
// TestSphericalHarmonicReal checks the real harmonics against the
|
||
// complex ones under the standard convention: both branches carry the
|
||
// (−1)^m phase, Y_real = √2·(−1)^m·Re(Y_l^m) for m > 0 and
|
||
// √2·(−1)^m·Im(Y_l^|m|) for m < 0.
|
||
func TestSphericalHarmonicReal(t *testing.T) {
|
||
theta := mustFloats(t, []float64{0.4, 1.2, 2.5})
|
||
phi := mustFloats(t, []float64{0.3, -0.8, 2.1})
|
||
pos, err := SphericalHarmonicReal(2, 1, theta, phi)
|
||
if err != nil {
|
||
t.Fatalf("SphericalHarmonicReal(+1): %v", err)
|
||
}
|
||
pos2, err := SphericalHarmonicReal(3, 2, theta, phi)
|
||
if err != nil {
|
||
t.Fatalf("SphericalHarmonicReal(+2): %v", err)
|
||
}
|
||
neg, err := SphericalHarmonicReal(2, -1, theta, phi)
|
||
if err != nil {
|
||
t.Fatalf("SphericalHarmonicReal(−1): %v", err)
|
||
}
|
||
cx, err := SphericalHarmonic(2, 1, theta, phi)
|
||
if err != nil {
|
||
t.Fatalf("SphericalHarmonic: %v", err)
|
||
}
|
||
cx2, err := SphericalHarmonic(3, 2, theta, phi)
|
||
if err != nil {
|
||
t.Fatalf("SphericalHarmonic(3, 2): %v", err)
|
||
}
|
||
root2 := math.Sqrt2
|
||
for i := range theta.Len() {
|
||
z := cx.ComplexAt(i)
|
||
// m = 1 is odd, so the standard real form negates the real
|
||
// part; m = 2, even, keeps it.
|
||
if math.Abs(pos.FloatAt(i)+root2*real(z)) > 1e-12 {
|
||
t.Fatalf("real m=1 [%d] = %v, want %v", i, pos.FloatAt(i), -root2*real(z))
|
||
}
|
||
z2 := cx2.ComplexAt(i)
|
||
if math.Abs(pos2.FloatAt(i)-root2*real(z2)) > 1e-12 {
|
||
t.Fatalf("real m=2 [%d] = %v, want %v", i, pos2.FloatAt(i), root2*real(z2))
|
||
}
|
||
if math.Abs(neg.FloatAt(i)+root2*imag(z)) > 1e-12 {
|
||
t.Fatalf("real m=−1 [%d] = %v, want %v", i, neg.FloatAt(i), -root2*imag(z))
|
||
}
|
||
}
|
||
// Pinned standard values against the external Cartesian forms:
|
||
// Y_real(1, 1) = √(3/(4π))·x/r, so it is +√(3/(4π)) at (θ, φ) =
|
||
// (π/2, 0) and 0 at (π/2, π/2); Y_real(1, −1) = √(3/(4π))·y/r is
|
||
// +√(3/(4π)) at (π/2, π/2) and 0 at (π/2, 0); Y_real(2, 1) =
|
||
// √(15/(4π))·xz/r² at (π/4, 0) is √(15/(4π))·cos(π/4)·sin(π/4).
|
||
eq := mustFloats(t, []float64{math.Pi / 2})
|
||
halfPi := mustFloats(t, []float64{math.Pi / 2})
|
||
zero := mustFloats(t, []float64{0})
|
||
y11, err := SphericalHarmonicReal(1, -1, eq, halfPi)
|
||
if err != nil {
|
||
t.Fatalf("SphericalHarmonicReal(1, −1): %v", err)
|
||
}
|
||
if want := math.Sqrt(3 / (4 * math.Pi)); math.Abs(y11.FloatAt(0)-want) > 1e-12 {
|
||
t.Fatalf("Y_real(1, −1)(π/2, π/2) = %.16g, want %.16g", y11.FloatAt(0), want)
|
||
}
|
||
y11z, err := SphericalHarmonicReal(1, -1, eq, zero)
|
||
if err != nil {
|
||
t.Fatalf("SphericalHarmonicReal(1, −1) at φ=0: %v", err)
|
||
}
|
||
if v := y11z.FloatAt(0); math.Abs(v) > 1e-12 {
|
||
t.Fatalf("Y_real(1, −1)(π/2, 0) = %v, want 0", v)
|
||
}
|
||
y11p, err := SphericalHarmonicReal(1, 1, eq, zero)
|
||
if err != nil {
|
||
t.Fatalf("SphericalHarmonicReal(1, 1): %v", err)
|
||
}
|
||
if want := math.Sqrt(3 / (4 * math.Pi)); math.Abs(y11p.FloatAt(0)-want) > 1e-12 {
|
||
t.Fatalf("Y_real(1, 1)(π/2, 0) = %.16g, want +%.16g", y11p.FloatAt(0), want)
|
||
}
|
||
y11q, err := SphericalHarmonicReal(1, 1, eq, halfPi)
|
||
if err != nil {
|
||
t.Fatalf("SphericalHarmonicReal(1, 1) at φ=π/2: %v", err)
|
||
}
|
||
if v := y11q.FloatAt(0); math.Abs(v) > 1e-12 {
|
||
t.Fatalf("Y_real(1, 1)(π/2, π/2) = %v, want 0", v)
|
||
}
|
||
qr := mustFloats(t, []float64{math.Pi / 4})
|
||
y21, err := SphericalHarmonicReal(2, 1, qr, zero)
|
||
if err != nil {
|
||
t.Fatalf("SphericalHarmonicReal(2, 1): %v", err)
|
||
}
|
||
if want := math.Sqrt(15/(4*math.Pi)) * math.Cos(math.Pi/4) * math.Sin(math.Pi/4); math.Abs(y21.FloatAt(0)-want) > 1e-12 {
|
||
t.Fatalf("Y_real(2, 1)(π/4, 0) = %.16g, want %.16g", y21.FloatAt(0), want)
|
||
}
|
||
}
|
||
|
||
// TestAiry checks the Airy pair against tabulated values at x = 0,
|
||
// ±1 and the NaN contract outside the series window.
|
||
func TestAiry(t *testing.T) {
|
||
x := mustFloats(t, []float64{0, 1, -1, 3.5})
|
||
ai, bi, err := Airy(x)
|
||
if err != nil {
|
||
t.Fatalf("Airy: %v", err)
|
||
}
|
||
wantAI := []float64{0.3550280538878172, 0.13529241631288142, 0.5355608832928901, 0.002584098786988065}
|
||
wantBI := []float64{0.6149266274460007, 1.2074235949528713, 0.1039973894969447, 33.05550674355921}
|
||
for i := range 4 {
|
||
if math.Abs(ai.FloatAt(i)-wantAI[i]) > 1e-12*(1+math.Abs(wantAI[i])) {
|
||
t.Fatalf("Ai(%v) = %.16g, want %.16g", x.FloatAt(i), ai.FloatAt(i), wantAI[i])
|
||
}
|
||
if math.Abs(bi.FloatAt(i)-wantBI[i]) > 5e-10*(1+math.Abs(wantBI[i])) {
|
||
t.Fatalf("Bi(%v) = %.16g, want %.16g", x.FloatAt(i), bi.FloatAt(i), wantBI[i])
|
||
}
|
||
}
|
||
// Outside the series window the contract is NaN, not silence.
|
||
out := mustFloats(t, []float64{9})
|
||
ai9, bi9, err := Airy(out)
|
||
if err != nil {
|
||
t.Fatalf("Airy(9): %v", err)
|
||
}
|
||
if !math.IsNaN(ai9.FloatAt(0)) || !math.IsNaN(bi9.FloatAt(0)) {
|
||
t.Fatalf("outside |x| ≤ 8 expected NaN, got %v and %v", ai9.FloatAt(0), bi9.FloatAt(0))
|
||
}
|
||
}
|