diff options
| -rw-r--r-- | src/solar.c | 206 | ||||
| -rw-r--r-- | src/solar.h | 32 |
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 */ |