#ifndef __JASTRONOMYTOOLKIT__
#define __JASTRONOMYTOOLKIT__

#include <cmath>

#include "JAstronomy/JAstronomy.hh"

#include "JMath/JConstants.hh"
#include "JTools/JGrid.hh"
#include "JTools/JHistogram1D_t.hh"
#include "JTools/JFunction1D_t.hh"


/**
 * \author mdejong
 */

namespace JASTRONOMY {}
namespace JPP { using namespace JASTRONOMY; }

namespace JASTRONOMY {

  using JTOOLS::JGridPolint1Function1D_t;


  /**
   * Auxiliary class for source tracking.
   * This class constitutes a function object corresponding to the normalised probability
   *
   *     \f[ \frac{dP}{d\Omega}(cos\theta) \f] 
   *
   * for a source to be at an observed zenith angle, \f$ \theta \f$, in local coordinates.
   * A simulation of a diffuse flux of neutrinos can thus be used to determine 
   * the number of detected events from an astrophysical source with the given telescope.
   */
  struct JStarTrek :
    public JGridPolint1Function1D_t
  {
    /**
     * Constructor.
     *
     * \param  telescope            telescope
     * \param  source               astrophysical source
     * \param  number_of_bins       number of bins
     * \param  number_of_samples    number of samples
     */
    JStarTrek(const JAstronomy&      telescope,
	      const JSourceLocation& source,
	      const int              number_of_bins    =   401,
	      const int              number_of_samples = 20000)
    {
      using namespace JPP;

      const double epsilon = 0.005;
      const double W       = 1.0 / (double) (number_of_samples * 2*PI);   // dP/dphi

      JGridHistogram1D_t histogram(make_grid(number_of_bins, -1.0 - epsilon, +1.0 + epsilon));
      
      for (int i = 0; i != number_of_samples; ++i) {
	
	const double t1  = (double) i * NUMBER_OF_SECONDS_PER_SEDERIAL_DAY / (double) number_of_samples;	
	const double ct  = cos(telescope.getDirectionOfNeutrino(t1, source).getZenith());
	
	histogram.fill(ct, W);
      }
      
      makePDF(histogram, *this);                                          // dP/dOmega
      
      this->setExceptionHandler(new JGridPolint1Function1D_t::JDefaultResult(0.0));
      this->compile();
    }
  };
}

#endif