Skip to content

Instantly share code, notes, and snippets.

@DouglasAllen
Last active February 11, 2017 03:37
Show Gist options
  • Save DouglasAllen/0fb978595c9c7c65e649d09e29730471 to your computer and use it in GitHub Desktop.
Save DouglasAllen/0fb978595c9c7c65e649d09e29730471 to your computer and use it in GitHub Desktop.
Qt Sunrise author unknown appears to be based upon the Wikipedia sunrise equation
#include <QtCore>
#include "sunrise.h"
Sunrise::Sunrise(double latitude_, double longitude_, double elevation_)
: latitude(latitude_)
, longitude(-longitude_)
, elevation(elevation_)
{
}
int Sunrise::julian(const QDate &d)
{
QDateTime dt(d, QTime(12, 1, 0, 0));
return dt.toUTC().date().toJulianDay();
}
QTime Sunrise::sunrise(const QDate &d)
{
if (!d.isValid())
return QTime();
double j = sunrise(julian(d));
return date(j).time();
}
QTime Sunrise::noon(const QDate &d)
{
if (!d.isValid())
return QTime();
double j = noon(julian(d));
return date(j).time();
}
QTime Sunrise::sunset(const QDate &d)
{
if (!d.isValid())
return QTime();
double j = sunset(julian(d));
return date(j).time();
}
QDateTime Sunrise::date(double julian)
{
if (julian < 0.0)
return QDateTime();
// The day number is the integer part of the date
int julianDays = qFloor(julian);
QDate d = QDate::fromJulianDay(julianDays);
// The fraction is the time of day
double julianMSecs = (julian - static_cast<double>(julianDays)) * 86400.0 * 1000.0;
// Julian days start at noon (12:00 UTC)
return QDateTime(d, QTime(12, 0, 0, 0), Qt::UTC).addMSecs(qRound(julianMSecs)).toLocalTime();
}
/*
* sunrise.h
*
* Sunrise Time Class
* written by Joerg Dentler
*
*/
//
// Version 1.0.0 14-12-19 initial design
#ifndef __SUNRISE_H__
#define __SUNRISE_H__
/////////////////////////////////////////////////////////////////////////////
// Permission is granted to anyone to use this software for any
// purpose and to redistribute it in any way, subject to the following
// restrictions:
//
// 1. The author is not responsible for the consequences of use of
// this software, no matter how awful, even if they arise
// from defects in it.
//
// 2. The origin of this software must not be misrepresented, either
// by explicit claim or by omission.
//
// 3. Altered versions must be plainly marked as such, and must not
// be misrepresented (by explicit claim or omission) as being
// the original software.
//
// 4. This notice must not be removed or altered.
/////////////////////////////////////////////////////////////////////////////
#include <QDateTime>
#include <math.h>
#include <stdlib.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
class Sunrise
{
public:
/**
* create sunrise calculator for a certain map point
*
* latitude
* north is +
* south is -
*
* longitude
* east is +
* west is -
*
* eleveation in meters above ground
*
* examples:
*
* Berlin
* 52.516N 13.271E
*
* Boa Vista
* 16.141N 22.904W
*
*/
Sunrise(double latitude_, double longitude_, double elevation_ = 0);
/*
* calculation functions
* param d is local date
* returns local time format
*/
QTime sunrise(const QDate &d);
QTime noon(const QDate &d);
QTime sunset(const QDate &d);
private:
// Convert radian angle to degrees
double dRadToDeg(double dAngleRad)
{
return (180.0 * dAngleRad / M_PI);
}
// Convert degree angle to radians
double dDegToRad(double dAngleDeg)
{
return (M_PI * dAngleDeg / 180.0);
}
double normalizedAngle(double a)
{
while (!(a < 360.0))
a -= 360.0;
while (!(a > -360.0))
a += 360.0;
return a;
}
double julianCycle(int iJulianDay)
{
double n = double(iJulianDay) - 2451545.0009 - longitude / 360.0;
return floor(n + 0.5);
}
double solarNoonApprox(int iJulianDay)
{
return 2451545.0009 + longitude / 360.0 + julianCycle(iJulianDay);
}
double solarMean(double noon)
{
double m = 357.5291 + 0.98560028 * (noon - 2451545.0);
return normalizedAngle(m);
}
double center(double mean)
{
double m1 = normalizedAngle(mean * 1.0);
double m2 = normalizedAngle(mean * 2.0);
double m3 = normalizedAngle(mean * 3.0);
m1 = sin(dDegToRad(m1));
m2 = sin(dDegToRad(m2));
m3 = sin(dDegToRad(m3));
return 1.9148*m1 + 0.0200*m2 + 0.0003*m3;
}
double gamma(double m)
{
double c = center(m);
return normalizedAngle(m + 102.9372 + c + 180.0);
}
double solarTransit(double noon)
{
double m = solarMean(noon);
double g = gamma(m);
m = sin(dDegToRad(m));
g = normalizedAngle(g * 2.0);
g = sin(dDegToRad(g));
return noon + 0.0053*m - 0.0069*g;
}
bool omega(double &o, double noon)
{
double m = solarMean(noon);
double g = gamma(m);
double l = dDegToRad(latitude);
double sd = sin(dDegToRad(g)) * sin(dDegToRad(23.45));
double cd = cos(asin(sd));
double el = 0;
if (elevation > 0.0) {
el = -2.076;
el *= sqrt(elevation) / 60.0;
}
o = sin(dDegToRad(-0.83 + el)) - sin(l)*sd;
o /= cos(l) * cd;
if (o > 1.0 || o < -1.0)
return false;
o = dRadToDeg(acos(o));
return true;
}
double sunset(int iJulianDay)
{
double noon = solarNoonApprox(iJulianDay);
double m = solarMean(noon);
double g = gamma(m);
m = sin(dDegToRad(m));
g = normalizedAngle(g * 2.0);
g = sin(dDegToRad(g));
double o;
if (!omega(o, noon))
return -1.0;
o += longitude;
o /= 360.0;
return 2451545.0009 + o + julianCycle(iJulianDay) + 0.0053*m - 0.0069*g;
}
double noon(int iJulianDay)
{
double noon = solarNoonApprox(iJulianDay);
return solarTransit(noon);
}
double sunrise(int iJulianDay)
{
double transit = noon(iJulianDay);
double ss = sunset(iJulianDay);
return ss < 0.0 ? -1.0 : transit - (ss - transit);
}
QDateTime date(double julian);
int julian(const QDate &d);
const double latitude;
const double longitude;
const double elevation;
};
#endif
#include "sunrise.h"
#include <iostream>
int main(int argc, char **argv)
{
QDate d;
if (argc != 2)
d = QDate::currentDate();
else
d = QDate::fromString(argv[1], Qt::ISODate);
using namespace std;
if (!d.isValid()) {
cerr << "Invalid date" << endl;
return 1;
}
cout << d.toString().toStdString() << endl;
Sunrise sunrise(41.943, -88.75); //
cout << "sunrise: " << sunrise.sunrise(d).toString().toStdString() << endl;
cout << "noon: " << sunrise.noon(d).toString().toStdString() << endl;
cout << "sunset: " << sunrise.sunset(d).toString().toStdString() << endl;
return 0;
}
/*
*
"Fr. Dez 19 2014"
./sunrise 2014-12-19
19. A: 08:10
U: 16:03
Sonnenaufgang 08:12
Transit 12:06
Sonnenuntergang 16:01
Algorithm
Calculating Sunrise and Sunset
So you want to calculate sunrise and sunset. Simple (maybe).
If this doesn't make sense to you without real numbers, you might want to check out my in your face sunrise/sunset calculator.
If you want to know where all of the constants come from, check out the source of this information.
Variables:
Jdate = Julian date
lw = West Longitude (75W = 75, 45E = -45)
ln = North Latitude (35N = 35, 25S = -25)
M = Mean Solar Anomaly
n = Julian cycle since Jan 1, 2000
C = Equation of center
λ = Ecliptical longitude of the sun
δ = Declination of the sun
H = Hour Angle (half the arc length of the sun)
Jtransit = Julian date of solar noon on cycle n
Jrise = Julian date of sunrise on cycle n
Jset = Julian date of sunset on cycle n
First, start by calculating the number of days since January 1, 2000. Add that number to 2451545 (the Julian day of January 1, 2000). This will be variable Jdate.
The next step is to calculate the Julian cycle. This is not equal to the days since Jan 1, 2000. Depending on your longitude, this may be a different number.
n* = (Jdate - 2451545 - 0.0009) - (lw/360)
n = round(n*)
Now, it is time to approximate the Julian date of solar noon. This is just an approximation so that we can make some intermediate calculations before we calculate the actual Julian date of solar noon.
J* = 2451545 + 0.0009 + (lw/360) + n
Using the approximate value, calculate the mean solar anomaly. This will get a very close value to the actual mean solar anomaly.
M = [357.5291 + 0.98560028 * (J* - 2451545)] mod 360
Calculate the equation of center
C = (1.9148 * sin(M)) + (0.0200 * sin(2 * M)) + (0.0003 * sin(3 * M))
Now, using C and M, calculate the ecliptical longitude of the sun.
λ = (M + 102.9372 + C + 180) mod 360
Now there is enough data to calculate an accurate Julian date for solar noon.
Jtransit = J* + (0.0053 * sin(M)) - (0.0069 * sin(2 * λ))
To calculate the hour angle we need to find the declination of the sun
δ = arcsin( sin(λ) * sin(23.45) )
Now, calculate the hour angle, which corresponds to half of the arc length of the sun at this latitude at this declination of the sun
H = arccos( [sin(-0.83) - sin(ln) * sin(δ)] / [cos(ln) * cos(δ)] )
Note: If H is undefined, then there is either no sunrise (in winter) or no sunset (in summer) for the supplied latitude.
Okay, time to go back through the approximation again, this time we use H in the calculation
J** = 2451545 + 0.0009 + ((H + lw)/360) + n
The values of M and λ from above don't really change from solar noon to sunset, so there is no need to recalculate them before calculating sunset.
Jset = J** + (0.0053 * sin(M)) - (0.0069 * sin(2 * λ))
Instead of going through that mess again, assume that solar noon is half-way between sunrise and sunset (valid for latitudes < 60) and approximate sunrise.
Jrise = Jtransit - (Jset - Jtransit)
At this point you should convert the Julian dates to something more readable like a regular date. This will be left as an exercise to the reader.
Source: http://aa.quae.nl/en/reken/zonpositie.html
*/
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment