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


/**
 * \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<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 = (*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<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;
}