#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <set>
#include <map>

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JTools/JFunction1D_t.hh"
#include "JTools/JFunctionalMap_t.hh"
#include "JPhysics/JTransformer.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/Antares.hh"
#include "JPhysics/KM3NeT.hh"


/**
 * Scaling of absorption and scattering length.
 */
double  absorptionLengthFactor;
double  scatteringLengthFactor;


inline double getAbsorptionLength(const double lambda)
{
  return absorptionLengthFactor * NAMESPACE::getAbsorptionLength(lambda);
}


inline double getScatteringLength(const double lambda)
{
  return scatteringLengthFactor * NAMESPACE::getScatteringLength(lambda);
}


/**
 * Use attenutation.
 */
bool    useAttenuation;


/**
 * Get attenuation length.
 *
 * \param  lambda   wavelength of photon [nm]
 * \return          attenuation length [m]
 */
inline double getAttenuationLength(const double lambda)
{
  static const double f = 0.17;

  const double l_abs = getAbsorptionLength(lambda);
  const double ls    = getScatteringLength(lambda) / f;

  double y = 1.0 / l_abs;

  if (useAttenuation)
    y +=  1.0 / ls;

  return 1.0 / y;
}



/**
 * Main program to create interpolation tables for PDFs.
 */
int main(int argc, char **argv)
{
  using namespace std;

  string              outputFile;
  int                 numberOfPoints;
  double              epsilon;
  int                 function;
  int                 debug;
  
  try {

    JParser<> zap;

    zap['o'] = make_field(outputFile);
    zap['n'] = make_field(numberOfPoints)         = 25;
    zap['e'] = make_field(epsilon)                = 1.0e-10;
    zap['A'] = make_field(absorptionLengthFactor) = 1.0;
    zap['S'] = make_field(scatteringLengthFactor) = 1.0;
    zap['F'] = make_field(function)               = 1, 2, 3, 4;
    zap['U'] = make_field(useAttenuation);
    zap['d'] = make_field(debug)                  = 0;

    zap['F'] = JPARSER::not_initialised();

    if (zap.read(argc, argv) != 0)
      return 1;
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  using namespace JPHYSICS;
  using namespace JTOOLS;

  typedef double (JPDF::*fcn)(const double,
			      const double,
			      const double,
			      const double) const;


  // set global parameters

  const double P_atm = 240.0;           // ambient pressure [Atm]
  const double wmin  = 320.0;           // minimal wavelength [nm]
  const double wmax  = 610.0;           // maximal wavelength [nm]


  const JPDF_C
    pdf_c(NAMESPACE::getPhotocathodeArea(),
	  NAMESPACE::getQE,
	  NAMESPACE::getAngularAcceptance,
	  getAttenuationLength,
	  getScatteringLength,
	  NAMESPACE::p00075,
	  P_atm,
	  wmin,
	  wmax,
	  numberOfPoints,
	  epsilon);


  typedef JSplineFunction1D_t                                     JFunction1D_t;
  typedef JMultipleMap<3, JPolint1FunctionalMap>::type_list       JMapList_t;
  typedef JPDFTable<JFunction1D_t, JMapList_t>                    JPDF_t;

  typedef JFunctionTransformer<3, JFunction1D_t::argument_type>   JFunction3DTransformer_t;
  typedef JArray<3, JFunction1D_t::argument_type>                 JArray_t;

  JPDF_t pdf;


  NOTICE("building multi-dimensional function object <" << function << ">... " << flush);
  
  const double n[]  = { pdf_c.getIndexOfRefractionPhase(wmax),
			pdf_c.getIndexOfRefractionPhase(wmin) };
  
  const double ng[] = { pdf_c.getIndexOfRefractionGroup(wmax),
			pdf_c.getIndexOfRefractionGroup(wmin) };
  
  const double kmin = (ng[0] * n[0] - 1.0) / sqrt(n[0]*n[0] - 1.0);
  const double kmax = (ng[1] * n[1] - 1.0) / sqrt(n[1]*n[1] - 1.0);
  
  const double cmin = sqrt(ng[0]*ng[0] - 1.0);   // dt/dz = 0

  map<int, pair<fcn,JFunction3DTransformer_t> > zmap;

  zmap[1] = make_pair((fcn) &JPDF::getDirectLightFromMuon,         JFunction3DTransformer_t(21.5, 2, kmin, kmax, NAMESPACE::getAngularAcceptance, 0.001));
  zmap[2] = make_pair((fcn) &JPDF::getScatteredLightFromMuon,      JFunction3DTransformer_t(35.0, 2, cmin, 0.0,  NAMESPACE::getAngularAcceptance, 0.10));
  zmap[3] = make_pair((fcn) &JPDF::getDirectLightFromEMshowers,    JFunction3DTransformer_t(21.5, 2, cmin, 0.0,  NAMESPACE::getAngularAcceptance, 0.10));
  zmap[4] = make_pair((fcn) &JPDF::getScatteredLightFromEMshowers, JFunction3DTransformer_t(35.0, 2, cmin, 0.0,  NAMESPACE::getAngularAcceptance, 0.10));

  if (zmap.find(function) == zmap.end()) {
    FATAL("illegal function specifier" << endl);
  }

  fcn                      f           = zmap[function].first;    // PDF
  JFunction3DTransformer_t transformer = zmap[function].second;   // transformer


  set<double> R;        // distance [m]

  R.insert(  0.10);
  R.insert(  0.30);
  R.insert(  0.50);
  R.insert(  1.00);
  R.insert(  2.00);
  R.insert(  3.00);
  R.insert(  4.00);
  R.insert(  5.00);
  R.insert(  6.00);
  R.insert(  7.00);
  R.insert(  8.00);
  R.insert(  9.00);
  R.insert( 10.00);
  R.insert( 11.00);
  R.insert( 12.00);
  R.insert( 13.00);
  R.insert( 14.00);
  R.insert( 15.00);
  R.insert( 16.00);
  R.insert( 17.00);
  R.insert( 18.00);
  R.insert( 19.00);
  R.insert( 20.00);
  R.insert( 22.00);
  R.insert( 24.00);
  R.insert( 26.00);
  R.insert( 28.00);
  R.insert( 30.00);
  R.insert( 32.00);
  R.insert( 34.00);
  R.insert( 36.00);
  R.insert( 38.00);
  R.insert( 40.00);
  R.insert( 42.00);
  R.insert( 44.00);
  R.insert( 46.00);
  R.insert( 48.00);
  R.insert( 50.00);
  R.insert( 55.00);
  R.insert( 60.00);
  R.insert( 65.00);
  R.insert( 70.00);
  R.insert( 75.00);
  R.insert( 80.00);
  R.insert( 85.00);
  R.insert( 90.00);
  R.insert( 95.00);
  R.insert(100.00);
  R.insert(110.00);
  R.insert(120.00);
  R.insert(130.00);
  R.insert(140.00);
  R.insert(150.00);

  if (function == 3 || function == 4) {
    R.insert(170.00);
    R.insert(190.00);
    R.insert(210.00);
  }

  set<double> X;

  if (function == 1) {
 
    for (double buffer[] = { -0.01, -0.005, 0.0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, -1.0 }, *x = buffer; *x != -1.0; ++x) {
      X.insert(0.0 + *x);
      X.insert(1.0 - *x);
    }

    for (double x = 0.01; x < 0.1; x += 0.0025) {
      X.insert(0.0 + x);
      X.insert(1.0 - x);
    }

    for (double x = 0.10; x < 0.5; x += 0.010) {
      X.insert(0.0 + x);
      X.insert(1.0 - x);
    }
     
  } else {

    X.insert(  0.00);
    X.insert(  0.01);
    X.insert(  0.02);
    X.insert(  0.03);
    X.insert(  0.04);
    X.insert(  0.05);
    X.insert(  0.06);
    X.insert(  0.07);
    X.insert(  0.08);
    X.insert(  0.09);
    X.insert(  0.10);
    X.insert(  0.12);
    X.insert(  0.15);
    X.insert(  0.20);
    X.insert(  0.25);
    X.insert(  0.30);
    X.insert(  0.40);
    X.insert(  0.50);
    X.insert(  0.60);
    X.insert(  0.70);
    X.insert(  0.80);
    X.insert(  0.90);
    X.insert(  1.00);
    X.insert(  1.10);
    X.insert(  1.20);
    X.insert(  1.30);
    X.insert(  1.40);
    X.insert(  1.50);
    X.insert(  1.60);
    X.insert(  1.70);
    X.insert(  1.80);
    X.insert(  1.90);
    X.insert(  2.00);
    X.insert(  2.20);
    X.insert(  2.40);
    X.insert(  2.60);
    X.insert(  2.80);
    X.insert(  3.00);
    X.insert(  3.25);
    X.insert(  3.50);
    X.insert(  3.75);
    X.insert(  4.00);
    X.insert(  4.25);
    X.insert(  4.50);
    X.insert(  4.75);
    X.insert(  5.0);
    X.insert(  6.0);
    X.insert(  7.0);
    X.insert(  8.0);
    X.insert(  9.0);
    X.insert( 10.0);
    X.insert( 15.0);
    X.insert( 20.0);
    X.insert( 25.0);
    X.insert( 30.0);
    X.insert( 40.0);
    X.insert( 50.0);
    X.insert( 60.0);
    X.insert( 70.0);
    X.insert( 80.0);
    X.insert( 90.0);
    X.insert(100.0);
    X.insert(120.0);
    X.insert(140.0);
    X.insert(160.0);
    X.insert(180.0);
    X.insert(200.0);
    X.insert(250.0);
    X.insert(300.0);
    X.insert(350.0);
    X.insert(400.0);
    X.insert(450.0);
    X.insert(500.0);
    X.insert(600.0);
    X.insert(700.0);
    X.insert(800.0);
    X.insert(900.0);
    X.insert(999.0);
  }

    
  const double grid  =  5.0;                        // [deg]

  const double alpha = 2.0 * sqrt(1.0 - cos(grid * PI / 180.0));  // azimuth angle unit step size


  for (set<double>::const_iterator r = R.begin(); r != R.end(); ++r) {

    const double R_m = *r;

    const unsigned int number_of_theta_points = max(2u, (unsigned int) (180.0/(1.4 * grid)));
    
    for (double theta = 0.0; theta <= PI + epsilon; theta += PI/number_of_theta_points) {
      
      const unsigned int number_of_phi_points = max(2u, (unsigned int) (PI * sin(theta) / alpha));
      
      for (double phi = 0.0; phi <= PI + epsilon; phi += PI/number_of_phi_points) {
	
	JFunction1D_t& f1 = pdf[R_m][theta][phi];
	
	const JArray_t array(R_m, theta, phi);

	double t_old = transformer.getXn(array, *X.begin());
	double y_old = 0.0;

	for (set<double>::const_iterator x = X.begin(); x != X.end(); ++x) {

	  const double t = transformer.getXn(array, *x);
	  const double y = (pdf_c.*f)(R_m, theta, phi, t);

	  if (y != 0.0) {

	    if (*x < 0.0) {
	      WARNING("dt < 0 " << *x << ' ' << R_m << ' ' << t << ' ' << y << endl);
	    }

	    if (y_old == 0.0)
	      f1[t_old] = y_old;
              
	    f1[t] = y;
	    
	  } else {
            
	    if (y_old != 0.0)
	      f1[t] = y;
	  }
	  
	  t_old = t;
	  y_old = y;
        }
      }
    }
  }

  pdf.transform(transformer);
  pdf.compile();

  NOTICE("OK" << endl);

  try {

    NOTICE("storing output to file " << outputFile << "... " << flush);

    pdf.store(outputFile.c_str());

    NOTICE("OK" << endl);
  }
  catch(const JException& error) {
    FATAL(error.what() << endl);
  }
}