#include "sun.h"

class Sun::JdTime
{
public:
	static constexpr double JULIAN_OFFSET = 68.184 / 86400.0;
	static constexpr double JULIAN_DATEATY2K = 2451544.500000;

	double julianDate;
	double daysSinceY2K;
	double siderialUtcTime;
	double localHour;
	double utcHour;
	int    day;
	int    year;

	JdTime();
	void update();
	std::time_t toStdTime();
	void fromStdTime(std::time_t time);
};

Sun::JdTime::JdTime()
{
	update();
}

void Sun::JdTime::update()
{
	std::time_t time = std::time(nullptr);
	fromStdTime(time);
}


void Sun::JdTime::fromStdTime(std::time_t time)
{
	std::tm     ltime = *std::localtime( &time );

	localHour = ltime.tm_hour + ltime.tm_min/60.0 + ltime.tm_sec/3600.0;
	utcHour = localHour - ltime.tm_gmtoff/3600.0;
	day = ltime.tm_yday;
	year = ltime.tm_year + 1900;

	std::tm millennium = {0,0,0,1,0,100};
	std::time_t millenniumTime = std::mktime(&millennium) - timezone;

	daysSinceY2K = (time - millenniumTime)/(60*60*24.0);
	//std::cout<<"days since y2k: "<<daysSinceY2K<<std::endl;
	julianDate = daysSinceY2K + JULIAN_DATEATY2K - JULIAN_OFFSET;
	//std::cout<<"julianDate: "<<julianDate<<std::endl;

	siderialUtcTime = 6.697374558 + 0.06570982441908 * (long)daysSinceY2K + 1.00273790935*utcHour;
	while(siderialUtcTime < 0)  siderialUtcTime+=24;
	siderialUtcTime = fmod(siderialUtcTime, 24);
}

std::time_t Sun::JdTime::toStdTime()
{
	std::tm millennium = {0,0,0,0,0,100};
	std::time_t millenniumTime = std::mktime(&millennium)- timezone;

	millenniumTime = millenniumTime + (julianDate-JULIAN_DATEATY2K+JULIAN_OFFSET)*(24*60*60);
	return millenniumTime;
}

Sun::Sun(double latitude, double longitude, double altitude): latitude_(latitude), longetude_(longitude),
	altitude_(altitude)
{}

double Sun::nextMeanSolarNoonJD(const JdTime& time)
{
	return (long)time.daysSinceY2K + 2451545.0 + 0.0009 - (longetude_ / 360);
}

double Sun::meanSolarAnomaly(double meanSolarNoon)
{
	return fmod(357.5291 + 0.98560028 *  (meanSolarNoon-2451545), 360);
}

double Sun::eqOfCenter(double meanSolarAnomaly)
{
	return 1.9148*sin(meanSolarAnomaly*TO_RADS) + 0.0200*sin(2*meanSolarAnomaly*TO_RADS) + 0.0003*sin(
	           3*meanSolarAnomaly*TO_RADS);
}

double Sun::eclipticLongitude(double eqOfCenter, double meanSolarAnomaly)
{
	return fmod(meanSolarAnomaly + eqOfCenter + ARGUMENT_OF_PERIHELION + 180, 360); //hmmm
}

double Sun::equationOfTime(double meanSolarAnomaly, double eclipticLongitude) //wrong
{
	return (-7.659*sin(meanSolarAnomaly*TO_RADS)+9.863*sin(2*(meanSolarAnomaly+eclipticLongitude)*TO_RADS))/60;
}

double Sun::localSolarTime(const JdTime& time, double equationOfTime)
{
	return time.utcHour + equationOfTime + (4/60.0)*longetude_;
}

double Sun::solarDeclination(double eclipticLongitude)
{
	return TO_DEGS*asin(sin(eclipticLongitude*TO_RADS)*sin(PLANETARY_AXIAL_TILT*TO_RADS));
}

double Sun::hourAngle(double localSolarTime)
{
	return 360/24 * (localSolarTime - 12);
}

