//____________________________________________________________________________
/*!

\class   GAstro

\brief   A class for astronomical calculation
         Astronomical algorithms implemented in this library are based on
          J. Meeus, Astronomical Algorithms, Willmann-Bell, Incorporated, 1991.
         Equatorial coordinates refer to 2000.0 epoch

\author  Carla Distefano <distefano_c@lns.infn.it>
         LNS-INFN, Catania

\created February 1, 2016

\cpright  Copyright (c) 2016-2019, The KM3NeT Collaboration
          For the full text of the license see $GSEAGEN/LICENSE
*/
//____________________________________________________________________________


#include <fstream>
#include <cmath>
#include <iostream>
#include <fstream>
#include <string>
#include <cstdio>
#include <sstream>
#include <stdlib.h>
#include <cmath>
#include <cstring>


#ifdef _GENIEVERSIONGTEQ3__
#include "Framework/Conventions/Constants.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Utils/StringUtils.h"
#else
#include "Conventions/Constants.h"
#include "Messenger/Messenger.h"
#include "Utils/StringUtils.h"
#endif

#include "SeaNuDrivers/GAstro.h"

using namespace std;
using namespace genie;
using namespace genie::constants;


const double  Deg2Rad = kPi/180.;
const double  Rad2Hour= 12./kPi;


GAstro::GAstro(double Latitude, double Longitude){
// The longitude is positive eastwards from the meridian of Greenwich, and negative to the West.

  fLatitude = Latitude;
  fLongitude = Longitude;
  fUTMConvergeAngle = 0.;
  fUseUTMSystem = false;

  LOG("GAstro", pINFO)<< "Initialization of GAstro with Latitude = "<<fLatitude<<", Longitude = "<<fLongitude<<" rad";
}
//___________________________________________________________________________
GAstro::GAstro(double Latitude, double Longitude, double UTMConvergeAngle, bool UseUTMSystem){
// The longitude is positive eastwards from the meridian of Greenwich, and negative to the west.

  fLatitude = Latitude;
  fLongitude = Longitude;
  fUTMConvergeAngle = UTMConvergeAngle;
  fUseUTMSystem = UseUTMSystem;

  LOG("GAstro", pINFO)<< "Initialization of GAstro with Latitude = "<<fLatitude<<", Longitude = "<<fLongitude<<" rad";
}
//___________________________________________________________________________
GAstro::~GAstro(){  }
//___________________________________________________________________________
void GAstro::SetCoordinates(double Latitude, double Longitude){
// The longitude is positive westwards from the meridian of Greenwich, and negative to the east.

  fLatitude = Latitude;
  fLongitude = Longitude;

  return;
}
//___________________________________________________________________________
double GAstro::GST(double MJD){  // Greenwich Sidereal Time as a function of the Modified Julian Day
  double gst,T;
  T=(MJD-51544.5)/36525.0;
  gst=(280.46061837+13185000.77*T+(pow(T,2)*(0.000387933-T/38710000.0)))/15.;
  gst=fmod(gst,24.0); //normalization between 0 e 24 hours
  return gst/Rad2Hour; // in radians
}
//___________________________________________________________________________
double GAstro::LST(double MJD){ // Local Sidereal Time as a function of the Modified Julina Day and the observer's longitude  // dire come si deve calcolare la longitudine
  double lst,gst;
  gst=this->GST(MJD);
  lst=gst+fLongitude;
  lst=fmod(lst,2*kPi);          //normalization between 0 e 24 hours
  return lst;                   // in radians
}
//___________________________________________________________________________
void GAstro::Eq2Hor(double lst, double RA, double Dec, double * Azimuth, double * altitude){ // transforms equatorial coordinates into horizontal ones
  double HA=lst-RA; //hour angle in radians
  *altitude =  asin(sin(fLatitude)*sin(Dec)+cos(fLatitude)*cos(Dec)*cos(HA));
  *Azimuth=atan2(sin(HA),(cos(HA)*sin(fLatitude)-tan(Dec)*cos(fLatitude)))+kPi; // from North
  *Azimuth=fmod(*Azimuth+2*kPi,2*kPi);
  return;
}
//___________________________________________________________________________
void GAstro::Hor2Eq(double Altitude, double Azimuth, double lst, double * RA, double * Dec){ // transforms horizontal into equatorial ones
  Azimuth-=kPi;  // azimth is from north, now from south
  *Dec=asin(sin(fLatitude)*sin(Altitude)-cos(fLatitude)*cos(Altitude)*cos(Azimuth));
  double HA=atan2(sin(Azimuth),cos(Azimuth)*sin(fLatitude)+tan(Altitude)*cos(fLatitude));
  *RA=lst-HA;
  *RA=fmod(*RA+2*kPi,2*kPi);
  return;
}

