#ifndef __JAANET__JHONDAFLUXINTERPOLATOR__
#define __JAANET__JHONDAFLUXINTERPOLATOR__

#include <iostream>
#include <string>
#include <regex>

#include "Jeep/JPrint.hh"
#include "Jeep/JMessage.hh"

#include "JLang/JClonable.hh"
#include "JLang/JException.hh"

#include "JTools/JRange.hh"
#include "JTools/JPolint.hh"
#include "JTools/JMapList.hh"
#include "JTools/JCollection.hh"
#include "JTools/JGridCollection.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JMultiFunction.hh"
#include "JTools/JFunctionalMap_t.hh"

#include "JAAnet/JDiffuseFlux.hh"


/**
 * \author bjung
 */

namespace JAANET {}
namespace JPP { using namespace JAANET; }

namespace JAANET {

  using JEEP::JMessage;

  using JLANG::JClonable;

  using JTOOLS::JRange;
  using JTOOLS::JArray;
  
  using JTOOLS::JMAPLIST;
  
  using JTOOLS::JMultiFunction;
  using JTOOLS::JConstantFunction1D;
  using JTOOLS::JPolint1FunctionalMap;

  
  static const char* const HONDA_ENERGY_KEYWORD                         = "Enu(GeV)";
  static const char* const HONDA_COSINE_ZENITH_ANGLE_KEYWORD            = "cosZ";
  static const char* const HONDA_AZIMUTH_ANGLE_KEYWORD                  = "phi_Az";
  
  static const char* const HONDA_MUON_NEUTRINO_FLUX_KEYWORD             = "NuMu";
  static const char* const HONDA_MUON_ANTINEUTRINO_FLUX_KEYWORD         = "NuMubar";
  static const char* const HONDA_ELECTRON_NEUTRINO_FLUX_KEYWORD         = "NuE";
  static const char* const HONDA_ELECTRON_ANTINEUTRINO_FLUX_KEYWORD     = "NuEbar";

  static const char* const HONDA_BINSPECS_SEPARATOR                     = "--";
  static const char        HONDA_BINSPECS_BEGIN                         = '[';
  static const char        HONDA_BINSPECS_END                           = ']';
  

