#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include <string>
#include <chrono>
#include <ctime>
#include <map>
#include <cmath>
#include <math.h>

#include "JAstronomy/JAstronomy.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"

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

/**
 * \file
 *                                                                                                                                                                    
 * Test program for JAstronomy.
 *
 * \author azegarelli
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace JASTRONOMY;

  double   utc_time_s;
  double   phi; 
  double   theta;
  string   detector;
  double   ra_j2000_input;
  double   dec_j2000_input;
  int      debug;

  try {

    JParser<> zap("Test program to test ra and dec j2000 coordinates.");

    zap['U'] = make_field(utc_time_s,           "UTC time in seconds");
    zap['P'] = make_field(phi,                  "phi angle in deg");
    zap['T'] = make_field(theta,                "theta angle in deg");
    zap['a'] = make_field(detector,             "detector");
    zap['R'] = make_field(ra_j2000_input,       "RAJ2000 input (known values)");
    zap['D'] = make_field(dec_j2000_input,      "DECJ2000 input (known values)");
    zap['d'] = make_field(debug)                = 3;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  const double t1 = utc_time_s + JASTRONOMY::MJD_EPOCH * 3600. * 24;  //number of seconds since MJD  
  
  // change theta's unit to degree;
  const double zenith =  180. - theta;               // flip zenith  for neutrino source direction (from detector to source)
  const double azimuth = fmod(phi + 180., 360.);     // flip azimuth for neutrino source direction (from detector to source)
  
  const double theta_rad  = theta  * M_PI/180;
  const double phi_rad = phi * M_PI/180;

  const JGEOMETRY3D::JAngle3D& neutrino_dir = JGEOMETRY3D::JAngle3D(theta_rad,phi_rad);

  int degreesNorth;
  int minutesNorth;
  int degreesEast;
  int minutesEast;

  if(detector=="ARCA") { degreesNorth = 36; minutesNorth = 17; degreesEast  = 15; minutesEast  = 58; }
  if(detector=="ORCA") { degreesNorth = 42; minutesNorth = 48; degreesEast  = 06; minutesEast  = 02; }

  const JGeographicalLocation &location = JGeographicalLocation(degreesNorth,minutesNorth,degreesEast,minutesEast);
  JAstronomy j_astro(location);

  const JASTRONOMY::JSourceLocation& source_location = j_astro.getLocationOfSourceFromZenithAzimuth(t1, neutrino_dir);

  double dec = (180/M_PI) * source_location.getDeclination();
  double ra  = (180/M_PI) * source_location.getRightAscension();
  double dec_j2000 = (180/M_PI) * source_location.getDeclinationJ2000(t1);
  double ra_j2000  = (180/M_PI) * source_location.getRightAscensionJ2000(t1);

  const JASTRONOMY::JGalacticCoordinates& gal_coordinates =  j_astro.getGalacticCoordinatesOfSource(t1, neutrino_dir);
  const double galactic_longitude = (180/M_PI) * gal_coordinates.getGalacticLongitude();
  const double galactic_latitude  = (180/M_PI) * gal_coordinates.getGalacticLatitude();

  DEBUG("Input (theta, phi) : (" << theta << ", " << phi << ")" <<endl);
  DEBUG("Output (zen, azi), (dec, ra, dec_j2000, ra_j2000), (gal_long, gal_lat): "  <<endl);
  DEBUG("( " << zenith << " " << azimuth << "),  (" << dec << ", " << ra << ", " << dec_j2000 << ", " << ra_j2000 << "), (" << galactic_longitude << ", " << galactic_latitude << ")" <<endl);

  // accept values if the difference between the known values of ra and dec j2000 and the computed ones is lower than 0.01 deg
  ASSERT(abs(ra_j2000 - ra_j2000_input)<=0.01, "Test whether RAJ2000 "<< ra_j2000 <<" differs from the known values ("<< ra_j2000_input<<") less for less than 0.01 deg --> diff = "<<abs(ra_j2000 -ra_j2000_input)<<" deg");
  ASSERT(abs(dec_j2000 - dec_j2000_input)<=0.01, "Test whether DECJ2000 "<< dec_j2000 <<" differs from the known values ("<< dec_j2000_input<<") less for less than 0.01 deg --> diff = "<<abs(dec_j2000 -dec_j2000_input)<<" deg");

}