#include "Det.hh"

#include <map>
#include "foreach.hh"
#include "Table.hh"
#include "stringutil.hh"
#include "EventFile.hh"
#include <limits>
#include <streambuf>
using namespace stringutil;


/*! auxilliary function to deal with the fact that root ttimestamp uses ints and thus
    cannot deal with dates beyedon ~2038 */

void set_timestamp( double t, TTimeStamp& s )
{
  if ( t > std::numeric_limits<int>::max() )
  {
    s.SetSec( std::numeric_limits<int>::max() );
  }
  else
  {
    s.SetSec( int(t) );
    s.SetNanoSec( (t - int(t)) * 1e9);
  }
}



Det::Det(EventFile& evtfile)
{
  string detfile = split(string( evtfile.header["detector"]))[0];
  ::print ("constructing detector from", detfile );
  bool ok = read( detfile );
  if ( !ok) ::print ("FAILED");
}

int Det::get_utm_zone(double lat)
{
  return 1 + int( ( pi + lat ) / ( 6 * pi / 180 ) );
}

double Det::longitude_of_central_meridian(int utmzone)
{
  const double zone_width = 6 * pi / 180;
  return  -pi + (utmzone - 1) * zone_width + zone_width / 2;
}

vector<double> Det::compute_lonlat_from_utm( double utmX, double utmY, int zone )
{
  bool isNorthHemisphere = true; // todo

  double diflat = -0.00066286966871111111111111111111111111;
  double diflon = -0.0003868060578;

  double c_sa = 6378137.000000;
  double c_sb = 6356752.314245;
  double e2 = pow((pow(c_sa, 2) - pow(c_sb, 2)), 0.5) / c_sb;
  double e2cuadrada = pow(e2, 2);
  double c = pow(c_sa, 2) / c_sb;
  double x = utmX - 500000;
  double y = isNorthHemisphere ? utmY : utmY - 10000000;

  double s = ((zone * 6.0) - 183.0);
  double lat = y / (c_sa * 0.9996);
  double v = (c / pow(1 + (e2cuadrada * pow(cos(lat), 2)), 0.5)) * 0.9996;
  double a = x / v;
  double a1 = sin(2 * lat);
  double a2 = a1 * pow((cos(lat)), 2);
  double j2 = lat + (a1 / 2.0);
  double j4 = ((3 * j2) + a2) / 4.0;
  double j6 = ((5 * j4) + pow(a2 * (cos(lat)), 2)) / 3.0;
  double alfa = (3.0 / 4.0) * e2cuadrada;
  double beta = (5.0 / 3.0) * pow(alfa, 2);
  double gama = (35.0 / 27.0) * pow(alfa, 3);
  double bm = 0.9996 * c * (lat - alfa * j2 + beta * j4 - gama * j6);
  double b = (y - bm) / v;
  double epsi = ((e2cuadrada * pow(a, 2)) / 2.0) * pow((cos(lat)), 2);
  double eps = a * (1 - (epsi / 3.0));
  double nab = (b * (1 - epsi)) + lat;
  double senoheps = (exp(eps) - exp(-eps)) / 2.0;
  double delt  = atan(senoheps / (cos(nab) ) );
  double tao = atan(cos(delt) * tan(nab));

  double longitude = ((delt * (180.0 / pi)) + s) + diflon;
  double latitude = ((lat + (1 + e2cuadrada * pow(cos(lat), 2) - (3.0 / 2.0) * e2cuadrada * sin(lat) * cos(lat) * (tao - lat)) * (tao - lat)) * (180.0 / pi)) + diflat;

  return { longitude *pi / 180, latitude * pi / 180 };
}


double Det::compute_meridian_convergence_angle( double longitude,
    double latitude)
{
  double latitude_deg = latitude * deg;

  if ( latitude_deg > 84. || latitude_deg < -80.)
  {
    cout << "UTM coordinate system not defined for given Latitude: %f (max 84 deg N or min -80 deg S), setting meridian convergence angle to 0 deg.";
    return 0;
  }

  // detector position, longitude and latitude in rad
  double lambda  = longitude;
  double phi     = latitude;

  // find UTM zone and central meridian
  // longitude of the central meridian of UTM zone in rad

  double lambda0 = longitude_of_central_meridian( get_utm_zone( longitude ) );
  double omega   = lambda - lambda0;


  // parameters of the Earth ellipsoid
  double a  = 6378137.;      // semi-major axis in meters (WGS84)
  double e2 = 0.0066943800;  // eccentricity (WGS84)

  double rho = a * (1. - e2) / pow(1. - e2 * pow(sin(phi), 2), 3. / 2.);
  double nu  = a / sqrt(1 - e2 * pow(sin(phi), 2.));
  double psi = nu / rho;
  double t = tan(phi);

  // power series for convergence angle
  double gamma = sin(phi) * omega
                 - sin(phi) * pow(omega, 3.) / 3.  * pow(cos(phi), 2.) * (2.*pow(psi, 2.) - psi)
                 - sin(phi) * pow(omega, 5.) / 15. * pow(cos(phi), 4.) * (pow(psi, 4.) * (11. - 24.*pow(t, 2.)) -
                     pow(psi, 3.) * (11. - 36.*pow(t, 2.)) +
                     2.*pow(psi, 2.) * (1. - 7.*pow(t, 2.)) +
                     psi * pow(t, 2.)
                                                                         )
                 - sin(phi) * pow(omega, 7.) / 315. * pow(cos(phi), 6.) * (17. - 26.*pow(t, 2.) + 2.*pow(t, 4.));

  return gamma;
}


#include "io_util.hh"

/*! function to help with parsing comments in detx file */

string replace_newlines_in_quotes( string input, string quote = "\"", string rep ="<newline>")
{
	vector<string> v = split( input, quote );
	for (unsigned i = 1; i < v.size() ; i += 2 ) {
  	   v[i] = replace_all( v[i], "\n", rep  );
	}
  	
  	return join( quote, v );
}



bool read_detfile_detx( string filename, Det& det )
{
  // we have to allow for newlines that occur in quotes (" or ') in
  // comment lines. 
 
  string s = read_file( filename );
  string rep = "<newline>";
 
  s = replace_newlines_in_quotes( s, "\"", rep );
  s = replace_newlines_in_quotes( s, "'" , rep );

  vector<string> v = split(s,"\n");
  vector<string> vv;

  for( auto& x : v ) if ( !startswith( x, "#" ) ) vv.push_back( x );

  string ss = join( "\n", vv );
  ss = replace_all( ss, rep , "\n" );
  

  // now we can read the file...

  istringstream f(ss);

  Dom dom;
  Pmt pmt;

  int ndoms;
  string V;

  f >> det.id >> V;

  ::print( det.id, V );

  det.version = atoi( V.erase(0, 1).c_str() );

  if ( det.version > 1 )
  {
    double t0, t1;
    f >> t0 >> t1;

    set_timestamp( t0, det.validity_range_start );
    set_timestamp( t1, det.validity_range_stop  );

    string dum_utm;
    string utm_grid;

    f >>
      dum_utm                    >>
      det.utm_reference_elipsoid >>
      utm_grid                   >>
      det.utm_ref_easting        >>
      det.utm_ref_northing       >>
      det.utm_ref_z;

    if ( dum_utm != "UTM" || !f )
    {
      dbg("error reading detector file");
    }

    // utm_grid is of the form "33N", we need the 33 as int
    utm_grid.resize(  utm_grid.size() - 1 );
    // ::print ("utm_grid -> ", utm_grid );
    det.utm_zone = to<int> ( utm_grid );

    f >> ndoms;
  }
  else
  {
    ndoms = to<int>( V );
  }


  for (int idom = 0; idom < ndoms; idom++)
  {
    int npmts, dom_id;

    f >> dom_id ;
    Dom& dom = det.doms[ dom_id ];
    dom.id = dom_id;

    double q0, qx, qy, qz; // quaternions; ignore for now
    double module_t0;


    if ( det.version > 4 ) f >> dom.line_id >> dom.floor_id >> dom.pos >> q0 >> qx >> qy >> qz >> module_t0 >> dom.status >> npmts;
    else if ( det.version == 4 ) f >> dom.line_id >> dom.floor_id >> dom.pos >> q0 >> qx >> qy >> qz >> module_t0 >> npmts;
    else                         f >> dom.line_id >> dom.floor_id                                                 >> npmts;

    dom.pmts.resize( npmts );

    Vec cg;

    enumerate( channel, pmt, dom.pmts )
    {
      f >> pmt.id >> pmt.pos >> pmt.dir >> pmt.t0;
      if ( det.version >= 3)
      {
        f >> pmt.status;
      }

      pmt.dom_id         = dom.id;
      pmt.channel_id     = channel;
      det.pmts[ pmt.id ] = &pmt;
      cg += pmt.pos; // in older versions, we set the module pos to the cf of the pmts
    }
    if ( det.version < 4 ) dom.pos = cg/ dom.pmts.size() ; 
  }
  return true;
}


using namespace coolprinter;