double Sun::hourAngleAtSunset(double solarDeclination)
{
	return TO_DEGS*acos((sin((-0.83-(2.076*sqrt(altitude_)/60.0))*TO_RADS) - sin(solarDeclination*TO_RADS)*sin(
	                         latitude_*TO_RADS))/(cos(latitude_*TO_RADS)*cos(solarDeclination*TO_RADS)));
}

double Sun::altitude()
{
	JdTime time;
	double meanSolarNoonValue     = nextMeanSolarNoonJD(time);
	double meanSolarAnomalyValue  = meanSolarAnomaly(meanSolarNoonValue);
	double eclipticLongitudeValue = eclipticLongitude(eqOfCenter(meanSolarAnomalyValue), meanSolarAnomalyValue);
	double localSolarTimeValue    = localSolarTime(time, equationOfTime(meanSolarAnomalyValue, eclipticLongitudeValue));
	double declinationValue       = solarDeclination(eclipticLongitudeValue);

	double cosZenithAngle = sin(latitude_*TO_RADS)*sin(declinationValue*TO_RADS)+cos(latitude_*TO_RADS)*cos(
	                            declinationValue*TO_RADS)*cos(hourAngle(localSolarTimeValue)*TO_RADS);

	return TO_DEGS*asin(cosZenithAngle);
}

double Sun::maximumAltitude()
{
	JdTime time;
	double meanSolarNoonValue     = nextMeanSolarNoonJD(time);
	double meanSolarAnomalyValue  = meanSolarAnomaly(meanSolarNoonValue);
	double eclipticLongitudeValue = eclipticLongitude(eqOfCenter(meanSolarAnomalyValue), meanSolarAnomalyValue);
	double declinationValue       = solarDeclination(eclipticLongitudeValue);

	double cosZenithAngle = sin(latitude_*TO_RADS)*sin(declinationValue*TO_RADS)+cos(latitude_*TO_RADS)*cos(
	                            declinationValue*TO_RADS);

	return TO_DEGS*asin(cosZenithAngle);
}

double Sun::azimuth()
{
	return 0;
}


std::time_t Sun::riseTime()
{
	JdTime time;
	double meanSolarNoonValue     = nextMeanSolarNoonJD(time);
	double meanSolarAnomalyValue  = meanSolarAnomaly(meanSolarNoonValue);
	double eclipticLongitudeValue = eclipticLongitude(eqOfCenter(meanSolarAnomalyValue), meanSolarAnomalyValue);
	double declinationValue       = solarDeclination(eclipticLongitudeValue);
	double hourAngleValue         = hourAngleAtSunset(declinationValue);

	time.julianDate = meanSolarNoonValue + equationOfTime(meanSolarAnomalyValue,
	                  eclipticLongitudeValue) - hourAngleValue/360.0;

	return time.toStdTime();
}
std::time_t Sun::setTime()
{
	JdTime time;
	double meanSolarNoonValue     = nextMeanSolarNoonJD(time);
	double meanSolarAnomalyValue  = meanSolarAnomaly(meanSolarNoonValue);
	double eclipticLongitudeValue = eclipticLongitude(eqOfCenter(meanSolarAnomalyValue), meanSolarAnomalyValue);
	double declinationValue       = solarDeclination(eclipticLongitudeValue);
	double hourAngleValue         = hourAngleAtSunset(declinationValue);

	time.julianDate = meanSolarNoonValue + equationOfTime(meanSolarAnomalyValue,
	                  eclipticLongitudeValue) + hourAngleValue/360.0;

	return time.toStdTime();
}


double Sun::declination()
{
	JdTime time;
	double meanSolarNoonValue     = nextMeanSolarNoonJD(time);
	double meanSolarAnomalyValue  = meanSolarAnomaly(meanSolarNoonValue);
	return solarDeclination(eclipticLongitude(eqOfCenter(meanSolarAnomalyValue), meanSolarAnomalyValue));
}