#include "Astro.hh"
#include "Evt.hh"
#include "Trk.hh"
#include "Det.hh"
#include "TTimeStamp.h"

#include "stringutil.hh"
#include "sla.hh"


/* return declination in the form deg, min, sec */
const char* EquatorialCoords::dec_str() const
{
  int sign = dec()<0? -1 : 1;

  double d = sign * dec_deg();

  int degrees    = int(d);
  int arcmins    = int( (d-degrees) * 60 );
  double arcsecs = (d-degrees - arcmins/60.0 ) * 3600;
		   
  return Form("%d°%d\'%4.2f\"", sign*degrees, arcmins, arcsecs );
}

/* return right assension in the form hours, min, sec */
const char* EquatorialCoords::ra_str() const
{
  double h = ra_deg() / 15.0;  
  int hour = int( h );
  int mins = abs(int( (h-hour)*60 ));
  double secs = abs((h-hour-mins/60.0)*3600);
  
  return Form("%dh %dm %4.2fs", hour, mins, secs );
}

const char* EquatorialCoords::__str__() const
{
  return Form("%s, %s", ra_str(), dec_str() );
}

EquatorialCoords EquatorialCoords::correct_to_j2000( const TTimeStamp& t ,
						     bool reverse /*=false*/ ) const
{
  // --- get the rotation matrix ----
  // NOTE: (from seatray Astro): Julian epoch of J2000 time definition
  // is 2000, corrsponding by definition to 2451545.0 TT (unmodified)
  // Julian date!
  
  double epoch = 2000;
  double mjd   = t.AsJulianDate() - 2400000.5;

  auto m = sla::prenut ( epoch, mjd );

  Mat MT; MT.set(m.RMATPN);
  Mat M = MT.T();  // transpose 'cause of fortan matrix convention in RMATPN
 
  
  // if reverse = false, this function converts from ra,dec at t to J2000,
  // if reverse = true, it goes the other way (i.e. if you have a catalog
  // position and want to compute a track, use reverse = true)

  Vec v = asVec();
  Vec v2 = reverse? M * v : M.inv() * v;
  return { wrap(v2.phi()), pi/2-v2.theta()};
}


EquatorialCoords::EquatorialCoords( const Det& det,
				    const Vec& track_direction,
				    const TTimeStamp& t,
				    bool j2000 /*= true*/ )
{
  Vec d = -track_direction; // track points back to its origin
  
  SlaCoords s = utm_to_sla( d.phi(), d.theta(), det.meridian_convergence_angle );
  
  auto equatorial = sla::dh2e( s.azimuth(),
			       s.elevation(),
			       det.latitude );
			       
  // rightascension = utm sidereal time - hour angle
  
  // TTimeStamp::AsLMST could be corrected for the difference between UTC and UT1
  // (leap seconds), which would improve the accuracy by <1 s.
  double lstRad = t.AsLMST( det.longitude * deg ) / 24 * 2 * pi;
  double raRad = lstRad - equatorial.HA;
    
  set( raRad, equatorial.DEC );

  if ( j2000 ) *this = correct_to_j2000( t );
}

Vec EquatorialCoords::track_direction( const Det& det,
				       const TTimeStamp& t,
				       bool j2000 /*= true*/ )
{
  auto nowcoords = j2000 ? correct_to_j2000( t, true ): *this; // nb: reverse = true

  
  double lst = t.AsLMST( det.longitude * deg ) / 24 * 2 * pi; // in radians

  double hour_angle =  lst - nowcoords.ra();
  auto hor = sla::de2h( hour_angle , nowcoords.dec(), det.latitude );
  Vec r = sla_to_utm( hor.AZ, hor.EL, det.meridian_convergence_angle );
  
  // minus because it points back to the source
  return -r;
}