// 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))
	}
}