#include #include #include #include "JTools/JMapList.hh" #include "JTools/JFunction1D_t.hh" #include "JTools/JFunctionalMap_t.hh" #include "JTools/JMultiFunction.hh" #include "JTools/JMultipleMap.hh" #include "JTools/JMultiGrid.hh" #include "JTools/JToolsToolkit.hh" #include "JTools/JFunctionalMap.hh" #include "Jeep/JPrint.hh" #include "Jeep/JTimer.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { /** * Auxiliary class to determine volume of N-dimensional sphere. */ template struct JVolume { /** * Get volume. * * \param R radius * \return volume */ static double get(const double R) { return 2*JTOOLS::PI*R*R * JVolume::get(R) / N; } }; /** * Specialisation of template class JVolume for 2D sphere. */ template<> struct JVolume<2> { /** * Get volume. * * \param R radius * \return volume */ static double get(const double R) { return JTOOLS::PI*R*R; } }; /** * Specialisation of template class JVolume for 1D sphere. */ template<> struct JVolume<1> { /** * Get volume. * * \param R radius * \return volume */ static double get(const double R) { return 2*R; } }; } /** * \file * * Example program to test integration of sphere in any number of dimensions. * \author mdejong */ int main(int argc, char **argv) { using namespace std; unsigned int numberOfBins; double precision; int debug; try { JParser<> zap("Example program to test integration of sphere in any number of dimensions."); zap['N'] = make_field(numberOfBins) = 11; zap['e'] = make_field(precision) = 1.0e-2; zap['d'] = make_field(debug) = 2; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } using namespace JPP; const int N = 7; const double R = 1.0; typedef JGridPolint1Function1D_t JFunction1D_t; typedef JMultipleMap::typelist JMaplist_t; typedef JMultiFunction JMultiFunction_t; const JGrid grid(numberOfBins, -R * 1.05, +R * 1.05); JMultiFunction_t gs; gs.configure(make_multigrid(grid)); for (JMultiFunction_t::super_iterator i = gs.super_begin(); i != gs.super_end(); ++i) { const double x = (*i).getKey().getLength(); for (int __i = 0; __i != grid.getSize(); ++__i) { const double y = grid.getX(__i); const double z = R*R - x*x - y*y; if (z > 0.0) (*i).getValue()[y] = 2.0 * sqrt(z); else (*i).getValue()[y] = 0.0; } } gs.compile(); const double U = JVolume::get(R); JTimer timer("integrator"); timer.start(); const double W = getIntegral(gs); timer.stop(); NOTICE("Sphere " << N << "D" << endl); NOTICE("Volume (real) " << SCIENTIFIC(12,4) << U << endl); NOTICE("Volume (calc) " << SCIENTIFIC(12,4) << W << endl); NOTICE(timer); ASSERT(fabs(U - W) < precision * U); return 0; }