bool read_detfile_det( string filename, Det& det )
{
  cout << filename << endl;
  map<string, vector<string> > M;

#define DEF( s ) M[ split(s)[0]] = slice(split(s),1);
  DEF("string            id x y z toffset twist twistrate");
  DEF("om                id pmttype serialnumber address");
  DEF("om_cluster        id string_id z n");
  DEF("om_position       id x y z theta phi");
  DEF("om_cluster_data   id x y z theta phi psi");
  DEF("lcm_logic         id string_id lcm_id n");
#undef DEF

  string s;
  ifstream f(filename.c_str());

  typedef map<string, string> Obj;      // object (=key value pairs)
  map<string, map<string, Obj > > objs; // map by tag and id

  while (getline(f, s ))
  {
    vector<string> v = split(s, ":");

    if (v.size() < 2 ) continue;

    string tag            = lower(v[0]);
    vector<string> values = split(v[1]);

    if ( M[tag].size() == 0 ) continue;

    Obj obj;

    foreach ( name_value, hzip( M[tag], values ) )
    {
      obj[ name_value[0] ] = name_value[1];
    }

    if ( tag == "lcm_logic")
    {
      obj["idlist"] = join(" ", slice( values, 4 ));
    }

    objs[tag][ obj["id"] ] = obj;
  }


  foreach_map( logicid, lcm_logic, objs["lcm_logic"] )
  {
    Obj& cluster    = objs["om_cluster"][ logicid ];
    Obj& line       = objs["string"][cluster["string_id"] ];

    int domid = to<int> ( lcm_logic["id"] );
    Dom& dom = det.doms[domid];

    dom.pmts.resize(31);
    dom.id       = domid;
    dom.line_id  = to<int> ( line["id"]      );
    dom.floor_id = 0;
    dom.pos.set( to<double>( line["x"]), to<double>( line["y"]), to<double>( line["z"] ) + to<double>(cluster["z"]));

    //int ichan = 0;
    string x = lcm_logic["idlist"];
    string xx = trim(x);

    foreach (id, split(xx) )
    {
      string pos_addr = objs["om"][ id ]["address"];
      //Obj& pos        = objs["om_position"][ pos_addr ];
      Obj& cluster    = objs["om_cluster"][ logicid ];
      Obj& line       = objs["string"][cluster["string_id"] ];

      int domid = to<int> ( lcm_logic["id"] );
      Dom& dom = det.doms[domid];

      dom.pmts.resize(31);
      dom.id       = domid;
      dom.line_id  = to<int> ( line["id"]      );
      dom.floor_id = 0;
      dom.pos.set( to<double>( line["x"]), to<double>( line["y"]), to<double>( line["z"] ) + to<double>(cluster["z"]));

      int ichan = 0;
      string x = lcm_logic["idlist"];
      string xx = stringutil::trim(x);

      foreach (id, split(xx) )
      {
        string pos_addr = objs["om"][ id ]["address"];
        Obj& pos        = objs["om_position"][ pos_addr ];
        Obj& cluster    = objs["om_cluster"][ logicid ];
        Obj& line       = objs["string"][cluster["string_id"] ];

        double x = to<double>( pos["x"] ) + to<double>( line["x"] );
        double y = to<double>( pos["y"] ) + to<double>( line["y"] );
        double z = to<double>( pos["z"] ) + to<double>( line["z"] ) + to<double>(cluster["z"]);

        // double theta = to<double> (pos["theta"]);
        // double phi   = to<double> (pos["phi"]);

        Pmt& pmt = dom.pmts[ichan];
        pmt.id         = to<int>(id);
        pmt.dom_id = dom.id;
        pmt.pos.set( x, y, z);
        pmt.dir.set_angles( to<double> (pos["theta"]), to<double> (pos["phi"]));
        pmt.t0 = 0;
        pmt.channel_id = ichan;

        ichan++;
      }
    }
  }

  det.init_pmts();
  return true;
}




#ifdef HAVEJPP

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JLang/JException.hh"
#include "Jeep/JMessage.hh"
#include "JDetector/JPMTRouter.hh"

// JDetector, JPMTRouter and load function are either in namespace
// JDETECTOR or JSIRENE depending on jpp version. The following
// uglyness makes it work for either.

namespace JDETECTOR {} // both must always
namespace JSIRENE {}   // exist
using namespace JDETECTOR; // let c++ figure out where it finds the
using namespace JSIRENE;   // required stuff.

void read_sirene( string detfilename, Det& det)
{
  JDetector jdet;

  try
  {
    cout << " loading " << endl;
    load( detfilename.c_str(), jdet );
  }
  catch (const JLANG::JException& error)
  {
    cout << "error" << endl;
    FATAL(error);
  }

  JPMTRouter router( jdet );

  foreach ( jdom, jdet )
  {
    int dom_id = jdom.getID ();

    foreach (jpmt, jdom)
    {
      Pmt pmt;// = det.pmts[pmt.id];
      pmt.id = jpmt.getID();
      pmt.pos.set( jpmt.getPosition().getX(),   jpmt.getPosition().getY(),   jpmt.getPosition().getZ()   );
      pmt.dir.set( jpmt.getDirection().getDX(), jpmt.getDirection().getDY(), jpmt.getDirection().getDZ() );

      int module_index = router.getAddress( pmt.id ).first;

      int id2 = jdet[ module_index ].getID ();

      cout << dom_id << " " << id2 << endl;
      cout << dom_id << " " << router.getAddress( pmt.id ).first << " " << router.getAddress( pmt.id ).second  << endl;
    }
  }
}
#else

void read_sirene( string detfilename, Det& det)
{
  cout << " read_sirene not available" << endl;
}

#endif


void read(Det& det, const JDAQSummaryslice& ss )
{
  det.zero_rates();
  foreach ( sf, ss ) // JDAQSummarysclice derives from vector<JDAQSummaryFrame>
  {
    Dom& dom = det.doms[ sf.getModuleID() ];
    foreach( pmt, dom.pmts ) pmt.rate = sf.getRate( pmt.channel_id );
  }
}