#include #include #include #include #include "TRandom3.h" #include "JLang/JException.hh" #include "JMath/JTrigonometric.hh" #include "JMath/JConstants.hh" #include "JTools/JFunction1D_t.hh" #include "JTools/JGrid.hh" #include "JTools/JQuantile.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { /** * Test method for interpolation with derivatives. * * Accuracy of nth derivative of N degree polynome * is equal to n-jth derivative of N-j degree polynome. * * \param argc number of command line arguments * \param argv list of command line arguments * \return status */ template int do_main(int argc, char **argv) { using namespace std; using namespace JPP; unsigned int numberOfEvents; unsigned int numberOfBins; int debug; try { JParser<> zap("Example program to test 1D interpolation of a polynome."); zap['n'] = make_field(numberOfEvents); zap['N'] = make_field(numberOfBins) = 21; zap['d'] = make_field(debug) = 3; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } if (N < 1) { FATAL("Number of degrees " << N << " < 1." << endl); } using namespace JPP; // Analytical function and its derivatives. vector f1(1, JTrigonometric(sin)); for (int i = 1; i != N + 1; ++i) { f1.push_back(f1.rbegin()->getDerivative()); } typedef JResultPolynome JResult_t; typedef JPolintFunction1D, JCollection, JResult_t> JFunction1D_t; JFunction1D_t polint; const double xmin = 0.0; const double xmax = PI; //const double dx = ((N+1)/2) * (xmax - xmin) / (double) (numberOfBins - 2*((N+1)/2)); const double dx = 0.0 * PI; const JGrid grid(numberOfBins + 1, xmin - dx, xmax + dx); polint.configure(grid, f1[0]); polint.compile(); polint.setExceptionHandler(new typename JFunction1D_t::JDefaultResult(JMATH::zero)); if (numberOfEvents != 0) { JQuantile Q[N+1]; for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) { ostringstream os; os << "f^" << i; Q[i].setTitle(os.str()); } for (unsigned int i = 0; i != numberOfEvents; ++i) { const double x = gRandom->Uniform(xmin, xmax); const JResult_t result = polint(x); for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) { Q[i].put(f1[i](x) - result.y[i]); } } cout << "\nInterpolation with " << N << " degrees polynomial:" << endl; for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) { Q[i].print(cout, i == 0); } } return 0; } } /** * \file * * Example program to test 1D interpolation of a trigonometric function. * \author mdejong */ int main(int argc, char **argv) { if (do_main<2>(argc, argv) != 0) { return 1; } if (do_main<3>(argc, argv) != 0) { return 1; } if (do_main<4>(argc, argv) != 0) { return 1; } if (do_main<5>(argc, argv) != 0) { return 1; } }