aboutsummaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to '')
-rw-r--r--src/solar.c319
-rw-r--r--src/solar.h100
2 files changed, 419 insertions, 0 deletions
diff --git a/src/solar.c b/src/solar.c
new file mode 100644
index 0000000..38278a0
--- /dev/null
+++ b/src/solar.c
@@ -0,0 +1,319 @@
+/**
+ * Copyright © 2016 Mattias Andrée <maandree@member.fsf.org>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ */
+#include <math.h>
+#include <time.h>
+#include <errno.h>
+
+
+
+/**
+ * Get current Julian Centuries time (100 Julian days since J2000.)
+ *
+ * @return The current Julian Centuries time.
+ *
+ * @throws 0 On success.
+ * @throws Any error specified for clock_gettime(3) on error.
+ */
+static double
+julian_centuries()
+{
+ struct timespec now;
+ double tm;
+#if defined(CLOCK_REALTIME_COARSE)
+ if (clock_gettime(CLOCK_REALTIME_COARSE, &now))
+#else
+ if (clock_gettime(CLOCK_REALTIME, &now))
+#endif
+ return 0.0;
+ tm = (double)(now.tv_nsec);
+ tm /= 1000000000.0;
+ tm += (double)(now.tv_sec);
+ tm = (tm / 86400.0 + 2440587.5 - 2451545.0) / 36525.0;
+ return errno = 0, tm;
+}
+
+/**
+ * Convert a Julian Centuries timestamp to a Julian Day timestamp.
+ *
+ * @param tm The time in Julian Centuries
+ * @return The time in Julian Days
+ */
+static inline double
+julian_centuries_to_julian_day(double tm)
+{
+ return tm * 36525.0 + 2451545.0;
+}
+
+
+/**
+ * Convert an angle (or otherwise) from degrees to radians.
+ *
+ * @param deg The angle in degrees.
+ * @param The angle in radians.
+ */
+static inline double
+radians(double deg)
+{
+ return deg * (double)M_PI / 180.0;
+}
+
+/**
+ * Convert an angle (or otherwise) from radians to degrees.
+ *
+ * @param rad The angle in radians.
+ * @param The angle in degrees.
+ */
+static inline double
+degrees(double rad)
+{
+ return rad * 180.0 / (double)M_PI;
+}
+
+
+/**
+ * Calculates the Sun's elevation from the solar hour angle
+ *
+ * @param longitude The longitude in degrees eastwards.
+ * from Greenwich, negative for westwards.
+ * @param declination The declination, in radians.
+ * @param hour_angle The solar hour angle, in radians.
+ * @return The Sun's elevation, in radians.
+ */
+static inline double
+elevation_from_hour_angle(double latitude, double declination, double hour_angle)
+{
+ double rc = cos(radians(latitude));
+ rc *= cos(hour_angle) * cos(declination);
+ rc += sin(radians(latitude)) * sin(declination);
+ return asin(rc);
+}
+
+/**
+ * Calculates the Sun's geometric mean longitude.
+ *
+ * @param tm The time in Julian Centuries.
+ * @return The Sun's geometric mean longitude in radians.
+ */
+static inline double
+sun_geometric_mean_longitude(double tm)
+{
+ double rc = pow(0.0003032 * tm, 2.0) + 36000.76983 * tm + 280.46646;
+ return radians(fmod(rc, 360.0));
+ /*
+ CANNIBALISERS:
+ The result of this function should always be positive, this
+ means that after division modulo 360 but before `radians`,
+ you will need to add 360 if the value is negative. This can
+ only happen if `tm` is negative, which can only happen for date
+ times before 2000-(01)Jan-01 12:00:00 UTC par division modulo
+ implementations with the signess of at least the left operand.
+ More precively, it happens between circa 1970-(01)Jan-11
+ 16:09:02 UTC and circa 374702470660351740 seconds before
+ January 1, 1970 00:00 UTC, which is so far back in time
+ it cannot be reliable pinned down to the right year, but it
+ is without a shadow of a doubt looooong before the Earth
+ was formed, is right up there with the age of the Milky Way
+ and the universe itself.
+ */
+}
+
+/**
+ * Calculates the Sun's geometric mean anomaly.
+ *
+ * @param tm The time in Julian Centuries.
+ * @return The Sun's geometric mean anomaly in radians.
+ */
+static inline double
+sun_geometric_mean_anomaly(double tm)
+{
+ return radians(pow(-0.0001537 * tm, 2.0) + 35999.05029 * tm + 357.52911);
+}
+
+/**
+ * Calculates the Earth's orbit eccentricity.
+ *
+ * @param tm The time in Julian Centuries.
+ * @return The Earth's orbit eccentricity.
+ */
+static inline double
+earth_orbit_eccentricity(double tm)
+{
+ return pow(-0.0000001267 * tm, 2.0) - 0.000042037 * tm + 0.016708634;
+}
+
+/**
+ * Calculates the Sun's equation of the centre, the difference
+ * between the true anomaly and the mean anomaly.
+ *
+ * @param tm The time in Julian Centuries.
+ * @return The Sun's equation of the centre, in radians.
+ */
+static inline double
+sun_equation_of_centre(double tm)
+{
+ double a = sun_geometric_mean_anomaly(tm), rc;
+ rc = sin(1.0 * a) * (pow(-0.000014 * tm, 2.0) - 0.004817 * tm + 1.914602);
+ rc += sin(2.0 * a) * (-0.000101 * tm + 0.019993);
+ rc += sin(3.0 * a) * 0.000289;
+ return radians(rc);
+}
+
+/**
+ * Calculates the Sun's real longitudinal position.
+ *
+ * @param tm The time in Julian Centuries.
+ * @return The longitude, in radians.
+ */
+static inline double
+sun_real_longitude(double tm)
+{
+ double rc = sun_geometric_mean_longitude(tm);
+ return rc + sun_equation_of_centre(tm);
+}
+
+/**
+ * Calculates the Sun's apparent longitudinal position.
+ *
+ * @param tm The time in Julian Centuries.
+ * @return The longitude, in radians.
+ */
+static inline double
+sun_apparent_longitude(double tm)
+{
+ double rc = degrees(sun_real_longitude(tm)) - 0.00569;
+ rc -= 0.00478 * sin(radians(-1934.136 * tm + 125.04));
+ return radians(rc);
+}
+
+/**
+ * Calculates the mean ecliptic obliquity of the Sun's
+ * apparent motion without variation correction.
+ *
+ * @param tm The time in Julian Centuries.
+ * @return The uncorrected mean obliquity, in radians.
+ */
+static double
+mean_ecliptic_obliquity(double tm)
+{
+ double rc = pow(0.001813 * tm, 3.0) - pow(0.00059 * tm, 2.0) - 46.815 * tm + 21.448;
+ rc = 26 + rc / 60;
+ rc = 23 + rc / 60;
+ return radians(rc);
+}
+
+/**
+ * Calculates the mean ecliptic obliquity of the Sun's
+ * parent motion with variation correction.
+ *
+ * @param tm The time in Julian Centuries.
+ * @return The mean obliquity, in radians.
+ */
+static double
+corrected_mean_ecliptic_obliquity(double tm)
+{
+ double rc = -1934.136 * tm + 125.04;
+ rc = 0.00256 * cos(radians(rc));
+ rc += degrees(mean_ecliptic_obliquity(tm));
+ return radians(rc);
+}
+
+/**
+ * Calculates the Sun's declination.
+ *
+ * @param tm The time in Julian Centuries.
+ * @return The Sun's declination, in radian.
+ */
+static inline double
+solar_declination(double tm)
+{
+ double rc = sin(corrected_mean_ecliptic_obliquity(tm));
+ rc *= sin(sun_apparent_longitude(tm));
+ return asin(rc);
+}
+
+/**
+ * Calculates the equation of time, the discrepancy
+ * between apparent and mean solar time.
+ *
+ * @param tm The time in Julian Centuries.
+ * @return The equation of time, in degrees.
+ */
+static inline double
+equation_of_time(double tm)
+{
+ double l, e, m, y, rc;
+ l = sun_geometric_mean_longitude(tm);
+ e = earth_orbit_eccentricity(tm);
+ m = sun_geometric_mean_anomaly(tm);
+ y = corrected_mean_ecliptic_obliquity(tm);
+ y = pow(tan(y / 2.0), 2.0);
+ rc = y * sin(2.0 * l);
+ rc += (4.0 * y * cos(2.0 * l) - 2.0) * e * sin(m);
+ rc -= pow(0.5 * y, 2.0) * sin(4.0 * l);
+ rc -= pow(1.25 * e, 2.0) * sin(2.0 * m);
+ return 4.0 * degrees(rc);
+}
+
+/**
+ * Calculates the Sun's elevation as apparent.
+ * from a geographical position.
+ *
+ * @param tm The time in Julian Centuries.
+ * @param latitude The latitude in degrees northwards from
+ * the equator, negative for southwards.
+ * @param longitude The longitude in degrees eastwards from
+ * Greenwich, negative for westwards.
+ * @return The Sun's apparent elevation at the specified time as seen
+ * from the specified position, measured in radians.
+ */
+static inline double
+solar_elevation_from_time(double tm, double latitude, double longitude)
+{
+ double rc = julian_centuries_to_julian_day(tm);
+ rc = (rc - round(rc) - 0.5) * 1440;
+ rc = 720.0 - rc - equation_of_time(tm);
+ rc = radians(rc / 4.0 - longitude);
+ return elevation_from_hour_angle(latitude, solar_declination(tm), rc);
+}
+
+
+/**
+ * Calculates the Sun's elevation as apparent.
+ * from a geographical position.
+ *
+ * @param latitude The latitude in degrees northwards from
+ * the equator, negative for southwards.
+ * @param longitude The longitude in degrees eastwards from
+ * Greenwich, negative for westwards.
+ * @return The Sun's apparent elevation as seen, right now,
+ * from the specified position, measured in degrees.
+ *
+ * @throws 0 On success.
+ * @throws Any error specified for clock_gettime(3) on error.
+ */
+double
+solar_elevation(double latitude, double longitude)
+{
+ double tm = julian_centuries();
+ return errno ? -1 : degrees(solar_elevation_from_time(rm, latitude, longitude));
+}
+
diff --git a/src/solar.h b/src/solar.h
new file mode 100644
index 0000000..7eec736
--- /dev/null
+++ b/src/solar.h
@@ -0,0 +1,100 @@
+/**
+ * Copyright © 2016 Mattias Andrée <maandree@member.fsf.org>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ */
+
+
+/**
+ * Approximate apparent size of the Sun in degrees.
+ */
+#define SOLAR_APPARENT_RADIUS (32.0 / 60.0)
+
+
+/**
+ * The Sun's elevation at sunset and sunrise, measured in degrees.
+ */
+#define SOLAR_ELEVATION_SUNSET_SUNRISE (0.0)
+
+/**
+ * The Sun's elevation at civil dusk and civil dawn, measured in degrees.
+ */
+#define SOLAR_ELEVATION_CIVIL_DUSK_DAWN (-6.0)
+
+/**
+ * The Sun's elevation at nautical dusk and nautical dawn, measured in degrees.
+ */
+#define SOLAR_ELEVATION_NAUTICAL_DUSK_DAWN (-12.0)
+
+/**
+ * The Sun's elevation at astronomical dusk and astronomical dawn, measured in degrees.
+ */
+#define SOLAR_ELEVATION_ASTRONOMICAL_DUSK_DAWN (-18.0)
+
+
+/**
+ * Test whether it is twilight.
+ *
+ * @param ELEV:double The current elevation.
+ * @return 1 if is twilight, 0 otherwise.
+ */
+#define SOLAR_IS_TWILIGHT(ELEV) ((-18.0 <= (ELEV)) && ((ELEV) <= 0.0))
+
+/**
+ * Test whether it is civil twilight.
+ *
+ * @param ELEV:double The current elevation.
+ * @return 1 if is civil twilight, 0 otherwise.
+ */
+#define SOLAR_IS_CIVIL_TWILIGHT(ELEV) ((-6.0 <= (ELEV)) && ((ELEV) <= 0.0))
+
+/**
+ * Test whether it is nautical twilight.
+ *
+ * @param ELEV:double The current elevation.
+ * @return 1 if is nautical twilight, 0 otherwise.
+ */
+#define SOLAR_IS_NAUTICAL_TWILIGHT(ELEV) ((-12.0 <= (ELEV)) && ((ELEV) <= -6.0))
+
+/**
+ * Test whether it is astronomical twilight.
+ *
+ * @param ELEV:double The current elevation.
+ * @return 1 if is astronomical twilight, 0 otherwise.
+ */
+#define SOLAR_IS_ASTRONOMICAL_TWILIGHT(ELEV) ((-18.0 <= (ELEV)) && ((ELEV) <= -12.0))
+
+
+
+/**
+ * Calculates the Sun's elevation as apparent.
+ * from a geographical position.
+ *
+ * @param latitude The latitude in degrees northwards from
+ * the equator, negative for southwards.
+ * @param longitude The longitude in degrees eastwards from
+ * Greenwich, negative for westwards.
+ * @return The Sun's apparent elevation as seen, right now,
+ * from the specified position, measured in degrees.
+ *
+ * @throws 0 On success.
+ * @throws Any error specified for clock_gettime(3) on error.
+ */
+double solar_elevation(double latitude, double longitude);
+
generated by cgit v1.2.3-70-g09d2 (git 2.47.0) at 2024-12-01 10:12:01 +0000