diff options
Diffstat (limited to 'src/solar.c')
| -rw-r--r-- | src/solar.c | 180 |
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; +} |