//___________________________________________________________________________
void GAstro::Gal2Eq(string Equinox, double GalLat, double GalLong, double * RA, double * Dec){// transforms galactic coordinates into equatorial ones

  double PoleRA=0.,PoleDec=0.,PosAngle=0.;

  if(Equinox.compare("J2000")==0){
    PoleRA = 192.859508;
    PoleDec = 27.128336;
    PosAngle = 122.932;
  }
  else if(Equinox.compare("B1950")==0){
    PoleRA = 192.25;
    PoleDec = 27.4;
    PosAngle = 123.;
  }
  else {
    LOG("GAstro", pFATAL)<< "Unknown Equinox, it must be J2000 or B1950, exiting...";
    gAbortingInErr = true;
    exit(1);
  }


  PoleRA *= Deg2Rad;
  PoleDec *= Deg2Rad;
  PosAngle *= Deg2Rad;


  *Dec = asin(sin(PoleDec)*sin(GalLat)+cos(PoleDec)*cos(GalLat)*cos(PosAngle-GalLong));
  *RA = atan2(cos(GalLat)*sin(PosAngle-GalLong),cos(PoleDec)*sin(GalLat)-sin(PoleDec)*cos(GalLat)*cos(PosAngle-GalLong))+PoleRA;
  *RA=fmod(*RA+2*kPi,2*kPi);


  return;
}
//___________________________________________________________________________
void GAstro::Eq2Gal(string Equinox, double RA, double Dec, double * GalLat, double * GalLong){// transforms equatorial coordinates into galactic ones

  double PoleRA=0.,PoleDec=0.,PosAngle=0.;

  if(Equinox.compare("J2000")==0){
    PoleRA = 192.859508;
    PoleDec = 27.128336;
    PosAngle = 122.932;
  }
  else if(Equinox.compare("B1950")==0){
    PoleRA = 192.25;
    PoleDec = 27.4;
    PosAngle = 123.;
  }
  else {
    LOG("GAstro", pFATAL)<< "Unknown Equinox, it must be J2000 or B1950, exiting...";
    gAbortingInErr = true;
    exit(1);
  }

  PoleRA *= Deg2Rad;
  PoleDec *= Deg2Rad;
  PosAngle *= Deg2Rad;

  *GalLat = asin(sin(PoleDec)*sin(Dec)+cos(PoleDec)*cos(Dec)*cos(RA-PoleRA));
  *GalLong = PosAngle-atan2(cos(Dec)*sin(RA-PoleRA),cos(PoleDec)*sin(Dec)-sin(PoleDec)*cos(Dec)*cos(RA-PoleRA));
  *GalLong=fmod(*GalLong+2*kPi,2*kPi);

  return;
}

