Sunrise and Sunset in C++

I wanted to calculate sunrise and sunset at my location for a recent project I was working on where I wouldn’t necessarily have access to internet connectivity. I did some quick searching and found a Wikipedia page with equations and thought it would be straightforward to put those equations into code. Somehow, I was not able to get the routines working correctly.

I came across a NOAA page that used equations in a different form, and also had links to an excel spreadsheet I could download and play with. I was able to reverse engineer the spreadsheet into C++ code to get a working routine that calculates the sunrise and sunset to an accuracy that is good enough for my needs.

I’m sharing my code here. If you see anything I’ve done wrong, please let me know. If you find it useful, please let me know about that as well.

/////////////////////////////////////////////////////////////////////////////
// From NOAA Spreadsheet https://gml.noaa.gov/grad/solcalc/calcdetails.html
bool getSunriseSunset(time_t& Sunrise, time_t& Sunset, const time_t& TheTime, const double Latitude, double Longitude)
{
	bool rval = false;
	struct tm LocalTime;
	if (0 != localtime_r(&TheTime, &LocalTime))
	{
		// if we don't have a valid latitude or longitude, declare sunrise to be midnight, and sunset one second before midnight
		if ((Latitude == 0) || (Longitude == 0))
		{
			LocalTime.tm_hour = 0;
			LocalTime.tm_min = 0;
			LocalTime.tm_sec = 0;
			Sunrise = mktime(&LocalTime);
			Sunset = Sunrise + 24*60*60 - 1;
		}
		else
		{
			double JulianDay = Time2JulianDate(TheTime); // F
			double JulianCentury = (JulianDay - 2451545) / 36525;	// G
			double GeomMeanLongSun = fmod(280.46646 + JulianCentury * (36000.76983 + JulianCentury * 0.0003032), 360);	// I
			double GeomMeanAnomSun = 357.52911 + JulianCentury * (35999.05029 - 0.0001537 * JulianCentury);	// J
			double EccentEarthOrbit = 0.016708634 - JulianCentury * (0.000042037 + 0.0000001267 * JulianCentury);	// K
			double SunEqOfCtr = sin(radians(GeomMeanAnomSun)) * (1.914602 - JulianCentury * (0.004817 + 0.000014 * JulianCentury)) + sin(radians(2 * GeomMeanAnomSun)) * (0.019993 - 0.000101 * JulianCentury) + sin(radians(3 * GeomMeanAnomSun)) * 0.000289; // L
			double SunTrueLong = GeomMeanLongSun + SunEqOfCtr;	// M
			double SunAppLong = SunTrueLong - 0.00569 - 0.00478 * sin(radians(125.04 - 1934.136 * JulianCentury));	// P
			double MeanObliqEcliptic = 23 + (26 + ((21.448 - JulianCentury * (46.815 + JulianCentury * (0.00059 - JulianCentury * 0.001813)))) / 60) / 60;	// Q
			double ObliqCorr = MeanObliqEcliptic + 0.00256 * cos(radians(125.04 - 1934.136 * JulianCentury));	// R
			double SunDeclin = degrees(asin(sin(radians(ObliqCorr)) * sin(radians(SunAppLong))));	// T
			double var_y = tan(radians(ObliqCorr / 2)) * tan(radians(ObliqCorr / 2));	// U
			double EquationOfTime = 4 * degrees(var_y * sin(2 * radians(GeomMeanLongSun)) - 2 * EccentEarthOrbit * sin(radians(GeomMeanAnomSun)) + 4 * EccentEarthOrbit * var_y * sin(radians(GeomMeanAnomSun)) * sin(2 * radians(GeomMeanLongSun)) - 0.5 * var_y * var_y * sin(4 * radians(GeomMeanLongSun)) - 1.25 * EccentEarthOrbit * EccentEarthOrbit * sin(2 * radians(GeomMeanAnomSun))); // V
			double HASunriseDeg = degrees(acos(cos(radians(90.833)) / (cos(radians(Latitude)) * cos(radians(SunDeclin))) - tan(radians(Latitude)) * tan(radians(SunDeclin)))); // W
			double SolarNoon = (720 - 4 * Longitude - EquationOfTime + LocalTime.tm_gmtoff / 60) / 1440; // X
			double SunriseTime = SolarNoon - HASunriseDeg * 4 / 1440;	// Y
			double SunsetTime = SolarNoon + HASunriseDeg * 4 / 1440;	// Z
			LocalTime.tm_hour = 0;
			LocalTime.tm_min = 0;
			LocalTime.tm_sec = 0;
			time_t Midnight = mktime(&LocalTime);
			Sunrise = Midnight + SunriseTime * 86400;
			Sunset = Midnight + SunsetTime * 86400;
		}
		rval = true;
	}
	return(rval);
}

