JSphereND.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/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<int N>
  struct JVolume {
    /**
     * Get volume.
     *
     * \param  R       radius
     * \return         volume
     */
    static double get(const double R) 
    {
      return 2*JTOOLS::PI*R*R * JVolume<N-2>::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;
    }
  };
 
 
  /**
   * Get length squared.
   *
   * \param  key           key
   * \return               square of length
   */
  template<class JKey_t>
  inline double getLengthSquared(const JTOOLS::JMultiKey<1, JKey_t>& key)
  {
    return key.first*key.first;
  }
 
 
  /**
   * Get length squared.
   *
   * \param  key           key
   * \return               square of length
   */
  template<unsigned int N, class JKey_t>
  inline double getLengthSquared(const JTOOLS::JMultiKey<N, JKey_t>& key)
  {
    return key.first*key.first + getLengthSquared(key.second);
  }
}
 
 
/**
 * \file
 *
 * Example program to test integration of sphere in any number of dimensions.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  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);
  }
 
 
  const int    N = 7;
  const double R = 1.0;
 
  typedef JGridPolint1Function1D_t                                 JFunction1D_t;
  typedef JMultipleMap<N-2, JPolint1FunctionalGridMap>::typelist   JMaplist_t;
  typedef JMultiFunction<JFunction1D_t, JMaplist_t>                JMultiFunction_t;
  
 
  const JGrid<double> grid(numberOfBins, -R * 1.05, +R * 1.05);
 
  JMultiFunction_t gs;
 
  gs.configure(make_multigrid<N-2>(grid));
    
  for (JMultiFunction_t::super_iterator i = gs.super_begin(); i != gs.super_end(); ++i) {
 
    const double x = sqrt(getLengthSquared(i.getKey()));
 
    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<N>::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;
}