#ifndef ASTROHHINCD
#define ASTROHHINCD

class Det;
class Vec;
class TTimeStamp;
#include <vector>
#include <string>
#include "Mat.hh"
#include "constants.hh"



/* bring angle between 0 and 2 pi */
inline double wrap( double angle )
{
  int n = int(angle / (2 * pi)) - (angle < 0);
  angle -= n * 2 * pi;
  return angle;
}


inline double siderial_day() {
  return 86164.1; // seconds
}





/*! Equatorial Coordinates, right assencion and declination */
struct EquatorialCoords
{
  double _ra, _dec;

  EquatorialCoords() {}
  EquatorialCoords(double ra, double dec) : _ra(ra), _dec(dec) {}


  /*! construct the equatorial coordinates. The Det object holds the geogrphical
     information of the detector. Timestamp must be in UTC */

  EquatorialCoords( const Det& det,
                    const Vec& track_direction,
                    const TTimeStamp& t,
                    bool j2000 = true );


  /*! set right ascension and declination (both in radians) */

  void set( double ra, double dec) { _ra = ra; _dec = dec ;}


  /*! set right ascension. hours should be between 0 and 23, minutes
    between 0 and 59 etc.  */

  void set_ra( unsigned hours, unsigned minutes, double seconds )
  {
    _ra = hours + minutes / 60.0 + seconds / 3600.0;
    _ra *= 15 / deg;
  }

  /*! Set using degrees, arc-minutes and arc-seconds. Note you need to
      specify the sign seperately. */
  void set_dec( unsigned degrees, unsigned arcminutes, double arcseconds, int sign = 1 )
  {
    _dec = degrees + arcminutes / 60.0 + arcseconds / 3600.0;
    _dec *= sign / deg;
  }

  void set_dec( double radians) { _dec = radians; };
  void set_ra( double radians) { _ra = radians; };

  double ra()     const { return _ra;}
  double ra_deg() const { return _ra * deg;}
  double dec()    const { return _dec;}
  double dec_deg()const { return _dec * deg;}

  /*! return a track-diretion in UTM/detector coordinates */
  Vec track_direction( const Det& det,
                       const TTimeStamp& t,
                       bool j2000 = true );

  /*! return declination in the form deg, min, sec */
  const char* dec_str() const;

  /*! return right assension in the form hours, min, sec */
  const char* ra_str() const;

  const char* __str__() const;

  /* Undo nutation and precession since J2000 so that ra and dec refer to that epoch.
     The reverse flag is used from going from J2000 to current(t) epoch.  */
  EquatorialCoords correct_to_j2000( const TTimeStamp& t ,
                                     bool reverse = false ) const;

  /*! return the angle between two points in the sky */

  double distance( const EquatorialCoords& other ) const
  {
    return angle_between( asVec() , other.asVec());
  }

  /* return a 3d-vector representing the spherical coordinates (for e.g. computing the angle between two objects) */
  Vec asVec() const
  {
    Vec v; v.set_angles( pi / 2 - dec() , ra() );
    return v;
  }

  /*! complement of asVec */
  void fromVec(const Vec& v) 
  {
    set_dec ( pi/2 - v.theta() );
    set_ra  ( v.phi() );
  }

  ClassDefNV(EquatorialCoords, 1)
};

// as of ROOT 6.16?, pyroot seems to use printValue to print
namespace cling
{
inline std::string printValue(const EquatorialCoords* t)
{
  return t->__str__();
}
}


/*! Coordinates for use in slalib. In this system, azimuth=0 points to
  (geodetic) North, and azimuth = pi/2 is east. */

struct SlaCoords
{
  double _azi, _ele;

  SlaCoords() {}
  SlaCoords(double azimuth, double elevation) : _azi(azimuth), _ele(elevation) {}

  double azimuth()       const { return _azi;     }
  double elevation()     const { return _ele;     }
  double azimuth_deg()   const { return _azi * deg; }
  double elevation_deg() const { return _ele * deg; }
};


inline Vec sla_to_utm( double azimuth,
                       double elevation,
                       double meridian_convergence_angle )
{
  double az  = wrap( pi / 2 - azimuth + meridian_convergence_angle );
  double zen = pi / 2 - elevation;
  Vec r; r.set_angles(zen, az);
  return r;
}

inline SlaCoords utm_to_sla( double azimuth,
                             double zenith,
                             double meridian_convergence_angle )
{
  // for slalib, the sign convention for azimuth is north zero,
  // east +π/2.

  double az = wrap( pi / 2 - azimuth + meridian_convergence_angle );
  double elevation = pi / 2. - zenith;
  return {az, elevation};
}

#include "TMultiGraph.h"
#include "TGraph.h"
#include "TPolyLine.h"
#include "TLatex.h"
#include "TMarker.h"
#include "rootutil.hh"



/*! Facility to plot skymaps in ROOT. */

struct SkyMap : public TObject
{
  int ngridlines_x;
  int ngridlines_y;
  string mode;
  bool input_degrees;
  TMultiGraph* grid;
  vector<TObject*> objects;
  TLatex labelfont;
  EquatorialCoords center;


