// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT

// Command pendulum computes the exact period of a simple pendulum at
// large amplitude through the complete elliptic integral of the first
// kind, and shows how far the small-angle formula drifts once the
// release angle stops being small. The period is
//
//	T = 4·sqrt(L/g)·K(sin²(θ₀/2)),
//
// where K is EllipticK with the m = k² parameter convention.
//
// Usage: go run ./examples/pendulum
package main

import (
	"fmt"
	"log"
	"math"

	"sourcedock.dev/petrbalvin/tensor"
)

// kComplete evaluates EllipticK at a single parameter.
func kComplete(m float64) float64 {
	arr, err := tensor.FromFloats([]float64{m}, 1)
	if err != nil {
		log.Fatal(err)
	}
	k, err := tensor.EllipticK(arr)
	if err != nil {
		log.Fatal(err)
	}
	v, _ := tensor.FloatAt(k, 0)
	return v
}

func main() {
	const (
		length = 1.0 // metres
		grav   = 9.80665
	)
	small := 2 * math.Pi * math.Sqrt(length/grav)

	fmt.Println("release angle    exact period    small-angle period    drift")
	for _, deg := range []float64{5, 15, 30, 45, 60, 90, 120, 170} {
		theta := deg * math.Pi / 180
		m := math.Sin(theta/2) * math.Sin(theta/2)
		period := 4 * math.Sqrt(length/grav) * kComplete(m)
		drift := (period/small - 1) * 100
		fmt.Printf("%10.0f°    %12.6f s    %14.6f s    %+6.2f %%\n",
			deg, period, small, drift)
	}

	// The inverse problem: which release angle doubles the small-angle
	// period? Bisection on the angle, the period being monotone in it.
	target := 2 * small
	lo, hi := 0.0, math.Pi
	angle := 0.0
	for range 80 {
		mid := (lo + hi) / 2
		m := math.Sin(mid/2) * math.Sin(mid/2)
		if 4*math.Sqrt(length/grav)*kComplete(m) < target {
			lo = mid
		} else {
			hi = mid
		}
		angle = mid
	}
	fmt.Printf("\na release angle of %.2f° doubles the period (%.4f s)\n",
		angle*180/math.Pi, target)
}