aboutsummaryrefslogtreecommitdiffstats
path: root/solar.c
diff options
context:
space:
mode:
authorJon Lund Steffensen <jonlst@gmail.com>2009-12-23 17:33:17 +0100
committerJon Lund Steffensen <jonlst@gmail.com>2009-12-23 17:33:17 +0100
commit03dcd5d78a8cd5a9706f174f32cdfe8649272904 (patch)
tree2e84bc0e471d47c9ad1c3b1ca2d8e07d3ff65b59 /solar.c
parentMove RandR code to separate file. (diff)
downloadredshift-ng-03dcd5d78a8cd5a9706f174f32cdfe8649272904.tar.gz
redshift-ng-03dcd5d78a8cd5a9706f174f32cdfe8649272904.tar.bz2
redshift-ng-03dcd5d78a8cd5a9706f174f32cdfe8649272904.tar.xz
Move source and headers to src dir.
Diffstat (limited to 'solar.c')
-rw-r--r--solar.c318
1 files changed, 0 insertions, 318 deletions
diff --git a/solar.c b/solar.c
deleted file mode 100644
index 2d66fff..0000000
--- a/solar.c
+++ /dev/null
@@ -1,318 +0,0 @@
-/* solar.c -- Solar position source
- This file is part of Redshift.
-
- Redshift 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, either version 3 of the License, or
- (at your option) any later version.
-
- Redshift 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 General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with Redshift. If not, see <http://www.gnu.org/licenses/>.
-
- Copyright (c) 2009 Jon Lund Steffensen <jonlst@gmail.com>
-*/
-
-/* Ported from javascript code by U.S. Department of Commerce,
- National Oceanic & Atmospheric Administration:
- http://www.srrb.noaa.gov/highlights/sunrise/calcdetails.html
- It is based on equations from "Astronomical Algorithms" by
- Jean Meeus. */
-
-#include <math.h>
-#include <time.h>
-
-#include "solar.h"
-
-#define RAD(x) ((x)*(M_PI/180))
-#define DEG(x) ((x)*(180/M_PI))
-
-
-/* Angels of various times of day. */
-static const double time_angle[] = {
- [SOLAR_TIME_ASTRO_DAWN] = RAD(-90.0 + SOLAR_ASTRO_TWILIGHT_ELEV),
- [SOLAR_TIME_NAUT_DAWN] = RAD(-90.0 + SOLAR_NAUT_TWILIGHT_ELEV),
- [SOLAR_TIME_CIVIL_DAWN] = RAD(-90.0 + SOLAR_CIVIL_TWILIGHT_ELEV),
- [SOLAR_TIME_SUNRISE] = RAD(-90.0 + SOLAR_DAYTIME_ELEV),
- [SOLAR_TIME_NOON] = RAD(0.0),
- [SOLAR_TIME_SUNSET] = RAD(90.0 - SOLAR_DAYTIME_ELEV),
- [SOLAR_TIME_CIVIL_DUSK] = RAD(90.0 - SOLAR_CIVIL_TWILIGHT_ELEV),
- [SOLAR_TIME_NAUT_DUSK] = RAD(90.0 - SOLAR_NAUT_TWILIGHT_ELEV),
- [SOLAR_TIME_ASTRO_DUSK] = RAD(90.0 - SOLAR_ASTRO_TWILIGHT_ELEV)
-};
-
-
-/* Unix time from Julian day */
-static time_t
-unix_time_from_jd(double jd)
-{
- return 86400.0*(jd - 2440587.5);
-}
-
-/* Julian day from unix time */
-static double
-jd_from_unix_time(time_t t)
-{
- return (t / 86400.0) + 2440587.5;
-}
-
-/* Julian centuries since J2000.0 from Julian day */
-static double
-jcent_from_jd(double jd)
-{
- return (jd - 2451545.0) / 36525.0;
-}
-
-/* Julian day from Julian centuries since J2000.0 */
-static double
-jd_from_jcent(double t)
-{
- return 36525.0*t + 2451545.0;
-}
-
-/* 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));
-}
-
-/* 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));
-}
-
-/* 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);
-}
-
-/* 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);
-}
-
-/* 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));
-}
-
-/* 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);
-}
-
-/* Hour angle at the location for the given angular elevation.
- lat: Latitude of location in degrees
- decl: Declination in radians
- elev: Angular elevation angle in radians
- Return: Hour angle in radians */
-static double
-hour_angle_from_elevation(double lat, double decl, double elev)
-{
- double omega = acos((cos(fabs(elev)) - sin(RAD(lat))*sin(decl))/
- (cos(RAD(lat))*cos(decl)));
- return copysign(omega, -elev);
-}
-
-/* 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));
-}
-
-/* Time of apparent solar noon of location on earth.
- t: Julian centuries since J2000.0
- lon: Longitude of location in degrees
- Return: Time difference from mean solar midnigth in minutes */
-static double
-time_of_solar_noon(double t, double lon)
-{
- /* First pass uses approximate solar noon to
- calculate equation of time. */
- double t_noon = jcent_from_jd(jd_from_jcent(t) - lon/360.0);
- double eq_time = equation_of_time(t_noon);
- double sol_noon = 720 - 4*lon - eq_time;
-
- /* Recalculate using new solar noon. */
- t_noon = jcent_from_jd(jd_from_jcent(t) - 0.5 + sol_noon/1440.0);
- eq_time = equation_of_time(t_noon);
- sol_noon = 720 - 4*lon - eq_time;
-
- /* No need to do more iterations */
- return sol_noon;
-}
-
-/* Time of given apparent solar angular elevation of location on earth.
- t: Julian centuries since J2000.0
- t_noon: Apparent solar noon in Julian centuries since J2000.0
- lat: Latitude of location in degrees
- lon: Longtitude of location in degrees
- elev: Solar angular elevation in radians
- Return: Time difference from mean solar midnight in minutes */
-static double
-time_of_solar_elevation(double t, double t_noon,
- double lat, double lon, double elev)
-{
- /* First pass uses approximate sunrise to
- calculate equation of time. */
- double eq_time = equation_of_time(t_noon);
- double sol_decl = solar_declination(t_noon);
- double ha = hour_angle_from_elevation(lat, sol_decl, elev);
- double sol_offset = 720 - 4*(lon + DEG(ha)) - eq_time;
-
- /* Recalculate using new sunrise. */
- double t_rise = jcent_from_jd(jd_from_jcent(t) + sol_offset/1440.0);
- eq_time = equation_of_time(t_rise);
- sol_decl = solar_declination(t_rise);
- ha = hour_angle_from_elevation(lat, sol_decl, elev);
- sol_offset = 720 - 4*(lon + DEG(ha)) - eq_time;
-
- /* No need to do more iterations */
- return sol_offset;
-}
-
-/* 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);
-}
-
-double
-solar_elevation(time_t date, double lat, double lon)
-{
- double jd = jd_from_unix_time(date);
- return DEG(solar_elevation_from_time(jcent_from_jd(jd), lat, lon));
-}
-
-void
-solar_table_fill(time_t date, double lat, double lon, time_t *table)
-{
- /* Calculate Julian day */
- double jd = jd_from_unix_time(date);
-
- /* Calculate Julian day number */
- double jdn = round(jd);
- double t = jcent_from_jd(jdn);
-
- /* Calculate apparent solar noon */
- double sol_noon = time_of_solar_noon(t, lon);
- double j_noon = jdn - 0.5 + sol_noon/1440.0;
- double t_noon = jcent_from_jd(j_noon);
- table[SOLAR_TIME_NOON] = unix_time_from_jd(j_noon);
-
- /* Calculate solar midnight */
- table[SOLAR_TIME_MIDNIGHT] = unix_time_from_jd(j_noon + 0.5);
-
- /* Calulate absoute time of other phenomena */
- for (int i = 2; i < SOLAR_TIME_MAX; i++) {
- double angle = time_angle[i];
- double offset =
- time_of_solar_elevation(t, t_noon, lat, lon, angle);
- table[i] = unix_time_from_jd(jdn - 0.5 + offset/1440.0);
- }
-}