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;
}