JTrigonometric1D.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <sstream>
#include <iomanip>
 
#include "TRandom3.h"
 
#include "JLang/JException.hh"
#include "JMath/JTrigonometric.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JGrid.hh"
#include "JTools/JQuantile.hh"
#include "JTools/JConstants.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
namespace {
 
  /**
   * Test method for interpolation with derivatives.
   *
   * Accuracy of  nth  derivative of  N  degree polynome
   * is equal to n-jth derivative of N-j degree polynome.
   */
  template<int N>
  int do_main(int argc, char **argv)
  {
    using namespace std;
    using namespace JPP;
 
    unsigned int  numberOfEvents;
    unsigned int  numberOfBins;
    int           debug;
 
    try {
 
      JParser<> zap("Example program to test 1D interpolation of a polynome.");
 
      zap['n'] = make_field(numberOfEvents);
      zap['N'] = make_field(numberOfBins)     = 21;
      zap['d'] = make_field(debug)            =  3;
 
    zap(argc, argv);
    }
    catch(const exception &error) {
      FATAL(error.what() << endl);
    }
 
    if (N < 1) {
      FATAL("Number of degrees " << N << " < 1." << endl);
    }
 
    
    using namespace JPP;
 
    // Analytical function and its derivatives.
    
    vector<JTrigonometric> f1(1, JTrigonometric(sin));
 
    for (int i = 1; i != N + 1; ++i) {
      f1.push_back(f1.rbegin()->getDerivative());
    }
 
 
    typedef JResultPolynome<N, double>                                                JResult_t;
    typedef JPolintFunction1D<N, JElement2D<double, double>, JCollection, JResult_t>  JFunction1D_t;
 
    JFunction1D_t  polint;
 
    const double xmin = 0.0;
    const double xmax = PI;
 
    //const double dx   = ((N+1)/2) * (xmax - xmin) / (double) (numberOfBins - 2*((N+1)/2));
    const double dx   = 0.0 * PI;
    
    const JGrid<double> grid(numberOfBins + 1, xmin - dx, xmax + dx);
 
    polint.configure(grid, f1[0]);
  
    polint.compile();
  
    polint.setExceptionHandler(new typename JFunction1D_t::JDefaultResult(JMATH::zero));
 
 
    if (numberOfEvents != 0) {
 
      JQuantile Q[N+1];
 
      for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
 
        ostringstream os;
 
        os << "f^" << i;
 
        Q[i].setTitle(os.str());
      }
 
      for (unsigned int i = 0; i != numberOfEvents; ++i) {
 
        const double    x      = gRandom->Uniform(xmin, xmax);
        const JResult_t result = polint(x);
 
        for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
          Q[i].put(f1[i](x) - result.y[i]);
        }
      }
 
      cout << "\nInterpolation with " << N << " degrees polynomial:" << endl;
 
      for (int i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
        Q[i].print(cout, i == 0);
      }
    }
    
    return 0;
  }
}
 
/**
 * \file
 *
 * Example program to test 1D interpolation of a trigonometric function.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  if (do_main<2>(argc, argv) != 0) { return 1; }
  if (do_main<3>(argc, argv) != 0) { return 1; }
  if (do_main<4>(argc, argv) != 0) { return 1; }
  if (do_main<5>(argc, argv) != 0) { return 1; }
}