Add Bessel functions in order to use them in residual code. See "How can I use the Bessel function in the residual function?" at https://groups.google.com/d/msg/ceres-solver/Vh1gpqac8v0/NIK1EiWJCAAJ Change-Id: I3e80d9f9d1cadaf7177076e493ff46ace5233b76
diff --git a/include/ceres/jet.h b/include/ceres/jet.h index 0eb9b03..a21fd7a 100644 --- a/include/ceres/jet.h +++ b/include/ceres/jet.h
@@ -482,6 +482,41 @@ return Jet<T, N>(tanh_a, tmp * f.v); } +// Bessel functions of the first kind with integer order equal to 0, 1, n. +inline double BesselJ0(double x) { return j0(x); } +inline double BesselJ1(double x) { return j1(x); } +inline double BesselJn(int n, double x) { return jn(n, x); } + +// For the formulae of the derivatives of the Bessel functions see the book: +// Olver, Lozier, Boisvert, Clark, NIST Handbook of Mathematical Functions, +// Cambridge University Press 2010. +// +// Formulae are also available at http://dlmf.nist.gov + +// See formula http://dlmf.nist.gov/10.6#E3 +// j0(a + h) ~= j0(a) - j1(a) h +template <typename T, int N> inline +Jet<T, N> BesselJ0(const Jet<T, N>& f) { + return Jet<T, N>(BesselJ0(f.a), + -BesselJ1(f.a) * f.v); +} + +// See formula http://dlmf.nist.gov/10.6#E1 +// j1(a + h) ~= j1(a) + 0.5 ( j0(a) - j2(a) ) h +template <typename T, int N> inline +Jet<T, N> BesselJ1(const Jet<T, N>& f) { + return Jet<T, N>(BesselJ1(f.a), + T(0.5) * (BesselJ0(f.a) - BesselJn(2, f.a)) * f.v); +} + +// See formula http://dlmf.nist.gov/10.6#E1 +// j_n(a + h) ~= j_n(a) + 0.5 ( j_{n-1}(a) - j_{n+1}(a) ) h +template <typename T, int N> inline +Jet<T, N> BesselJn(int n, const Jet<T, N>& f) { + return Jet<T, N>(BesselJn(n, f.a), + T(0.5) * (BesselJn(n - 1, f.a) - BesselJn(n + 1, f.a)) * f.v); +} + // Jet Classification. It is not clear what the appropriate semantics are for // these classifications. This picks that IsFinite and isnormal are "all" // operations, i.e. all elements of the jet must be finite for the jet itself