aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/solar.c206
-rw-r--r--src/solar.h32
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 */