#include #include #include #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" /** * \file * * Example program to test 3D integration of sphere. * \author mdejong */ int main(int argc, char **argv) { using namespace std; unsigned int numberOfBins; int debug; try { JParser<> zap("Example program to test 3D integration of sphere."); zap['N'] = make_field(numberOfBins) = 101; zap['d'] = make_field(debug) = 3; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } using namespace JPP; const double R = 1.0; typedef JGridPolint3Function1D_t JFunction1D_t; typedef JMapList JMaplist_t; typedef JMultiFunction JMultiFunction_t; JMultiFunction_t ga; JMultiFunction_t gb; const double xmin = -R * 1.05; const double xmax = +R * 1.05; const double dx = (xmax - xmin) / (numberOfBins - 1); for (double x = xmin; x < xmax + 0.5*dx; x += dx) { for (double y = xmin; y < xmax + 0.5*dx; y += dx) { const double z = R*R - x*x - y*y; if (z > 0.0) { ga[x][y] = -sqrt(z); gb[x][y] = +sqrt(z); } else { ga[x][y] = 0.0; gb[x][y] = 0.0; } } } ga.compile(); gb.compile(); cout << "Volume (real) " << SCIENTIFIC(12,4) << 4.0*PI*R*R*R/3.0 << endl; cout << "Volume (calc) " << SCIENTIFIC(12,4) << getIntegral(gb) - getIntegral(ga) << endl; }