aboutsummaryrefslogtreecommitdiffstats
path: root/src/solar.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/solar.c')
-rw-r--r--src/solar.c180
1 files changed, 85 insertions, 95 deletions
diff --git a/src/solar.c b/src/solar.c
index 8a57507..697ecd6 100644
--- a/src/solar.c
+++ b/src/solar.c
@@ -13,14 +13,21 @@
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ * This file is a stripped down port of <https://github.com/maandree/solar-python>.
*/
+#include "libred.h"
#include <math.h>
#include <time.h>
#include <errno.h>
-#if !defined(CLOCK_REALTIME_COARSE)
+/* Select clock. */
+#if defined(DO_NOT_USE_COARSEC_CLOCK) || !defined(CLOCK_REALTIME_COARSE)
+# ifdef CLOCK_REALTIME_COARSE
+# undef CLOCK_REALTIME_COARSE
+# endif
# define CLOCK_REALTIME_COARSE CLOCK_REALTIME
#endif
@@ -34,15 +41,14 @@
* @throws 0 On success.
* @throws Any error specified for clock_gettime(3) on error.
*/
-static double
-julian_centuries()
+static double julian_centuries()
{
- struct timespec now;
- double tm;
- if (clock_gettime(CLOCK_REALTIME_COARSE, &now)) return 0.0;
- tm = (double)(now.tv_nsec) / 1000000000.0 + (double)(now.tv_sec);
- tm = (tm / 86400.0 + 2440587.5 - 2451545.0) / 36525.0;
- return errno = 0, tm;
+ struct timespec now;
+ double tm;
+ if (clock_gettime(CLOCK_REALTIME_COARSE, &now)) return 0.0;
+ tm = (double)(now.tv_nsec) / 1000000000.0 + (double)(now.tv_sec);
+ tm = (tm / 86400.0 + 2440587.5 - 2451545.0) / 36525.0;
+ return errno = 0, tm;
}
/**
@@ -51,10 +57,9 @@ julian_centuries()
* @param tm The time in Julian Centuries
* @return The time in Julian Days
*/
-static inline double
-julian_centuries_to_julian_day(double tm)
+static inline double julian_centuries_to_julian_day(double tm)
{
- return tm * 36525.0 + 2451545.0;
+ return tm * 36525.0 + 2451545.0;
}
@@ -64,10 +69,9 @@ julian_centuries_to_julian_day(double tm)
* @param deg The angle in degrees.
* @param The angle in radians.
*/
-static inline double
-radians(double deg)
+static inline double radians(double deg)
{
- return deg * (double)M_PI / 180.0;
+ return deg * (double)M_PI / 180.0;
}
/**
@@ -76,10 +80,9 @@ radians(double deg)
* @param rad The angle in radians.
* @param The angle in degrees.
*/
-static inline double
-degrees(double rad)
+static inline double degrees(double rad)
{
- return rad * 180.0 / (double)M_PI;
+ return rad * 180.0 / (double)M_PI;
}
@@ -92,13 +95,12 @@ degrees(double rad)
* @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)
+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);
+ double rc = cos(radians(latitude));
+ rc *= cos(hour_angle) * cos(declination);
+ rc += sin(radians(latitude)) * sin(declination);
+ return asin(rc);
}
/**
@@ -107,14 +109,13 @@ elevation_from_hour_angle(double latitude, double declination, double hour_angle
* @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)
+static inline double sun_geometric_mean_longitude(double tm)
{
- double rc = fmod(pow(0.0003032 * tm, 2.0) + 36000.76983 * tm + 280.46646, 360.0);
+ double rc = fmod(pow(0.0003032 * tm, 2.0) + 36000.76983 * tm + 280.46646, 360.0);
#if defined(TIMETRAVELLER)
- rc = rc < 0.0 ? (rc + 360.0) : rc;
+ rc = rc < 0.0 ? (rc + 360.0) : rc;
#endif
- return radians(rc);
+ return radians(rc);
}
/**
@@ -123,10 +124,9 @@ sun_geometric_mean_longitude(double tm)
* @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)
+static inline double sun_geometric_mean_anomaly(double tm)
{
- return radians(pow(-0.0001537 * tm, 2.0) + 35999.05029 * tm + 357.52911);
+ return radians(pow(-0.0001537 * tm, 2.0) + 35999.05029 * tm + 357.52911);
}
/**
@@ -135,10 +135,9 @@ sun_geometric_mean_anomaly(double tm)
* @param tm The time in Julian Centuries.
* @return The Earth's orbit eccentricity.
*/
-static inline double
-earth_orbit_eccentricity(double tm)
+static inline double earth_orbit_eccentricity(double tm)
{
- return pow(-0.0000001267 * tm, 2.0) - 0.000042037 * tm + 0.016708634;
+ return pow(-0.0000001267 * tm, 2.0) - 0.000042037 * tm + 0.016708634;
}
/**
@@ -148,14 +147,13 @@ earth_orbit_eccentricity(double tm)
* @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)
+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);
+ 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);
}
/**
@@ -164,10 +162,9 @@ sun_equation_of_centre(double tm)
* @param tm The time in Julian Centuries.
* @return The longitude, in radians.
*/
-static inline double
-sun_real_longitude(double tm)
+static inline double sun_real_longitude(double tm)
{
- return sun_geometric_mean_longitude(tm) + sun_equation_of_centre(tm);
+ return sun_geometric_mean_longitude(tm) + sun_equation_of_centre(tm);
}
/**
@@ -176,11 +173,10 @@ sun_real_longitude(double tm)
* @param tm The time in Julian Centuries.
* @return The longitude, in radians.
*/
-static inline double
-sun_apparent_longitude(double tm)
+static inline double sun_apparent_longitude(double tm)
{
- double rc = degrees(sun_real_longitude(tm)) - 0.00569;
- return radians(rc - 0.00478 * sin(radians(-1934.136 * tm + 125.04)));
+ double rc = degrees(sun_real_longitude(tm)) - 0.00569;
+ return radians(rc - 0.00478 * sin(radians(-1934.136 * tm + 125.04)));
}
/**
@@ -190,11 +186,10 @@ sun_apparent_longitude(double tm)
* @param tm The time in Julian Centuries.
* @return The uncorrected mean obliquity, in radians.
*/
-static double
-mean_ecliptic_obliquity(double tm)
+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;
- return radians(23.0 + (26.0 + rc / 60.0) / 60.0);
+ double rc = pow(0.001813 * tm, 3.0) - pow(0.00059 * tm, 2.0) - 46.815 * tm + 21.448;
+ return radians(23.0 + (26.0 + rc / 60.0) / 60.0);
}
/**
@@ -204,11 +199,10 @@ mean_ecliptic_obliquity(double tm)
* @param tm The time in Julian Centuries.
* @return The mean obliquity, in radians.
*/
-static double
-corrected_mean_ecliptic_obliquity(double tm)
+static double corrected_mean_ecliptic_obliquity(double tm)
{
- double rc = 0.00256 * cos(radians(-1934.136 * tm + 125.04));
- return radians(rc + degrees(mean_ecliptic_obliquity(tm)));
+ double rc = 0.00256 * cos(radians(-1934.136 * tm + 125.04));
+ return radians(rc + degrees(mean_ecliptic_obliquity(tm)));
}
/**
@@ -217,11 +211,10 @@ corrected_mean_ecliptic_obliquity(double tm)
* @param tm The time in Julian Centuries.
* @return The Sun's declination, in radian.
*/
-static inline double
-solar_declination(double tm)
+static inline double solar_declination(double tm)
{
- double rc = sin(corrected_mean_ecliptic_obliquity(tm));
- return asin(rc * sin(sun_apparent_longitude(tm)));
+ double rc = sin(corrected_mean_ecliptic_obliquity(tm));
+ return asin(rc * sin(sun_apparent_longitude(tm)));
}
/**
@@ -231,18 +224,17 @@ solar_declination(double tm)
* @param tm The time in Julian Centuries.
* @return The equation of time, in degrees.
*/
-static inline double
-equation_of_time(double tm)
+static inline double equation_of_time(double tm)
{
- double l = sun_geometric_mean_longitude(tm);
- double e = earth_orbit_eccentricity(tm);
- double m = sun_geometric_mean_anomaly(tm);
- double y = pow(tan(corrected_mean_ecliptic_obliquity(tm) / 2.0), 2.0);
- double 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);
+ double l = sun_geometric_mean_longitude(tm);
+ double e = earth_orbit_eccentricity(tm);
+ double m = sun_geometric_mean_anomaly(tm);
+ double y = pow(tan(corrected_mean_ecliptic_obliquity(tm) / 2.0), 2.0);
+ double 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);
}
/**
@@ -257,14 +249,13 @@ equation_of_time(double tm)
* @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)
+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);
+ 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);
}
@@ -282,11 +273,10 @@ solar_elevation_from_time(double tm, double latitude, double longitude)
* @throws 0 On success.
* @throws Any error specified for clock_gettime(3) on error.
*/
-double
-solar_elevation(double latitude, double longitude)
+double libred_solar_elevation(double latitude, double longitude)
{
- double tm = julian_centuries();
- return errno ? -1 : degrees(solar_elevation_from_time(tm, latitude, longitude));
+ double tm = julian_centuries();
+ return errno ? -1 : degrees(solar_elevation_from_time(tm, latitude, longitude));
}
@@ -295,18 +285,18 @@ solar_elevation(double latitude, double longitude)
*
* @return 0 on success, -1 on error.
*/
-#if defined(TIMETRAVELLER)
-int
-check_timetravel(void)
+int libred_check_timetravel(void)
{
- struct timespec now;
- if (clock_gettime(CLOCK_REALTIME, &now)) return -1;
- if (now.tv_nsec < (time_t)946728000L)
- fprintf(stderr, "We have detected that you are a time-traveller"
- "(or your clock is not configured correctly.)"
- "Please recompile with -DTIMETRAVELLER"
- "(or correct your clock.)"), exit(1);
- return 0;
-}
+#if !defined(TIMETRAVELLER)
+ struct timespec now;
+ if (clock_gettime(CLOCK_REALTIME, &now)) return -1;
+ if (now.tv_nsec < (time_t)946728000L)
+ fprintf(stderr,
+ "We have detected that you are a time-traveller"
+ "(or your clock is not configured correctly.)"
+ "Please recompile libred with -DTIMETRAVELLER"
+ "(or correct your clock.)"), exit(1);
#endif
+ return 0;
+}