JSphere3D.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>

#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<JPolint3FunctionalGridMap>               JMaplist_t;  
  typedef JMultiFunction<JFunction1D_t, JMaplist_t>         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;
}