I’ve included the individual routines I created when trying to duplicate the Wikipedia equations. If anyone can point out why they didn’t work for me, I’d appreciate that as well.

/////////////////////////////////////////////////////////////////////////////
double radians(const double degrees)
{
	return((degrees * M_PI) / 180.0);
}
double degrees(const double radians)
{
	return((radians * 180.0) / M_PI);
}
double Time2JulianDate(const time_t& TheTime)
{
	double JulianDay = 0;
	struct tm UTC;
	if (0 != gmtime_r(&TheTime, &UTC))
	{
		// https://en.wikipedia.org/wiki/Julian_day
		// JDN = (1461 × (Y + 4800 + (M ? 14)/12))/4 +(367 × (M ? 2 ? 12 × ((M ? 14)/12)))/12 ? (3 × ((Y + 4900 + (M - 14)/12)/100))/4 + D ? 32075
		JulianDay = (1461 * ((UTC.tm_year + 1900) + 4800 + ((UTC.tm_mon + 1) - 14) / 12)) / 4
			+ (367 * ((UTC.tm_mon + 1) - 2 - 12 * (((UTC.tm_mon + 1) - 14) / 12))) / 12
			- (3 * (((UTC.tm_year + 1900) + 4900 + ((UTC.tm_mon + 1) - 14) / 12) / 100)) / 4
			+ (UTC.tm_mday)
			- 32075;
		// JD = JDN + (hour-12)/24 + minute/1440 + second/86400
		double partialday = (static_cast<double>((UTC.tm_hour - 12)) / 24) + (static_cast<double>(UTC.tm_min) / 1440.0) + (static_cast<double>(UTC.tm_sec) / 86400.0);
		JulianDay += partialday;
	}
	return(JulianDay);
}
time_t JulianDate2Time(const double JulianDate)
{
	time_t TheTime = (JulianDate - 2440587.5) * 86400;
	return(TheTime);
}
double JulianDate2JulianDay(const double JulianDate)
{
	double n = JulianDate - 2451545.0 + 0.0008;
	return(n);
}
/////////////////////////////////////////////////////////////////////////////
// These equations all come from https://en.wikipedia.org/wiki/Sunrise_equation
double getMeanSolarTime(const double JulianDay, const double longitude)
{
	// an approximation of mean solar time expressed as a Julian day with the day fraction.
	double MeanSolarTime = JulianDay - (longitude / 360);
	return (MeanSolarTime);
}
double getSolarMeanAnomaly(const double MeanSolarTime)
{
	double SolarMeanAnomaly = fmod(357.5291 + 0.98560028 * MeanSolarTime, 360);
	return(SolarMeanAnomaly);
}
double getEquationOfTheCenter(const double SolarMeanAnomaly)
{
	double EquationOfTheCenter = 1.9148 * sin(radians(SolarMeanAnomaly)) + 0.0200 * sin(radians(2 * SolarMeanAnomaly)) + 0.0003 * sin(radians(3 * SolarMeanAnomaly));
	return(EquationOfTheCenter);
}
double getEclipticLongitude(const double SolarMeanAnomaly, const double EquationOfTheCenter)
{
	double EclipticLongitude = fmod(SolarMeanAnomaly + EquationOfTheCenter + 180 + 102.9372, 360);
	return(EclipticLongitude);
}
double getSolarTransit(const double MeanSolarTime, const double SolarMeanAnomaly, const double EclipticLongitude)
{
	// the Julian date for the local true solar transit (or solar noon).
	double SolarTransit = 2451545.0 + MeanSolarTime + 0.0053 * sin(radians(SolarMeanAnomaly)) - 0.0069 * sin(radians(2 * EclipticLongitude));
	return(SolarTransit);
}
double getDeclinationOfTheSun(const double EclipticLongitude)
{
	double DeclinationOfTheSun = sin(radians(EclipticLongitude)) * sin(radians(23.44));
	return(DeclinationOfTheSun);
}
double getHourAngle(const double Latitude, const double DeclinationOfTheSun)
{
	double HourAngle = (sin(radians(-0.83)) - sin(radians(Latitude)) * sin(radians(DeclinationOfTheSun))) / (cos(radians(Latitude)) * cos(radians(DeclinationOfTheSun)));
	return(HourAngle);
}
double getSunrise(const double SolarTransit, const double HourAngle)
{
	double Sunrise = SolarTransit - (HourAngle / 360);
	return(Sunrise);
}
double getSunset(const double SolarTransit, const double HourAngle)
{
	double Sunset = SolarTransit + (HourAngle / 360);
	return(Sunset);
}

