/* See LICENSE file for copyright and license details. */ #include "libred.h" #include #include #include #include #include #include #if __GNUC__ # pragma GCC diagnostic push # pragma GCC diagnostic ignored "-Wunsuffixed-float-constants" #endif /* 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 /** * Get current Julian Centuries time (100 Julian Days since J2000) * and Julian Day time * * @param tc_out Output parameter for the current Julian Centuries time * @param td_out Output parameter for the current Julian Day time * @return 0 on success, -1 on failure * @throws Any error specified for clock_gettime(3) on error */ static int julian_time(double *tc_out, double *td_out) { struct timespec now; double tu; if (clock_gettime(CLOCK_REALTIME_COARSE, &now)) return -1; tu = fma((double)now.tv_nsec, 0.000000001, (double)now.tv_sec); *td_out = tu / 86400.0 + 2440587.5; *tc_out = (*td_out - 2451545.0) / 36525.0; return 0; } /** * Convert an angle (or otherwise) from degrees to radians * * @param deg The angle in degrees * @param The angle in radians */ static double radians(double deg) { return (double)M_PI / 180.0 * deg; } /** * Convert an angle (or otherwise) from radians to degrees * * @param rad The angle in radians * @param The angle in degrees */ static double degrees(double rad) { return 180.0 / (double)M_PI * rad; } /** * Convert an angle (or otherwise) from degrees to radians * and, using fused multply–add, add some number of degrees * * @param deg The angle in degrees * @param aug The number of radians to add * @param The angle in radians, plus `aug` */ static double radians_plus(double deg, double aug) { return fma((double)M_PI / 180.0, deg, aug); } /** * Convert an angle (or otherwise) from radians to degrees * and, using fused multply–add, add some number of degrees * * @param rad The angle in radians * @param aug The number of degrees to add * @param The angle in degrees, plus `aug` */ static double degrees_plus(double rad, double aug) { return fma(180.0 / (double)M_PI, rad, aug); } /** * Calculates the Sun's elevation from the solar hour angle * * @param latitude The latitude in degrees northwards from * the equator, negative for southwards * @param declination The declination, in radians * @param hour_angle The solar hour angle, in radians * @return The Sun's elevation, in radians */ static double elevation_from_hour_angle(double latitude, double declination, double hour_angle) { double c, s; latitude = radians(latitude); c = cos(latitude) * cos(declination); s = sin(latitude) * sin(declination); return asin(fma(c, cos(hour_angle), s)); } /** * Calculates the Sun's geometric mean longitude * * @param t The time in Julian Centuries * @return The Sun's geometric mean longitude in radians */ static double sun_geometric_mean_longitude(double t) { return radians(fmod(fma(fma(0.0003032, t, 36000.76983), t, 280.46646), 360.0)); } /** * Calculates the Sun's geometric mean anomaly * * @param t The time in Julian Centuries * @return The Sun's geometric mean anomaly in radians */ static double sun_geometric_mean_anomaly(double t) { return radians(fmod(fma(fma(-0.0001537, t, 35999.05029), t, 357.52911), 360.0)); } /** * Calculates the Earth's orbit eccentricity * * @param t The time in Julian Centuries * @return The Earth's orbit eccentricity */ static double earth_orbit_eccentricity(double t) { return fma(fma(-0.0000001267, t, -0.000042037), t, 0.016708634); } /** * Calculates the Sun's equation of the centre, the difference * between the true anomaly and the mean anomaly * * @param t The time in Julian Centuries * @return The Sun's equation of the centre, in radians */ static double sun_equation_of_centre(double t) { double a = sun_geometric_mean_anomaly(t), r; r = sin(1.0 * a) * fma(fma(-0.000014, t, 0.004817), t, 1.914602); r = fma(sin(2.0 * a), fma(-0.000101, t, 0.019993), r); r = fma(sin(3.0 * a), 0.000289, r); return radians(r); } /** * Calculates the Sun's real longitudinal position * * @param t The time in Julian Centuries * @return The longitude, in radians */ static double sun_real_longitude(double t) { return sun_geometric_mean_longitude(t) + sun_equation_of_centre(t); } /** * Calculates the Sun's apparent longitudinal position * * @param t The time in Julian Centuries * @return The longitude, in radians */ static double sun_apparent_longitude(double t) { double r = degrees_plus(sun_real_longitude(t), -0.00569); double a = radians(fma(-1934.136, t, 125.04)); return radians(fma(-0.00478, sin(a), r)); } /** * Calculates the mean ecliptic obliquity of the Sun's * apparent motion without variation correction * * @param t The time in Julian Centuries * @return The uncorrected mean obliquity, in radians */ static double mean_ecliptic_obliquity(double t) { double r = fma(fma(fma(0.001813, t, -0.00059), t, -46.815), t, 21.448); return radians(23.0 + (26.0 + r / 60.0) / 60.0); } /** * Calculates the mean ecliptic obliquity of the Sun's * parent motion with variation correction * * @param t The time in Julian Centuries * @return The mean obliquity, in radians */ static double corrected_mean_ecliptic_obliquity(double t) { double r = cos(radians(fma(-1934.136, t, 125.04))); return radians_plus(0.00256 * r, mean_ecliptic_obliquity(t)); } /** * Calculates the Sun's declination * * @param t The time in Julian Centuries * @return The Sun's declination, in radian */ static double solar_declination(double t) { double r = sin(corrected_mean_ecliptic_obliquity(t)); return asin(r * sin(sun_apparent_longitude(t))); } /** * Calculates the equation of time, the discrepancy * between apparent and mean solar time * * @param t The time in Julian Centuries * @return The equation of time, in degrees */ static double equation_of_time(double t) { double l = sun_geometric_mean_longitude(t); double e = earth_orbit_eccentricity(t); double m = sun_geometric_mean_anomaly(t); double y = pow(tan(corrected_mean_ecliptic_obliquity(t) / 2.0), 2.0); double r, c, s; s = y * sin(2.0 * l); c = y * cos(2.0 * l); r = fma(fma(4.0, c, -2.0), e * sin(m), s); r = fma(pow(0.50 * y, 2.0), sin(-4.0 * l), r); r = fma(pow(1.25 * e, 2.0), sin(-2.0 * m), r); return 4.0 * degrees(r); } /** * Calculates the Sun's elevation as apparent * from a geographical position * * @param tc The time in Julian Centuries * @param td The time in Julian Days * @param latitude The latitude in degrees northwards from * the equator, negative for southwards * @param longitude The longitude in degrees eastwards from * Greenwich, negative for westwards * @return The Sun's apparent elevation at the specified time as seen * from the specified position, measured in radians */ static double solar_elevation_from_time(double tc, double td, double latitude, double longitude) { double r; td = td - round(td); r = fma(1440, td - 1, -equation_of_time(tc)); r = radians(fma(0.25, r, -longitude)); return elevation_from_hour_angle(latitude, solar_declination(tc), r); } /** * Calculates the Sun's elevation as apparent * from a geographical position * * @param latitude The latitude in degrees northwards from * the equator, negative for southwards * @param longitude The longitude in degrees eastwards from * Greenwich, negative for westwards * @param elevation Output parameter for the Sun's apparent elevation * as seen, right now, from the specified position, * measured in degrees * @return 0 on success, -1 on failure * @throws Any error specified for clock_gettime(3) on error */ double libred_solar_elevation(double latitude, double longitude, double *elevation) { double tc, td; if (julian_time(&tc, &td)) return -1; *elevation = degrees(solar_elevation_from_time(tc, td, latitude, longitude)); return 0; } /** * This function is obsolete */ int libred_check_timetravel(void) { return 0; }