#include #include #include #include #include #include "JMath/JPolynome.hh" #include "JTools/JMapList.hh" #include "JTools/JFunction1D_t.hh" #include "JTools/JFunctionalMap_t.hh" #include "JTools/JMultiFunction.hh" #include "JTools/JToolsToolkit.hh" #include "JTools/JFunctionalMap.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { using namespace JPP; /** * 3D polynomial function. */ struct JFunction3D { /** * Result. */ struct result_type { /** * Constructor. * * \param __f function value f * \param __fx derivative df/dx * \param __fy derivative df/dy * \param __fz derivative df/dz */ result_type(const double __f, const double __fx, const double __fy, const double __fz) : f (__f), fx(__fx), fy(__fy), fz(__fz) {} /** * Type conversion operator. * * \return function value */ operator double() const { return f; } const double f; const double fx; const double fy; const double fz; }; /** * Constructor. * * In each dimension, the function value and the derivative are given by: *
     *      f (x) = f1.getValue(x);
     *      f'(x) = f1.getDerivative(x);
     * 
* * \param f1 polynome */ JFunction3D(const JPolynome& f1) : __f1(f1) {} /** * Evaluation of function value and derivatives. * * \param x x-value * \param y y-value * \param z z-value * \return function and derivatives */ inline result_type operator()(const double x, const double y, const double z) const { return result_type(f (x) * f (y) * f (z), fp(x) * f (y) * f (z), f (x) * fp(y) * f (z), f (x) * f (y) * fp(z)); } /** * Evaluation of integral. * * \param xmin minimal x-value * \param xmax maximal x-value * \param ymin minimal y-value * \param ymax maximal y-value * \param zmin minimal z-value * \param zmax maximal z-value * \return integral */ inline double operator()(const double xmin, const double xmax, const double ymin, const double ymax, const double zmin, const double zmax) const { return ((__f1.getIntegral(xmax) - __f1.getIntegral(xmin)) * (__f1.getIntegral(ymax) - __f1.getIntegral(ymin)) * (__f1.getIntegral(zmax) - __f1.getIntegral(zmin))); } protected: /** * Evaluation of function value. * * \param x x-value * \return function value */ inline double f(const double x) const { return __f1.getValue(x); } /** * Evaluation of derivative. * * \param x x-value * \return derivative */ inline double fp(const double x) const { return __f1.getDerivative(x); } JPolynome __f1; }; /** * Test method for multi-dimensional interpolation with derivatives. * * \param argc number of command line arguments * \param argv list of command line arguments * \param title title of this template process * \return status */ template int do_main(int argc, char **argv, const char* title) { using namespace std; int numberOfBins; JPolynome f1; map precision; int debug; precision["spline"] = 1.0e-3; precision["polint"] = 1.0e-10; try { JParser<> zap("Example program to test multi-dimensional integration."); zap['N'] = make_field(numberOfBins) = 11; zap['P'] = make_field(f1) = JPolynome(1.0, 1.0, 1.0); zap['e'] = make_field(precision) = JPARSER::initialised(); zap['d'] = make_field(debug) = 3; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } if (numberOfBins <= 2) { FATAL("Number of bins " << numberOfBins << endl); } using namespace JPP; const JFunction3D f3(f1); const double xmin = -1.0; const double xmax = +1.0; const double dx = (xmax - xmin) / (numberOfBins - 1); JMultiFunction_t g3; for (double x = xmin; x < xmax + 0.5*dx; x += dx) { for (double y = xmin; y < xmax + 0.5*dx; y += dx) { for (double z = xmin; z < xmax + 0.5*dx; z += dx) { g3[x][y][z] = f3(x,y,z); } } } g3.compile(); //g3.setExceptionHandler(new JFunction1D_t::JDefaultResult(JMATH::zero)); if (debug >= debug_t) { cout << endl << setw(8) << left << title << ' ' << setw(12) << right << "real" << ' ' << setw(12) << right << "calculated" << endl; cout << setw(8) << left << "integral"; cout << ' ' << SCIENTIFIC(12,4) << f3(xmin, xmax, xmin, xmax, xmin, xmax); cout << ' ' << SCIENTIFIC(12,4) << getIntegral(g3); cout << endl; } ASSERT(fabs(f3(xmin, xmax, xmin, xmax, xmin, xmax) - getIntegral(g3)) <= precision[title] * f3(xmin, xmax, xmin, xmax, xmin, xmax)); return 0; } } /** * \file * * Example program to test multi-dimensional integration. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; { typedef JGridPolint3Function1D_t JFunction1D_t; typedef JMAPLIST::maplist JMaplist_t; typedef JMultiFunction JMultiFunction_t; if (do_main(argc, argv, "polint") != 0) return 1; } { typedef JGridSplineFunction1D_t JFunction1D_t; typedef JMAPLIST::maplist JMaplist_t; typedef JMultiFunction JMultiFunction_t; if (do_main(argc, argv, "spline") != 0) return 1; } }