Here’s the code snippet I was using when trying to get the Wikipedia equations working.

double JulianDay = JulianDate2JulianDay(Time2JulianDate(LoopStartTime));
double MeanSolarTime = getMeanSolarTime(JulianDay, Longitude);
double SolarMeanAnomaly = getSolarMeanAnomaly(MeanSolarTime);
double EquationOfTheCenter = getEquationOfTheCenter(SolarMeanAnomaly);
double EclipticLongitude = getEclipticLongitude(SolarMeanAnomaly, EquationOfTheCenter);
double SolarTransit = getSolarTransit(MeanSolarTime, SolarMeanAnomaly, EclipticLongitude);
double DeclinationOfTheSun = getDeclinationOfTheSun(EclipticLongitude);
double HourAngle = getHourAngle(Latitude, DeclinationOfTheSun);
double Sunrise = getSunrise(SolarTransit, HourAngle);
double Sunset = getSunset(SolarTransit, HourAngle);
std::cout.precision(std::numeric_limits<double>::max_digits10);
std::cout << "         Julian Date: " << Time2JulianDate(LoopStartTime) << std::endl;
std::cout << "           Unix Time: " << LoopStartTime << std::endl;
std::cout << "         Julian Date: " << timeToExcelLocal(JulianDate2Time(Time2JulianDate(LoopStartTime))) << std::endl;
std::cout << "            Latitude: " << Latitude << std::endl;
std::cout << "           Longitude: " << Longitude << std::endl;
std::cout << "          Julian Day: " << JulianDay << std::endl;
std::cout << "       MeanSolarTime: " << MeanSolarTime << std::endl;
std::cout << "    SolarMeanAnomaly: " << SolarMeanAnomaly << std::endl;
std::cout << " EquationOfTheCenter: " << EquationOfTheCenter << std::endl;
std::cout << "   EclipticLongitude: " << EclipticLongitude << std::endl;
std::cout << "        SolarTransit: " << timeToExcelLocal(JulianDate2Time(SolarTransit)) << std::endl;
std::cout << " DeclinationOfTheSun: " << DeclinationOfTheSun << std::endl;
std::cout << "           HourAngle: " << HourAngle << std::endl;
std::cout << "             Sunrise: " << timeToExcelLocal(JulianDate2Time(Sunrise)) << std::endl;
std::cout << "              Sunset: " << timeToExcelLocal(JulianDate2Time(Sunset)) << std::endl;

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s