// This is core/vnl/tests/test_bessel.cxx
#include <fstream>
#include "vnl/vnl_bessel.h"
#include "testlib/testlib_test.h"


static void
test_bessel()
{
  TEST_NEAR("J_0(0)=1.0", vnl_bessel(0, 0.0), 1.0, 1e-8);

  // First zero of Bessel function of order 0
  TEST_NEAR("J_0(2.4048)=0.0", vnl_bessel(0, 2.4048), 0.0, 1e-4);
  // Second zero of Bessel function of order 0
  TEST_NEAR("J_0(5.5201)=0.0", vnl_bessel(0, 5.5201), 0.0, 1e-4);

  TEST_NEAR("J_1(0)=0.0", vnl_bessel(1, 0.0), 0.0, 1e-5);

  // First zero of Bessel function of order 1
  TEST_NEAR("J_1(3.8315)=0.0", vnl_bessel(1, 3.8315), 0.0, 1e-4);

  TEST_NEAR("bessel0(1.23)", vnl_bessel0(1.23), vnl_bessel(0, 1.23), 1e-8);
  TEST_NEAR("bessel0(0.23)", vnl_bessel0(0.23), vnl_bessel(0, 0.23), 1e-8);
  TEST_NEAR("bessel0(0.001)", vnl_bessel0(0.001), 1.0, 1e-5);

  // Test consistency
  vnl_vector<double> J(3);
  vnl_bessel(2, 1.234, J);
  TEST_NEAR("bessel(0,1.234)", vnl_bessel(0, 1.234), J[0], 1e-8);
  TEST_NEAR("bessel(1,1.234)", vnl_bessel(1, 1.234), J[1], 1e-8);
  TEST_NEAR("bessel(2,1.234)", vnl_bessel(2, 1.234), J[2], 1e-8);
}

TESTMAIN(test_bessel);