//___________________________________________________________________________
void GAstro::CalcPrecession(double MJD0, double MJD, double RA0, double Dec0, double *RA, double *Dec){// calc correction to equtorial cooridnates due to the precession

  double JD0=MJD0+2400000.5;
  double JD=MJD+2400000.5;

  double T= (JD0-2451545.0)/36525.;
  double t= (JD-JD0)/36525.;

  double zeta = (2306.2181 + 1.39656*T - 0.000139*T*T)*t;
         zeta+= (0.30188 - 0.000344*T)*t*t;
         zeta+= 0.017998*t*t*t;

  double z = (2306.2181 + 1.39656*T - 0.000139*T*T)*t;
         z+= (1.09468 + 0.000066*T)*t*t;
         z+= 0.018203*t*t*t;

  double theta = (2004.3109 - 0.85330*T - 0.000217*T*T)*t;
         theta += -(0.42665 + 0.000217*T)*t*t;
         theta += -0.041833*t*t*t;


  zeta/=3600.;
  z/=3600;
  theta/=3600.;

  zeta *= Deg2Rad;
  z *= Deg2Rad;
  theta *= Deg2Rad;

  double A = cos(Dec0)*sin(RA0+zeta);
  double B = cos(theta)*cos(Dec0)*cos(RA0+zeta)-sin(theta)*sin(Dec0);
  double C = sin(theta)*cos(Dec0)*cos(RA0+zeta)+cos(theta)*sin(Dec0);

  *Dec=asin(C);
  *RA = atan2(A,B)+z;
  *RA=fmod(*RA+2*kPi,2*kPi);


  return;
}
//___________________________________________________________________________
void GAstro::CalcNutation(double MJD, double RA0, double Dec0, double *RA, double *Dec){

  double JD=MJD+2400000.5;
  double T= (JD-2451545.0)/36525.;
  
  double L = 280.46665+36000.7698*T+0.000303*T*T;
  double L1 = 218.3165+481267.8813*T- 0.001599*T*T;
  double M = 357.52772 + 35999.050340*T - 0.000154*T*T;
  double M1 = 134.96298 + 477198.867398*T + 0.008721*T*T;
  double Omega = 125.0443 - 1934.1363* T + 0.002075*T*T;
  
  Omega*=Deg2Rad;
  L*=Deg2Rad;
  L1*=Deg2Rad;
  M*=Deg2Rad;
  M1*=Deg2Rad;

  double DeltaPsi = - (17.1996 + 0.01742 * T) * sin(Omega);
	 DeltaPsi += + (1.3187 + 0.00016 * T) * sin (2*L);
	 DeltaPsi += - 0.2274 * sin (2*L1);
	 DeltaPsi += + 0.2062 * sin (2*Omega);
	 DeltaPsi += + (0.1426 - 0.00034 * T) * sin(M);
	 DeltaPsi += + 0.0712 * sin (M1);
	 DeltaPsi += - (0.0517 - 0.00012 * T) * sin (2*L + M);
	 DeltaPsi += - 0.0386 * sin (2*L1 - Omega);
	 DeltaPsi += - 0.0301 * sin (2*L1 + M1);
	 DeltaPsi += + 0.0217 * sin (2*L - M);
	 DeltaPsi += - 0.0158 * sin (2*L - 2*L1 + M1);
	 DeltaPsi += + 0.0129 * sin (2*L - Omega);
	 DeltaPsi += + 0.0123 * sin (2*L1 - M1);

  double DeltaEpsilon = + (9.2025 + 0.00089 * T) * cos(Omega);
         DeltaEpsilon += + (0.5736 - 0.00031 * T) * cos (2*L);
         DeltaEpsilon += + 0.0977 * cos (2*L1);
         DeltaEpsilon += - 0.0895 * cos (2*Omega);
         DeltaEpsilon += + 0.0224 * cos (2*L + M);
         DeltaEpsilon += + 0.0200 * cos (2*L1 - Omega);
         DeltaEpsilon += + 0.0129 * cos (2*L1 + M);
         DeltaEpsilon += - 0.0095 * cos (2*L - M);
         DeltaEpsilon += - 0.0070 * cos (2*L - Omega);


  DeltaPsi /= 3600.;		// in deg
  DeltaEpsilon /= 3600.;	// in deg
  
  double Epsilon0 = -46.8150*T -0.00059*T*T + 0.001813*T*T*T;
         Epsilon0 /=3600.;	// in deg
	 Epsilon0 += 23. + 26./60. + 21.448/3600.;
  
  double Epsilon = Epsilon0 + DeltaEpsilon;
  
  Epsilon *= Deg2Rad;
  DeltaPsi *= Deg2Rad;
  DeltaEpsilon *= Deg2Rad;
  
  
  double DeltaRA = (cos(Epsilon)+sin(Epsilon)*sin(RA0)*tan(Dec0))*DeltaPsi - (cos(RA0)*tan(Dec0))*DeltaEpsilon;
  double DeltaDec = (sin(Epsilon)*cos(RA0))*DeltaPsi+sin(RA0)*DeltaEpsilon;
  
  *RA=RA0+DeltaRA;
  *Dec=Dec0+DeltaDec;
  
  return;
}
//___________________________________________________________________________
void GAstro::Hor2Direction(double Azimuth, double altitude, double * vx, double * vy, double * vz){ // computes direction cosines corresponding to the source horizontal coordinates

  // source theta and phi angles in radians
  double AstroTheta= kPi/2.-altitude;
  double AstroPhi=2*kPi-Azimuth;

  if(fUseUTMSystem) AstroPhi += fUTMConvergeAngle + kPi/2.;


  // source direction cosines
  *vx=sin(AstroTheta)*cos(AstroPhi);
  *vy=sin(AstroTheta)*sin(AstroPhi);
  *vz=cos(AstroTheta);

  return;
}
//___________________________________________________________________________
void GAstro::Direction2Hor(double vx, double vy, double vz,double * Azimuth, double * Altitude){ // computes source horzontal coordinates from the corresponging direction cosines

  // source theta and phi angles in radians
  double AstroTheta= acos(vz);
  double AstroPhi=atan2(vy,vx);

  if(fUseUTMSystem) AstroPhi -= fUTMConvergeAngle + kPi/2.;


  *Altitude=kPi/2.-AstroTheta;
  *Azimuth=2*kPi-AstroPhi;
  *Azimuth=fmod(*Azimuth+2*kPi,2*kPi);

  return;
}
//___________________________________________________________________________
double GAstro::CalcMJD(int M, int day, int Dh, int Dm, int Ds, int Y){ //computes the Modified Julian Day

  double D;
  int A,B;
  double MJD;

  D=(double)day+((double)Dh+(double)Dm/60.+(double)Ds/3600.)/24.;

  if(M == 1 || M  == 2){
    Y=Y-1;
    M=M+12;
  }

  A = (int)(Y/100);
  B = 2 - A +(int)(A / 4);
  MJD = (int)(365.25*Y) + (int)(30.6001*(M+1)) + D + 1720994.5 + B;  // Julian day
  MJD-=2400000.5;

  return MJD;
}
//___________________________________________________________________________
double GAstro::CalcMJD(string Date){ //computes the Modified Julian Day
  // Date in format dd-mm-yyyy;hh:mm:ss

  int M, day, Dh, Dm, Ds, Y;
  bool FormatOk = true;

  if(Date.size()!=19)FormatOk=false;

  if(Date.substr(2,1)!="-")FormatOk=false;
  if(Date.substr(5,1)!="-")FormatOk=false;
  if(Date.substr(10,1)!=";")FormatOk=false;
  if(Date.substr(13,1)!=":")FormatOk=false;
  if(Date.substr(16,1)!=":")FormatOk=false;

  if(!FormatOk){
    LOG("GAstro", pFATAL)<< "Wrong format for Date "<<Date<<". It must be dd-mm-yyyy;hh:mm:ss, exiting...";
    gAbortingInErr = true;
    exit(1);
  }

  string::size_type semicolon = Date.find_first_of(";",0);
  string Date1 = Date.substr(0,semicolon);
  Date.erase(0,semicolon+1);

  vector<string> vDate = utils::str::Split(Date1, "-");
  assert(vDate.size() == 3);

  day  = atoi(vDate[0].c_str());
  M    = atoi(vDate[1].c_str());
  Y    = atoi(vDate[2].c_str());

  vDate = utils::str::Split(Date, ":");
  assert(vDate.size() == 3);

  Dh  = atoi(vDate[0].c_str());
  Dm  = atoi(vDate[1].c_str());
  Ds  = atoi(vDate[2].c_str());

  double MJD=CalcMJD(M, day, Dh, Dm, Ds, Y);

  return MJD;
}
//___________________________________________________________________________
double GAstro::CalcDate(double MJD, string * Date){

  double JD;
  double Z,F,A;
  double B,C,D,E;

  JD=MJD+2400000.5;     // Julian date
  JD+=0.5;

  F = modf (JD , &Z);

  if(Z<2299161)A=Z;
  else{
    double alpha,alpha1;
    modf( ((Z - 1867216.25 ) / 36524.25)  ,&alpha);
    modf(alpha/4.,&alpha1);
    A = Z + 1 + alpha - alpha1;
  }

  B=A+1524;

  modf((( B - 122.1 ) / 365.25), &C);

  modf(365.25 * C, &D);

  modf( (( B - D ) / 30.6001) ,&E);

  double PartInt,PartFrac;
  modf(30.6001 * E, &PartInt);

  double dayall= B - D + F - PartInt;

  PartFrac = modf(dayall,&PartInt);

  int day = (int)PartInt;

  int month;

  if((int)E<14)month=(int)E-1;
  else month=(int)E-13;

  int year=0;

  if(month>2)year=(int)C-4716;
  else year=(int)C-4715;


  double hh,mm,ss;

  PartFrac*=24.;
  PartFrac=modf(PartFrac,&hh);
  int hour=(int)hh;

  PartFrac*=60.;
  PartFrac=modf(PartFrac,&mm);
  int min=(int)mm;

  PartFrac*=60.;
  PartFrac=modf(PartFrac,&ss);
  int sec=(int)ss;

  char ch[20];
  sprintf(ch,"%2.2d-%2.2d-%4.4d;%2.2d:%2.2d:%2.2d",day,month,year,hour,min,sec);

  Date->clear();
  Date->append(ch);

  return PartFrac;

}
//___________________________________________________________________________
void GAstro::CalcTimeStamp(double MJD, unsigned int * TimeStampSec, unsigned int * TimeStampTick){
// computes the TimeStamp as number of seconds counted since 01.01.1970 00:00:00  UTC and 16 nanosecond-ticks.

  double MJD0=40587.;           // MJD corresponding to 01.01.1970 00:00:00 UTC

  double PartInt,PartFrac;

  PartFrac=modf(((MJD-MJD0)*24.*3600.),&PartInt);
  *TimeStampSec=(unsigned int)PartInt;

  PartFrac*=62500000;
  *TimeStampTick=(unsigned int)PartFrac;

  return;
}
//___________________________________________________________________________
double GAstro::CalcMJDFromTimeStamp(unsigned int TimeStampSec, unsigned int TimeStampTick){
// computes MJD from the TimeStamp as number of seconds counted since 01.01.1970 00:00:00  UTC and 16 nanosecond-ticks.

  double MJD0=40587.;           // MJD corresponding to 01.01.1970 00:00:00 UTC

  double MJD=(double)TimeStampSec+(double)TimeStampTick/62500000;

  MJD/=24.*3600.;
  MJD+=MJD0;

  return MJD;
}
//___________________________________________________________________________
void GAstro::CalcSrcDirMJD(double MJD, double RA0, double Dec0, double * vx, double * vy, double * vz){
// computes the direction cosines corresponding to the source as a function of the MJD, observer's coordinates and source equatorial coordinates

  double lst;
  double RA,Dec; // after correction due to precession (RA0 and Dec0 referring to 2000)
  double MJD0=51544.5;  // J2000.0
  lst=this->LST(MJD);
  this->CalcPrecession(MJD0, MJD, RA0, Dec0, &RA, &Dec);
  this->CalcNutation(MJD,RA,Dec,&RA,&Dec);
  
  this->CalcSrcDirLST(lst, RA, Dec, vx, vy, vz);
  
  
  return;
}
//___________________________________________________________________________
void GAstro::CalcSrcDirMJD(double MJD, double RA, double Dec, double * CosTheta, double * Phi){
// computes the direction angles corresponding to the source as a function of the MJD, observer's coordinates and source equatorial coordinates

  double vx,vy,vz;
  this->CalcSrcDirMJD(MJD, RA, Dec, &vx, &vy, &vz);
  *CosTheta=vz;
  *Phi=atan2(vy,vx);

  return;
}
//___________________________________________________________________________
void GAstro::CalcSrcDirLST(double lst, double RA, double Dec, double * vx, double * vy, double * vz){
// computes the direction cosines corresponding to the source as a function of the MJD, observer's coordinates and source equatorial coordinates

  double Azimuth, altitude;
  this->Eq2Hor(lst, RA, Dec, &Azimuth,&altitude);
  this->Hor2Direction(Azimuth,altitude, vx,vy,vz);
//  cout<<"source: "<<LST<<" "<<Azimuth<<" "<<altitude<<" "<<-*vx<<" "<<-*vy<<" "<<-*vz<<endl;


  return;
}
//___________________________________________________________________________
void GAstro::CalcSrcDirLST(double lst, double RA, double Dec, double * CosTheta, double * Phi){
// computes the direction angles corresponding to the source as a function of the MJD, observer's coordinates and source equatorial coordinates

  double vx,vy,vz;
  this->CalcSrcDirLST(lst, RA, Dec, &vx, &vy, &vz);
  *CosTheta=vz;
  *Phi=atan2(vy,vx);

  return;
}