  SkyMap( string mode = "eq", bool input_degrees = false , EquatorialCoords center = {0,0}) :
    mode(mode), input_degrees( input_degrees ), grid(nullptr), center(center)
  {
    labelfont.SetTextSize( 0.03 );
  }



  /*! add a grid with nx x ny lines */
  void SetGrid(int nx = 9, int ny = 11)
  {
    if (grid) delete grid;
    grid = sky_grid(nx, ny);
  }

  void zoom( EquatorialCoords coords, 
             double scalex, double scaley = -1 ) 
  {
    if (!grid) SetGrid();
    if ( scaley == -1 ) scaley = scalex;
    double xx,yy;
    transform( coords , xx, yy );
    grid->GetXaxis()->SetRangeUser( xx - scalex, xx + scalex );
    grid->GetYaxis()->SetRangeUser( yy - scaley, yy + scaley );

  }


  /*! Draw the skymap */
  void Draw()
  {
    if (!grid) SetGrid();
    grid->Draw("AL");
    grid->GetXaxis()->SetNdivisions(0);
    grid->GetYaxis()->SetNdivisions(0);

    if (mode == "eq" )
      {
        draw_scale_ra_hours();
        draw_scale_dec_degs();
      }

    for ( auto& o : objects ) o->Draw("same");
  }


  void transform( const EquatorialCoords& coords , double& xx, double& yy )
  {
    transform( coords.ra(), coords.dec(), xx, yy );
  }

  void transform(double x, double y, double& xx, double& yy )
  {
    if ( input_degrees )
    {
      x /= deg;
      y /= deg;
    }

    vector<double> m;
    if (mode == "raw")
    {
      m = hammer_projection(x, y);
    }

    if (mode == "eq" ) // 0 on the right 24h on the left
    {
      m = hammer_projection(pi - x, y);
    }

    xx = m[0]; yy = m[1];
    return;

  }


  // void add_circle( double x, double y, double radius);
  // void add_line( double x1, doubel y1, double x1, double x2)

  template< typename T>
  void transform_and_add( const T& m )
  {
    T&  mm = *(new T(m));
    double x, y;
    transform( m.GetX(), m.GetY(), x, y );
    mm.SetX(x); mm.SetY(y);
    objects.push_back( (TObject*) &mm );
  }

  template< typename T>
  void transform_and_add( const T& m , EquatorialCoords coords )
  {
    T&  mm = *(new T(m));
    double x, y;
    transform( coords.ra() , coords.dec() , x, y );
    mm.SetX(x); mm.SetY(y);
    objects.push_back( &mm );
  }

  /*! add a TMarker at the coordinates */
  void add (const TMarker& m , EquatorialCoords coords )
  {
    transform_and_add( m, coords );
  }

  /*! add a TMarker, using its position as coordinates. */
  void add( const TMarker& m )
  {
    transform_and_add( m );
  }

  /*! Add TLatex text. */
  void add( TLatex& l )
  {
    transform_and_add( l );
  }

  /*! Add TLatex text, at specified coordinates */
  void add( TLatex& l,  EquatorialCoords coords )
  {
    transform_and_add( l, coords );
  }

  /*! Add some text, at specified coordinates. */
  void add( const char* latex,
            EquatorialCoords coords )
  {
    TLatex m(0, 0, latex);
    add( m, coords );
  }

  /*! Add a circle at the specified coordinates. User can set the style
      of returned tgraph object and is responsible to delete it. */

  TGraph* add_circle( EquatorialCoords coords,
                      double radius, unsigned npoints = 100 )
  {
    Vec v;
    Mat R; R.set( coords.asVec() );

    vector<double> xx,yy;
    xx.resize( npoints + 1 );
    yy.resize( npoints + 1 );

    for (unsigned i=0; i< npoints+1; i++ ) {
      double phi = i / double(npoints) * 2* pi  ;
      v.set_angles( radius, phi );
      Vec p = R * v;
      transform ( p.phi(), pi/2 - p.theta() , xx[i], yy[i] );
    }

    auto g = new TGraph ( xx.size() , &xx[0], &yy[0] );
    objects.push_back( (TObject*) g );
    return g;
  }



  /*! Add scale indication in hours. */
  void draw_scale_ra_hours(int interval = 6 )
  {
    EquatorialCoords eq;

    for (int h = 0; h <= 24; h += interval )
    {
      eq.set_ra ( h, 0, 0 );
      eq.set_dec( 0.0     );

      TLatex m = *( labelfont.DrawLatex( 0, 0, Form("%dh", h) ));
      add( m, eq );
    }
  }

  /*! Add scale indication in hours. */
  void draw_scale_dec_degs( int interval = 15 ,
                            bool skip_zero = true,
                            double xcoord = 2 * pi )
  {
    for (int d = -90 + interval; d <= 90 - interval; d += interval )
    {
      if ( skip_zero && d == 0 ) continue;

      EquatorialCoords eq( xcoord, d / deg );
      TLatex m(0, 0, Form("%d^{o}", d) );
      if ( d < 0 ) m.SetTextAlign(33);
      else         m.SetTextAlign(31);
      m.SetTextSize( 0.03 );
      add( m, eq );
    }
  }


  ClassDef(SkyMap, 1);

};

#endif