src/solar.*: Add solar elevation calculation from Redshift

2 files changed, 238 insertions, 0 deletions

diff --git a/src/solar.c b/src/solar.c
new file mode 100644
index 0000000..d0f9d27
--- /dev/null
+++ b/src/solar.c
@@ -0,0 +1,206 @@
+/*
+ * Parts of this file have been taken from the Redshift project:
+ * https://github.com/jonls/redshift/
+ *
+ * Copyright 2010 Jon Lund Steffensen
+ * Copyright 2023 Robert Tari
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; version 3.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * Authors:
+ *   Jon Lund Steffensen <jonlst@gmail.com>
+ *   Robert Tari <robert@tari.in>
+ */
+
+#include <math.h>
+#include "solar.h"
+
+#define RAD(x)  ((x)*(M_PI/180))
+#define DEG(x)  ((x)*(180/M_PI))
+
+/* Julian centuries since J2000.0 from Julian day */
+static double
+jcent_from_jd(double jd)
+{
+    return (jd - 2451545.0) / 36525.0;
+}
+
+/* Angular elevation at the location for the given hour angle.
+   lat: Latitude of location in degrees
+   decl: Declination in radians
+   ha: Hour angle in radians
+   Return: Angular elevation in radians */
+static double
+elevation_from_hour_angle(double lat, double decl, double ha)
+{
+    return asin(cos(ha)*cos(RAD(lat))*cos(decl) +
+            sin(RAD(lat))*sin(decl));
+}
+
+/* Geometric mean anomaly of the sun.
+   t: Julian centuries since J2000.0
+   Return: Geometric mean anomaly in radians. */
+static double
+sun_geom_mean_anomaly(double t)
+{
+    return RAD(357.52911 + t*(35999.05029 - t*0.0001537));
+}
+
+/* Equation of center of the sun.
+   t: Julian centuries since J2000.0
+   Return: Center(?) in radians */
+static double
+sun_equation_of_center(double t)
+{
+    /* Use the first three terms of the equation. */
+    double m = sun_geom_mean_anomaly(t);
+    double c = sin(m)*(1.914602 - t*(0.004817 + 0.000014*t)) +
+        sin(2*m)*(0.019993 - 0.000101*t) +
+        sin(3*m)*0.000289;
+    return RAD(c);
+}
+
+/* Geometric mean longitude of the sun.
+   t: Julian centuries since J2000.0
+   Return: Geometric mean logitude in radians. */
+static double
+sun_geom_mean_lon(double t)
+{
+    /* FIXME returned value should always be positive */
+    return RAD(fmod(280.46646 + t*(36000.76983 + t*0.0003032), 360));
+}
+
+/* True longitude of the sun.
+   t: Julian centuries since J2000.0
+   Return: True longitude in radians */
+static double
+sun_true_lon(double t)
+{
+    double l_0 = sun_geom_mean_lon(t);
+    double c = sun_equation_of_center(t);
+    return l_0 + c;
+}
+
+/* Apparent longitude of the sun. (Right ascension).
+   t: Julian centuries since J2000.0
+   Return: Apparent longitude in radians */
+static double
+sun_apparent_lon(double t)
+{
+    double o = sun_true_lon(t);
+    return RAD(DEG(o) - 0.00569 - 0.00478*sin(RAD(125.04 - 1934.136*t)));
+}
+
+/* Mean obliquity of the ecliptic
+   t: Julian centuries since J2000.0
+   Return: Mean obliquity in radians */
+static double
+mean_ecliptic_obliquity(double t)
+{
+    double sec = 21.448 - t*(46.815 + t*(0.00059 - t*0.001813));
+    return RAD(23.0 + (26.0 + (sec/60.0))/60.0);
+}
+
+/* Corrected obliquity of the ecliptic.
+   t: Julian centuries since J2000.0
+   Return: Currected obliquity in radians */
+static double
+obliquity_corr(double t)
+{
+    double e_0 = mean_ecliptic_obliquity(t);
+    double omega = 125.04 - t*1934.136;
+    return RAD(DEG(e_0) + 0.00256*cos(RAD(omega)));
+}
+
+/* Declination of the sun.
+   t: Julian centuries since J2000.0
+   Return: Declination in radians */
+static double
+solar_declination(double t)
+{
+    double e = obliquity_corr(t);
+    double lambda = sun_apparent_lon(t);
+    return asin(sin(e)*sin(lambda));
+}
+
+/* Eccentricity of earth orbit.
+   t: Julian centuries since J2000.0
+   Return: Eccentricity (unitless). */
+static double
+earth_orbit_eccentricity(double t)
+{
+    return 0.016708634 - t*(0.000042037 + t*0.0000001267);
+}
+
+/* Difference between true solar time and mean solar time.
+   t: Julian centuries since J2000.0
+   Return: Difference in minutes */
+static double
+equation_of_time(double t)
+{
+    double epsilon = obliquity_corr(t);
+    double l_0 = sun_geom_mean_lon(t);
+    double e = earth_orbit_eccentricity(t);
+    double m = sun_geom_mean_anomaly(t);
+    double y = pow(tan(epsilon/2.0), 2.0);
+
+    double eq_time = y*sin(2*l_0) - 2*e*sin(m) +
+        4*e*y*sin(m)*cos(2*l_0) -
+        0.5*y*y*sin(4*l_0) -
+        1.25*e*e*sin(2*m);
+    return 4*DEG(eq_time);
+}
+
+/* Julian day from Julian centuries since J2000.0 */
+static double
+jd_from_jcent(double t)
+{
+    return 36525.0*t + 2451545.0;
+}
+
+/* Solar angular elevation at the given location and time.
+   t: Julian centuries since J2000.0
+   lat: Latitude of location
+   lon: Longitude of location
+   Return: Solar angular elevation in radians */
+static double
+solar_elevation_from_time(double t, double lat, double lon)
+{
+    /* Minutes from midnight */
+    double jd = jd_from_jcent(t);
+    double offset = (jd - round(jd) - 0.5)*1440.0;
+
+    double eq_time = equation_of_time(t);
+    double ha = RAD((720 - offset - eq_time)/4 - lon);
+    double decl = solar_declination(t);
+    return elevation_from_hour_angle(lat, decl, ha);
+}
+
+/* Julian day from unix epoch */
+static double
+jd_from_epoch(double t)
+{
+    return (t / 86400.0) + 2440587.5;
+}
+
+/* Solar angular elevation at the given location and time.
+   date: Seconds since unix epoch
+   lat: Latitude of location
+   lon: Longitude of location
+   Return: Solar angular elevation in degrees */
+double
+solar_elevation(double date, double lat, double lon)
+{
+    double jd = jd_from_epoch(date);
+    return DEG(solar_elevation_from_time(jcent_from_jd(jd), lat, lon));
+}
diff --git a/src/solar.h b/src/solar.h
new file mode 100644
index 0000000..b1faed9
--- /dev/null
+++ b/src/solar.h
@@ -0,0 +1,32 @@
+/*
+ * Parts of this file have been taken from the Redshift project:
+ * https://github.com/jonls/redshift/
+ *
+ * Copyright 2010 Jon Lund Steffensen
+ * Copyright 2023 Robert Tari
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; version 3.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * Authors:
+ *   Jon Lund Steffensen <jonlst@gmail.com>
+ *   Robert Tari <robert@tari.in>
+ */
+
+#ifndef REDSHIFT_SOLAR_H
+#define REDSHIFT_SOLAR_H
+
+#define SOLAR_CIVIL_TWILIGHT_ELEV    -6.0
+
+double solar_elevation(double date, double lat, double lon);
+
+#endif /* ! REDSHIFT_SOLAR_H */