  /**
   * Auxiliary data structure for reading Honda bin ranges.
   */
  struct JHondaBinRange :
    JRange<double>
  {
    typedef JRange<double>        JRange_t;

	
    /**
     * Default constructor.
     */
    JHondaBinRange()
    {}


    /**
     * Constructor.
     *
     * \param  minimum            minimum value
     * \param  maximum            maximum value
     */
    JHondaBinRange(const double minimum,
		   const double maximum) :
      JRange_t(minimum, maximum)
    {}


    /**
     * Get central value of this range.
     *
     * \return                    central value
     */
    double getCentralValue() const
    {
      return (getLowerLimit() + getUpperLimit()) / 2.0;
    }

	
    /**
     * Read bin range from input.
     *
     * \param  in                 input stream
     * \param  object             bin range
     * \return                    input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JHondaBinRange& object)
    {
      using namespace std;
	
      double min = JRange_t::getMinimum();
      double max = JRange_t::getMaximum();

      string buffer;
	  
      if (in >> min >> buffer >> max &&
	  min    <= max              &&
	  buffer == HONDA_BINSPECS_SEPARATOR) {
	static_cast<JRange_t&>(object) = JRange_t(min, max);
      }

      return in;
    }
  };
      

  /**
   * Auxiliary data structure for reading angular binning specifications of Honda flux table.
   */
  struct JHondaAngularBinSpecs
  {
    /**
     * Default constructor.
     */
    JHondaAngularBinSpecs() :
      coszRange(-1.0,   1.0),
      phiRange ( 0.0, 360.0)
    {}


    /**
     * Constructor.
     *
     * \param  minCosz            minimum cosine zenith angle
     * \param  maxCosz            maximum cosine zenith angle
     * \param  minPhi             minimum azimuth angle
     * \param  maxPhi             maximum azimuth angle
     */
    JHondaAngularBinSpecs(const double minCosz,
			  const double maxCosz,
			  const double minPhi,
			  const double maxPhi) :
      coszRange(minCosz, maxCosz),
      phiRange (minPhi,  maxPhi)
    {}
    

    /**
     * Read bin specifications from input.
     *
     * \param  in                 input stream
     * \param  object             bin specifications
     * \return                    input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JHondaAngularBinSpecs& binspecs)
    {
      static const JEquationParameters eqpars("=", ",", "./", "#");
	  
      JProperties properties(eqpars, 1);
	  
      properties[HONDA_COSINE_ZENITH_ANGLE_KEYWORD] = binspecs.coszRange;
      properties[HONDA_AZIMUTH_ANGLE_KEYWORD]       = binspecs.phiRange;

      in >> properties;

      return in >> properties;
    }


    JHondaBinRange coszRange; //!< Cosine zenith-angle range
    JHondaBinRange phiRange;  //!< Azimuthal angle range [deg]
  };
    
  
  /** 
   * Template definition for Honda flux table interpolator.
   */
  template<class                               JPhiFunction_t         = JConstantFunction1D<double, JArray<4, double> >,
	   template<class, class, class> class JCoszFunctionalMap_t   = JPolint1FunctionalMap,
	   template<class, class, class> class JEnergyFunctionalMap_t = JPolint1FunctionalMap>
  class JHondaFluxInterpolator final :
    public JMultiFunction                          <JPhiFunction_t, typename JMAPLIST<JEnergyFunctionalMap_t, JCoszFunctionalMap_t>::maplist>,
    public JClonable<JFlux, JHondaFluxInterpolator <JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >,
    public JMessage        <JHondaFluxInterpolator <JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >
  {
  public:

    typedef typename JMAPLIST<JEnergyFunctionalMap_t, JCoszFunctionalMap_t>::maplist                                      JFunctionalMaplist_t;
    
    typedef JHondaFluxInterpolator<JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t>                             interpolator_type;
    typedef JMultiFunction<JPhiFunction_t, JFunctionalMaplist_t>                                                            multifunction_type;
    
    enum { NUMBER_OF_DIMENSIONS = multifunction_type::NUMBER_OF_DIMENSIONS };

    typedef typename multifunction_type::abscissa_type                                                                           abscissa_type;
    typedef typename multifunction_type::argument_type                                                                           argument_type;
    typedef typename multifunction_type::result_type                                                                               result_type;
    typedef typename multifunction_type::value_type                                                                                 value_type;

    typedef typename multifunction_type::multimap_type                                                                           multimap_type;
    
    typedef typename multifunction_type::super_const_iterator                                                             super_const_iterator;
    typedef typename multifunction_type::super_iterator                                                                         super_iterator;
    typedef typename multifunction_type::function_type                                                                           function_type;

    typedef JRange<abscissa_type>                                                                                                     JRange_t;
    
    
    /**
     * Default constructor.
     */
    JHondaFluxInterpolator() :
      table(),
      coszRange(0.0, -1.0),
      phiRange (0.0, -1.0)
    {}


    /**
     * Constructor.
     *
     * \param  fileName             oscillation probability table filename
     */
    JHondaFluxInterpolator(const char* fileName) :
      JHondaFluxInterpolator()
    {
      this->load(fileName);
      this->check_validity();
    }
    

    /**
     * Load oscillation probability table from file.
     *
     * \param  fileName             oscillation probability table filename
     */
    void load(const char* const fileName)
    {
      using namespace std;
      using namespace JPP;

      try {
	
	NOTICE("Loading Honda flux table from file " << fileName << "... " << flush);

	multimap_type& zmap = static_cast<multimap_type&>(*this);

	ifstream ifs(fileName);
      
	JHondaAngularBinSpecs binspecs;
	
	for (string line; getline(ifs, line); ) {

	  double energy = 0.0;
	
	  result_type values;

	  istringstream iss1(line);

	  const double cosz = binspecs.coszRange.getCentralValue();
	  const double phi  = binspecs.phiRange .getCentralValue();

	  coszRange.include(cosz);
	  phiRange .include(phi);
	  
	  if (iss1 >> energy >> values) {
	    
	    const double log10E = log10(energy);

	    for (typename result_type::iterator i = values.begin(); i != values.end(); ++i) {
	      *i = log(*i);
	    }

	    zmap[log10E][cosz][phi] = values;
	    
	  } else {

	    const size_t p1 = line.find (HONDA_BINSPECS_BEGIN);
	    const size_t p2 = line.rfind(HONDA_BINSPECS_END);
	
	    if (p1 != string::npos && p2 != string::npos) {

	      const string binspecstr = line.substr(p1 + 1, p2 - p1 - 1);
	    
	      istringstream iss2(binspecstr);
	    
	      iss2 >> binspecs;
	    }
	  }
	}
      }
      catch(const exception& error) {
	THROW(JFileReadException, "JHondaFluxInterpolator::load(): Error reading file " << fileName);	
      }
    }


    /**
     * Check whether this interpolator is valid.
     *
     * \return                   true if valid; else false
     */
    bool is_valid() const override final
    {
      return this->size() > 0 && coszRange.getLength() > 0.0 && phiRange.getLength() > 0.0;
    }

    
    /**
     * Get diffuse flux for given particle PDG-identifier, energy, cosine zenith angle and azimuth angle.
     *
     * \param  type              PDG particle type
     * \param  log10E            logarithmic neutrino energy [GeV]
     * \param  costh             cosine zenith angle
     * \param  phi               azimuth angle (defined clockwise w.r.t. the North) [deg]
     * \return                   diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double getFactor(const int    type,
		     const double log10E,
		     const double costh,
		     const double phi) const
    {
      using namespace std;
      using namespace JPP;

      static const double epsilon_costh = 1.0e-3; //!< [--]
      static const double epsilon_phi   = 1.0;    //!< [deg]

      static const JRange_t mod(0.0, 360.0);
      
      const multifunction_type& function = static_cast<const multifunction_type&>(*this);

      const double _costh = min(max(costh, -1.0), 1.0);
      const double _phi   = mod.mod(180.0 - phi); //!< Honda azimuth angle is defined counterclockwise w.r.t. the South\n
                                                  //!< (c.f. arXiv:1102.2688), so a conversion is necessary.
      
      result_type lnFluxes;
      
      if (coszRange(_costh) && phiRange(_phi)) {
	
	lnFluxes = function(log10E, _costh, _phi);
	
      } else { // Perform linear extrapolation for edges of angular ranges
	
	if        (_costh > coszRange.getUpperLimit() && phiRange(_phi)) {

	  const abscissa_type x1 = coszRange.getUpperLimit();
	  const abscissa_type x2 = coszRange.getUpperLimit() - epsilon_costh;
	    
	  const result_type&  y1 = function(log10E, x1, _phi);	
	  const result_type&  y2 = function(log10E, x2, _phi);

	  lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1);
	  
	} else if (_costh < coszRange.getLowerLimit() && phiRange(_phi)) {
	
	  const abscissa_type x1 = coszRange.getLowerLimit();
	  const abscissa_type x2 = coszRange.getLowerLimit() + epsilon_costh;
	    
	  const result_type&  y1 = function(log10E, x1, _phi);	
	  const result_type&  y2 = function(log10E, x2, _phi);

	  lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1);
	
	} else if (coszRange(_costh) && _phi > phiRange.getUpperLimit()) {

	  const abscissa_type x1 = phiRange.getUpperLimit();
	  const abscissa_type x2 = phiRange.getUpperLimit() - epsilon_costh;
	    
	  const result_type&  y1 = function(log10E, _costh, x1);	
	  const result_type&  y2 = function(log10E, _costh, x2);

	  lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_phi - x1);
	
	} else if (coszRange(_costh) && _phi < phiRange.getLowerLimit()) {

	  const abscissa_type x1 = phiRange.getLowerLimit();
	  const abscissa_type x2 = phiRange.getLowerLimit() + epsilon_costh;
	    
	  const result_type&  y1 = function(log10E, _costh, x1);	
	  const result_type&  y2 = function(log10E, _costh, x2);

	  lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_phi - x1);
	
	} else if (_costh > coszRange.getUpperLimit() && _phi > phiRange.getUpperLimit()) {

	  const abscissa_type c1 = coszRange.getUpperLimit();
	  const abscissa_type c2 = coszRange.getUpperLimit() - epsilon_costh;
	  
	  const abscissa_type p1 = phiRange .getUpperLimit();
	  const abscissa_type p2 = phiRange .getUpperLimit() - epsilon_phi;

	  const result_type&  y1 = function(log10E, c1, p1);
	  const result_type&  y2 = function(log10E, c2, p2);

	  lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1);
	  
	} else if (_costh > coszRange.getUpperLimit() && _phi < phiRange.getUpperLimit()) {

	  const abscissa_type c1 = coszRange.getUpperLimit();
	  const abscissa_type c2 = coszRange.getUpperLimit() - epsilon_costh;
	  
	  const abscissa_type p1 = phiRange .getLowerLimit();
	  const abscissa_type p2 = phiRange .getLowerLimit() + epsilon_phi;

	  const result_type&  y1 = function(log10E, c1, p1);
	  const result_type&  y2 = function(log10E, c2, p2);

	  lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1);
	  
	} else if (_costh < coszRange.getLowerLimit() && phi > phiRange.getUpperLimit()) {

	  const abscissa_type c1 = coszRange.getLowerLimit();
	  const abscissa_type c2 = coszRange.getLowerLimit() + epsilon_costh;
	  
	  const abscissa_type p1 = phiRange .getUpperLimit();
	  const abscissa_type p2 = phiRange .getUpperLimit() - epsilon_phi;

	  const result_type&  y1 = function(log10E, c1, p1);
	  const result_type&  y2 = function(log10E, c2, p2);

	  lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1);
	  
	} else {

	  const abscissa_type c1 = coszRange.getLowerLimit();
	  const abscissa_type c2 = coszRange.getLowerLimit() + epsilon_costh;
	  
	  const abscissa_type p1 = phiRange .getLowerLimit();
	  const abscissa_type p2 = phiRange .getLowerLimit() + epsilon_phi;

	  const result_type&  y1 = function(log10E, c1, p1);
	  const result_type&  y2 = function(log10E, c2, p2);

	  lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1);
	}
      }

      switch(type) {
      case TRACK_TYPE_NUMU:
	return exp(lnFluxes[0]);
      case TRACK_TYPE_ANTINUMU:
	return exp(lnFluxes[1]);
      case TRACK_TYPE_NUE:
	return exp(lnFluxes[2]);
      case TRACK_TYPE_ANTINUE:
	return exp(lnFluxes[3]);
      default:
	return 0.0;
      }
    }


    /**
     * Get flux of given event.
     *
     * \param  type              PDG particle type
     * \param  log10E            logarithmic neutrino energy [GeV]
     * \param  costh             cosine zenith angle
     * \param  phi               azimuth angle (defined clockwise w.r.t. the North) [deg]
     * \return                   diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double getFlux(const int    type,
		   const double log10E,
		   const double costh,
		   const double phi) const
    {
      return getFactor(type, log10E, costh, phi);
    }


    /**
     * Get diffuse flux for given particle PDG-identifier, energy, cosine zenith angle and azimuth angle.
     *
     * \param  type              PDG particle type
     * \param  log10E            logarithmic neutrino energy [GeV]
     * \param  costh             cosine zenith angle
     * \param  phi               azimuth angle (defined clockwise w.r.t. the North) [deg]
     * \return                   diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double operator()(const int    type,
		      const double log10E,
		      const double costh,
		      const double phi) const
    {
      return getFactor(type, log10E, costh, phi);
    }


    /**
     * Get flux of given event.
     *
     * \param  evt               event
     * \return                   flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double getFactor(const Evt& evt) const override final
    {
      const Trk& neutrino = get_neutrino(evt);
      
      const double log10E = log10(neutrino.E);
      const double costh  = neutrino.dir.z;
      
      //!< The gSeaGen azimuth angle is defined clockwise w.r.t. the North (c.f. arXiv:2003.14040)
      const double phi = atan(neutrino.dir.y / neutrino.dir.x) * 180.0 / M_PI;

      return getFactor(neutrino.type, log10E, costh, phi);
    }


    /**
     * Get flux of given event.
     *
     * \param  event             event
     * \return                   flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double getFlux(const Evt& event) const
    {
      return getFactor(event);
    }


    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) override final
    {
      return JHondaFluxInterpolatorHelper(*this, eqpars);
    }


    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) const override final
    {
      return JHondaFluxInterpolatorHelper(*this, eqpars);
    }


    /**
     * Stream input.
     *
     * \param  in              input stream
     * \return                 input stream
     */
    std::istream& read(std::istream& in) override final
    {
      using namespace std;

      streampos pos = in.tellg();      

      if (!(in >> table)) {

	in.clear();
	in.seekg(pos);
	
	JProperties properties = getProperties();
	
	in >> properties;
      }

      load(table.c_str());
      
      this->check_validity();
      
      return in;
    }


  private:
    
    /**
     * Auxiliary class for I/O of Honda flux interpolator.
     */
    struct JHondaFluxInterpolatorHelper :
      public JProperties
    {
      /**
       * Constructor.
       *
       * \param  interpolator    Honda flux interpolator
       * \param  eqpars          equation parameters 
       */
      template<class JHondaFluxInterpolator_t>
      JHondaFluxInterpolatorHelper(const JHondaFluxInterpolator_t& interpolator,
				   const JEquationParameters&      eqpars) :
	JProperties(eqpars, 1)
      {
	(*this)[JEvtWeightFactor::getTypeKey()] = "Honda flux interpolator";

	this->insert(gmake_property(interpolator.table));
      }
    };

    
    std::string table;  //!< Table filename

    JRange_t coszRange; //!< Zenith  angle range
    JRange_t phiRange;  //!< Azimuth angle range [deg]
    
    using JMessage<interpolator_type>::debug; //!< Debug level
  };



  /** 
   * Template specialisation for interpolator of azimuth-averaged Honda flux table.
   */
  template<template<class, class, class> class JCoszFunctionalMap_t,
	   template<class, class, class> class JEnergyFunctionalMap_t>
  class JHondaFluxInterpolator<JConstantFunction1D<double, JArray<4, double> >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> final :
    public JMultiFunction                                 <JConstantFunction1D<double, JArray<4, double> >, typename JMAPLIST<JEnergyFunctionalMap_t, JCoszFunctionalMap_t>::maplist>,
    public JClonable<JDiffuseFlux, JHondaFluxInterpolator <JConstantFunction1D<double, JArray<4, double> >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >,
    public JMessage               <JHondaFluxInterpolator <JConstantFunction1D<double, JArray<4, double> >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >
  {
  public:
    
    typedef typename JMAPLIST<JEnergyFunctionalMap_t, JCoszFunctionalMap_t>::maplist                                      JFunctionalMaplist_t;
    
    typedef JConstantFunction1D<double, JArray<4, double> >                                                                     JPhiFunction_t;
    
    typedef JHondaFluxInterpolator<JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t>                             interpolator_type;
    typedef JMultiFunction<JPhiFunction_t, JFunctionalMaplist_t>                                                            multifunction_type;
    
    enum { NUMBER_OF_DIMENSIONS = multifunction_type::NUMBER_OF_DIMENSIONS };

    typedef typename multifunction_type::abscissa_type                                                                           abscissa_type;
    typedef typename multifunction_type::argument_type                                                                           argument_type;
    typedef typename multifunction_type::result_type                                                                               result_type;
    typedef typename multifunction_type::value_type                                                                                 value_type;

    typedef typename multifunction_type::multimap_type                                                                           multimap_type;

    typedef typename multifunction_type::super_const_iterator                                                             super_const_iterator;
    typedef typename multifunction_type::super_iterator                                                                         super_iterator;
    typedef typename multifunction_type::function_type                                                                           function_type;

    typedef JRange<abscissa_type>                                                                                                     JRange_t;
    
    
    /**
     * Default constructor.
     */
    JHondaFluxInterpolator() :
      table(),
      coszRange(0.0, -1.0)
    {}


    /**
     * Constructor.
     *
     * \param  fileName             oscillation probability table filename
     */
    JHondaFluxInterpolator(const char* fileName) :
      JHondaFluxInterpolator()
    {
      this->load(fileName);
      this->check_validity();
    }
    

    /**
     * Load oscillation probability table from file.
     *
     * \param  fileName             oscillation probability table filename
     */
    void load(const char* const fileName)
    {
      using namespace std;
      using namespace JPP;

      try {
	
	NOTICE("Loading Honda flux table from file " << fileName << "... " << flush);	

	multimap_type& zmap = static_cast<multimap_type&>(*this);

	ifstream ifs(fileName);
	
	JHondaAngularBinSpecs binspecs;
      
	for (string line; getline(ifs, line); ) {

	  double energy = 0.0;
	
	  result_type values;

	  istringstream iss1(line);

	  const double cosz = binspecs.coszRange.getCentralValue();
	  coszRange.include(cosz);

	  if (iss1 >> energy >> values) {

	    const double log10E = log10(energy);
	  
	    for (typename result_type::iterator i = values.begin(); i != values.end(); ++i) {
	      *i = log(*i);
	    }

	    zmap[log10E][cosz] = values;
	    
	  } else {

	    const size_t p1 = line.find (HONDA_BINSPECS_BEGIN);
	    const size_t p2 = line.rfind(HONDA_BINSPECS_END);
	
	    if (p1 != string::npos && p2 != string::npos) {

	      const string binspecstr = line.substr(p1 + 1, p2 - p1 - 1);
	    
	      istringstream iss2(binspecstr);
	    
	      iss2 >> binspecs;
	    }
	  }
	}
	
	NOTICE("OK" << endl);	
      }
      catch(const exception& error) {
	THROW(JFileReadException, "JHondaFluxInterpolator::load(): Error reading file " << fileName);	
      }
    }


    /**
     * Check whether this interpolator is valid.
     *
     * \return                   true if valid; else false
     */
    bool is_valid() const override final
    {
      return this->size() > 0 && coszRange.getLength();
    }
    
    
    /**
     * Get diffuse flux for given particle PDG-identifier, energy and zenith-angle.
     *
     * \param  type              PDG particle type
     * \param  log10E            logarithmic neutrino energy [GeV]
     * \param  costh             cosine zenith angle
     * \return                   diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double dNdEdOmega(const int    type,
		      const double log10E,
		      const double costh) const override final
    {
      using namespace std;
      using namespace JPP;

      static const double epsilon_costh = 1.0e-3;
      
      const multifunction_type& function = static_cast<const multifunction_type&>(*this);

      const double _costh = min(max(costh, -1.0), 1.0);

      result_type lnFluxes;
      
      if (coszRange(_costh)) {
	
	lnFluxes = function(log10E, _costh);
	
      } else { // Perform linear extrapolation for edges of cosine zenith angle range
	
	if        (_costh > coszRange.getUpperLimit()) {

	  const abscissa_type x1 = coszRange.getUpperLimit();
	  const abscissa_type x2 = coszRange.getUpperLimit() - epsilon_costh;
	    
	  const result_type&  y1 = function(log10E, x1);	
	  const result_type&  y2 = function(log10E, x2);

	  lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1);
	  
	} else {
	
	  const abscissa_type x1 = coszRange.getLowerLimit();
	  const abscissa_type x2 = coszRange.getLowerLimit() + epsilon_costh;
	    
	  const result_type&  y1 = function(log10E, x1);	
	  const result_type&  y2 = function(log10E, x2);

	  lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1);
	}
      }

      switch(type) {
      case TRACK_TYPE_NUMU:
	return exp(lnFluxes[0]);
      case TRACK_TYPE_ANTINUMU:
	return exp(lnFluxes[1]);
      case TRACK_TYPE_NUE:
	return exp(lnFluxes[2]);
      case TRACK_TYPE_ANTINUE:
	return exp(lnFluxes[3]);
      default:
	return 0.0;
      }
    }


    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) override final
    {
      return JHondaFluxInterpolatorHelper(*this, eqpars);
    }


    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) const override final
    {
      return JHondaFluxInterpolatorHelper(*this, eqpars);
    }
    
    
    /**
     * Stream input.
     *
     * \param  in              input stream
     * \return                 input stream
     */
    std::istream& read(std::istream& in) override final
    {
      using namespace std;

      streampos pos = in.tellg();      

      if (!(in >> table)) {

	in.clear();
	in.seekg(pos);
	
	JProperties properties = getProperties();
	
	in >> properties;
      }

      load(table.c_str());
      
      this->check_validity();
      
      return in;
    }

  private:
    
    /**
     * Auxiliary class for I/O of Honda flux interpolator.
     */
    struct JHondaFluxInterpolatorHelper :
      public JProperties
    {
      /**
       * Constructor.
       *
       * \param  interpolator    Honda flux interpolator
       * \param  eqpars          equation parameters 
       */
      template<class JHondaFluxInterpolator_t>
      JHondaFluxInterpolatorHelper(const JHondaFluxInterpolator_t& interpolator,
				   const JEquationParameters&      eqpars) :
	JProperties(eqpars, 1)
      {
	(*this)[JEvtWeightFactor::getTypeKey()] = "Honda flux interpolator";

	this->insert(gmake_property(interpolator.table));
      }
    };

    
    std::string table;  //!< Table filename
    
    JRange_t coszRange; //!< Cosine zenith angle range
    
    using JMessage<interpolator_type>::debug; //!< Debug level
  };


  /**
   * Alias template definition for 2-dimensional Honda flux interpolator.
   */
  template<template<class, class, class> class JCoszFunctionalMap_t   = JPolint1FunctionalMap,
	   template<class, class, class> class JEnergyFunctionalMap_t = JPolint1FunctionalMap>
  using JHondaFluxInterpolator2D = JHondaFluxInterpolator<JConstantFunction1D<double, JArray<4, double> >,
							  JCoszFunctionalMap_t,
							  JEnergyFunctionalMap_t>;
}


namespace JEEP {

  /**
   * JMessage template specialization for oscillation probability interpolators.
   */
  template<class                               JPhiFunction_t,
	   template<class, class, class> class JCoszFunctionalMap_t,
	   template<class, class, class> class JEnergyFunctionalMap_t>  
  struct JMessage<JAANET::JHondaFluxInterpolator<JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >
  {
    static int debug;
  };


  /**
   * Default verbosity for oscillation probability interpolators.
   */
  template<class                               JPhiFunction_t,
	   template<class, class, class> class JCoszFunctionalMap_t,
	   template<class, class, class> class JEnergyFunctionalMap_t>  
  int JMessage<JAANET::JHondaFluxInterpolator<JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >::debug = (int) notice_t;
  
}

#endif