Commit 53e4fee5 authored by Volker Krause's avatar Volker Krause

Add support for calculating civil twilight times

This is needed for example to implement the OSM opening hours specification
(https://wiki.openstreetmap.org/wiki/Key:opening_hours), as attempted in
https://invent.kde.org/vkrause/kopeninghours.

The existing sunrise code has been refactored to support this, and it
should now be very easy to add support for other sun position times such
as astronomical or nautical twilight as well, should that ever be needed.
parent c3e7aa49
......@@ -65,3 +65,48 @@ void SunriseTest::TestSunset()
QCOMPARE(utcSunset(QDate(2012, 7, 1), -14.60, 133.77), QTime(8, 47));
QCOMPARE(utcSunset(QDate(2012, 12, 31), -14.60, 133.77), QTime(9, 37));
}
void SunriseTest::TestDawn()
{
//NYC
QCOMPARE(utcDawn(QDate(2012, 1, 1), 40.72, -74.02), QTime(11, 49));
QCOMPARE(utcDawn(QDate(2012, 7, 1), 40.72, -74.02), QTime(8, 56));
QCOMPARE(utcDawn(QDate(2012, 12, 31), 40.72, -74.02), QTime(11, 49));
//LA
QCOMPARE(utcDawn(QDate(2012, 1, 1), 34.05, -118.23), QTime(14, 31));
QCOMPARE(utcDawn(QDate(2012, 7, 1), 34.05, -118.23), QTime(12, 16));
QCOMPARE(utcDawn(QDate(2012, 12, 31), 34.05, -118.23), QTime(14, 31));
//Berlin
QCOMPARE(utcDawn(QDate(2012, 1, 1), 51.40, 7.38), QTime(6, 56));
QCOMPARE(utcDawn(QDate(2012, 7, 1), 51.40, 7.38), QTime(2, 32));
QCOMPARE(utcDawn(QDate(2012, 12, 31), 51.40, 7.38), QTime(6, 56));
//Tasmania
QCOMPARE(utcDawn(QDate(2012, 1, 1), -14.60, 133.77), QTime(20, 16));
QCOMPARE(utcDawn(QDate(2012, 7, 1), -14.60, 133.77), QTime(21, 7));
QCOMPARE(utcDawn(QDate(2012, 12, 31), -14.60, 133.77), QTime(20, 15));
}
void SunriseTest::TestDusk()
{
QCOMPARE(utcDusk(QDate(2012, 1, 1), 40.72, -74.02), QTime(22, 10));
QCOMPARE(utcDusk(QDate(2012, 7, 1), 40.72, -74.02), QTime(1, 4));
QCOMPARE(utcDusk(QDate(2012, 12, 31), 40.72, -74.02), QTime(22, 9));
//LA
QCOMPARE(utcDusk(QDate(2012, 1, 1), 34.05, -118.23), QTime(1, 22));
QCOMPARE(utcDusk(QDate(2012, 7, 1), 34.05, -118.23), QTime(3, 37));
QCOMPARE(utcDusk(QDate(2012, 12, 31), 34.05, -118.23), QTime(1, 22));
//Berlin
QCOMPARE(utcDusk(QDate(2012, 1, 1), 51.40, 7.38), QTime(16, 12));
QCOMPARE(utcDusk(QDate(2012, 7, 1), 51.40, 7.38), QTime(20, 37));
QCOMPARE(utcDusk(QDate(2012, 12, 31), 51.40, 7.38), QTime(16, 12));
//Tasmania
QCOMPARE(utcDusk(QDate(2012, 1, 1), -14.60, 133.77), QTime(10, 1));
QCOMPARE(utcDusk(QDate(2012, 7, 1), -14.60, 133.77), QTime(9, 10));
QCOMPARE(utcDusk(QDate(2012, 12, 31), -14.60, 133.77), QTime(10, 1));
}
......@@ -17,6 +17,8 @@ class SunriseTest : public QObject
private Q_SLOTS:
void TestSunrise();
void TestSunset();
void TestDawn();
void TestDusk();
};
#endif
......@@ -144,18 +144,25 @@ static double calcEquationOfTime(double t)
return radToDeg(Etime) * 4.0; // in minutes of time
}
//calculate the hour angle of the sun at sunrise for the latitude (in radians)
static double calcHourAngleSunrise(double latitude, double solarDec)
// sun positions (in degree relative to the horizon) for certain events
static constexpr const double Sunrise = -0.833; // see https://en.wikipedia.org/wiki/Sunrise
static constexpr const double CivilTwilight = -6.0; // see https://en.wikipedia.org/wiki/Twilight
static constexpr const double Up = 1.0;
static constexpr const double Down = -1.0;
//calculate the hour angle of the sun at angle @p sunHeight for the latitude (in radians)
static double calcHourAngleSunrise(double latitude, double solarDec, double sunHeight)
{
double latRad = degToRad(latitude);
double sdRad = degToRad(solarDec);
double HAarg = (cos(degToRad(90.833)) /
double HAarg = (cos(degToRad(90.0 - sunHeight)) /
(cos(latRad) * cos(sdRad)) - tan(latRad) * tan(sdRad));
double HA = acos(HAarg);
return HA; // in radians (for sunset, use -HA)
}
QTime KHolidays::SunRiseSet::utcSunrise(const QDate &date, double latitude, double longitude)
static QTime calcSunEvent(const QDate &date, double latitude, double longitude, double sunHeight, double direction)
{
latitude = qMax(qMin(MaxLat, latitude), -MaxLat);
longitude = qMax(qMin(MaxLong, longitude), -MaxLong);
......@@ -163,7 +170,7 @@ QTime KHolidays::SunRiseSet::utcSunrise(const QDate &date, double latitude, doub
double t = calcTimeJulianCent(date.toJulianDay());
double eqTime = calcEquationOfTime(t);
double solarDec = calcSunDeclination(t);
double hourAngle = calcHourAngleSunrise(latitude, solarDec);
double hourAngle = direction * calcHourAngleSunrise(latitude, solarDec, sunHeight);
double delta = longitude + radToDeg(hourAngle);
QTime timeUTC(0, 0);
if (std::isnan(delta)) {
......@@ -175,22 +182,22 @@ QTime KHolidays::SunRiseSet::utcSunrise(const QDate &date, double latitude, doub
0);
}
QTime KHolidays::SunRiseSet::utcSunrise(const QDate &date, double latitude, double longitude)
{
return calcSunEvent(date, latitude, longitude, Sunrise, Up);
}
QTime KHolidays::SunRiseSet::utcSunset(const QDate &date, double latitude, double longitude)
{
latitude = qMax(qMin(MaxLat, latitude), -MaxLat);
longitude = qMax(qMin(MaxLong, longitude), -MaxLong);
return calcSunEvent(date, latitude, longitude, Sunrise, Down);
}
double t = calcTimeJulianCent(date.toJulianDay());
double eqTime = calcEquationOfTime(t);
double solarDec = calcSunDeclination(t);
double hourAngle = -calcHourAngleSunrise(latitude, solarDec);
double delta = longitude + radToDeg(hourAngle);
QTime timeUTC(0, 0);
if (std::isnan(delta)) {
return timeUTC;
}
timeUTC = timeUTC.addSecs((720 - (4.0 * delta) - eqTime) * 60);
return QTime(timeUTC.hour(),
timeUTC.second() > 29 ? timeUTC.minute() + 1 : timeUTC.minute(),
0);
QTime SunRiseSet::utcDawn(const QDate& date, double latitude, double longitude)
{
return calcSunEvent(date, latitude, longitude, CivilTwilight, Up);
}
QTime SunRiseSet::utcDusk(const QDate& date, double latitude, double longitude)
{
return calcSunEvent(date, latitude, longitude, CivilTwilight, Down);
}
......@@ -46,6 +46,32 @@ KHOLIDAYS_EXPORT QTime utcSunrise(const QDate &date, double latitude, double lon
*/
KHOLIDAYS_EXPORT QTime utcSunset(const QDate &date, double latitude, double longitude);
/**
Compute the civil dawn time (UTC) for a date and Earth location.
@param date is any valid QDate.
@param latitude is a floating point representing a valid latitude (-90.0, 90.0)
@param longitude is a floating point representing a valid longitude (-180.0, 180.0)
@see https://en.wikipedia.org/wiki/Twilight
@return the QTime of the sunrise in UTC.
@note the latitude and longitude are truncated as needed to fit into their proper range.
@since 5.77
*/
KHOLIDAYS_EXPORT QTime utcDawn(const QDate &date, double latitude, double longitude);
/**
Compute the civil dawn time (UTC) for a date and Earth location.
@param date is any valid QDate.
@param latitude is a floating point representing a valid latitude (-90.0, 90.0)
@param longitude is a floating point representing a valid longitude (-180.0, 180.0)
@see https://en.wikipedia.org/wiki/Twilight
@return the QTime of the sunset in UTC.
@note the latitude and longitude are truncated as needed to fit into their proper range.
@since 5.77
*/
KHOLIDAYS_EXPORT QTime utcDusk(const QDate &date, double latitude, double longitude);
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment