physics/astronomy.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // SPDX-License-Identifier: BSD-3-Clause | ||
| 2 | /// @file astronomy.cpp | ||
| 3 | /// @brief Synthetic stars, stellar photometry and sky orientation. | ||
| 4 | #include "ernie/physics/astronomy.h" | ||
| 5 | |||
| 6 | #include "ernie/physics/spectrum.h" | ||
| 7 | |||
| 8 | #include <glm/gtc/constants.hpp> | ||
| 9 | |||
| 10 | #include <algorithm> | ||
| 11 | #include <array> | ||
| 12 | #include <cmath> | ||
| 13 | |||
| 14 | namespace ernie { | ||
| 15 | |||
| 16 | namespace { | ||
| 17 | |||
| 18 | constexpr double kMinKelvin = 1000.0; | ||
| 19 | constexpr double kMaxKelvin = 50000.0; | ||
| 20 | constexpr std::size_t kTableSize = 256; | ||
| 21 | |||
| 22 | /// The illuminance in lux of a black body at `kelvin` whose spectral | ||
| 23 | /// irradiance at 555 nm is 1 W m^-2 nm^-1, by the midpoint rule at 1 | ||
| 24 | /// nm over 360 to 830 nm. | ||
| 25 | 256 | double LuxPerUnitAt555(double kelvin) { | |
| 26 | 256 | const double at555 = PlanckRadiance(555.0, kelvin); | |
| 27 | 256 | double sum = 0.0; | |
| 28 | 120576 | for (double nm = 360.5; nm < 830.0; nm += 1.0) { | |
| 29 | 120320 | sum += PlanckRadiance(nm, kelvin) / at555 * CieXyz(nm).y; | |
| 30 | } | ||
| 31 | 256 | return kLuminousEfficacy * sum; | |
| 32 | } | ||
| 33 | |||
| 34 | /// `LuxPerUnitAt555` at temperatures spaced evenly in log temperature | ||
| 35 | /// over 1000 to 50000 K. | ||
| 36 | 1880 | const std::array<double, kTableSize>& LuxTable() { | |
| 37 | 1 | static const std::array<double, kTableSize> table = [] { | |
| 38 | 1 | std::array<double, kTableSize> values{}; | |
| 39 | 257 | for (std::size_t i = 0; i < kTableSize; ++i) { | |
| 40 | 256 | const double t = static_cast<double>(i) / (kTableSize - 1); | |
| 41 | 256 | values[i] = LuxPerUnitAt555(kMinKelvin * std::pow(kMaxKelvin / kMinKelvin, t)); | |
| 42 | } | ||
| 43 | 1 | return values; | |
| 44 | 1880 | }(); | |
| 45 | 1880 | return table; | |
| 46 | } | ||
| 47 | |||
| 48 | } // namespace | ||
| 49 | |||
| 50 | ✗ | StarCatalogue SyntheticStars(std::size_t count, std::uint64_t seed, double faintest) { | |
| 51 | // SplitMix64 (Steele, Lea and Flood 2014), for numbers that | ||
| 52 | // repeat on every platform. | ||
| 53 | ✗ | std::uint64_t state = seed; | |
| 54 | ✗ | const auto uniform = [&state] { | |
| 55 | ✗ | std::uint64_t z = (state += 0x9E3779B97F4A7C15ull); | |
| 56 | ✗ | z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ull; | |
| 57 | ✗ | z = (z ^ (z >> 27)) * 0x94D049BB133111EBull; | |
| 58 | ✗ | z ^= z >> 31; | |
| 59 | ✗ | return (static_cast<double>(z >> 11) + 0.5) / 9007199254740992.0; | |
| 60 | ✗ | }; | |
| 61 | // Counts brighter than m grow as 10^(0.48 m), so magnitudes | ||
| 62 | // follow that by inversion between the brightest star's (about | ||
| 63 | // -1.5) and the faintest. | ||
| 64 | ✗ | constexpr double kGrowth = 0.48; | |
| 65 | ✗ | constexpr double kBrightest = -1.5; | |
| 66 | ✗ | const double low = std::pow(10.0, kGrowth * kBrightest); | |
| 67 | ✗ | const double high = std::pow(10.0, kGrowth * faintest); | |
| 68 | ✗ | StarCatalogue stars; | |
| 69 | ✗ | for (std::size_t i = 0; i < count; ++i) { | |
| 70 | ✗ | stars.rightAscension.push_back(static_cast<float>(360.0 * uniform())); | |
| 71 | ✗ | stars.declination.push_back( | |
| 72 | ✗ | static_cast<float>(glm::degrees(std::asin(2.0 * uniform() - 1.0)))); | |
| 73 | ✗ | stars.magnitude.push_back( | |
| 74 | ✗ | static_cast<float>(std::log10(low + (high - low) * uniform()) / kGrowth)); | |
| 75 | // A normal spread of B-V about 0.6 (Box-Muller), within the | ||
| 76 | // stars' range. | ||
| 77 | ✗ | const double normal = | |
| 78 | ✗ | std::sqrt(-2.0 * std::log(uniform())) * std::cos(2.0 * glm::pi<double>() * uniform()); | |
| 79 | ✗ | stars.colorIndex.push_back(static_cast<float>(std::clamp(0.6 + 0.45 * normal, -0.3, 2.0))); | |
| 80 | } | ||
| 81 | ✗ | return stars; | |
| 82 | ✗ | } | |
| 83 | |||
| 84 | 3 | double StarIlluminance(double magnitude, double zeroPoint) { | |
| 85 | 3 | return zeroPoint * std::pow(10.0, -0.4 * magnitude); | |
| 86 | } | ||
| 87 | |||
| 88 | 3 | double ColorIndexToTemperature(double colorIndex) { | |
| 89 | 3 | return 4600.0 * (1.0 / (0.92 * colorIndex + 1.7) + 1.0 / (0.92 * colorIndex + 0.62)); | |
| 90 | } | ||
| 91 | |||
| 92 | 1880 | double BlackbodySpectralIrradiance(double wavelengthNm, double kelvin, double lux) { | |
| 93 | 1880 | kelvin = std::clamp(kelvin, kMinKelvin, kMaxKelvin); | |
| 94 | 1880 | const double t = std::log(kelvin / kMinKelvin) / std::log(kMaxKelvin / kMinKelvin) * | |
| 95 | 1880 | static_cast<double>(kTableSize - 1); | |
| 96 | 1880 | const std::size_t i = std::min(static_cast<std::size_t>(t), kTableSize - 2); | |
| 97 | 1880 | const double f = t - static_cast<double>(i); | |
| 98 | 1880 | const std::array<double, kTableSize>& table = LuxTable(); | |
| 99 | // Interpolated in log space, where the table is smooth. | ||
| 100 | 1880 | const double luxPerUnit = std::exp(std::log(table[i]) * (1.0 - f) + std::log(table[i + 1]) * f); | |
| 101 | 1880 | return lux / luxPerUnit * PlanckRadiance(wavelengthNm, kelvin) / PlanckRadiance(555.0, kelvin); | |
| 102 | } | ||
| 103 | |||
| 104 | 3 | double JulianDate(int year, int month, int day, double hours) { | |
| 105 | 3 | if (month <= 2) { | |
| 106 | 1 | year -= 1; | |
| 107 | 1 | month += 12; | |
| 108 | } | ||
| 109 | 3 | const double a = std::floor(year / 100.0); | |
| 110 | 3 | const double b = 2.0 - a + std::floor(a / 4.0); | |
| 111 | 3 | return std::floor(365.25 * (year + 4716)) + std::floor(30.6001 * (month + 1)) + day + b - | |
| 112 | 3 | 1524.5 + hours / 24.0; | |
| 113 | } | ||
| 114 | |||
| 115 | 3 | double GreenwichMeanSiderealTime(double julianDate) { | |
| 116 | 3 | const double days = julianDate - 2451545.0; | |
| 117 | 3 | const double centuries = days / 36525.0; | |
| 118 | 3 | const double degrees = 280.46061837 + 360.98564736629 * days + | |
| 119 | 3 | 0.000387933 * centuries * centuries - | |
| 120 | 3 | centuries * centuries * centuries / 38710000.0; | |
| 121 | 3 | const double wrapped = std::fmod(degrees, 360.0); | |
| 122 | 3 | return wrapped < 0.0 ? wrapped + 360.0 : wrapped; | |
| 123 | } | ||
| 124 | |||
| 125 | 1 | glm::dmat4 CelestialToLocal(double latitude, double longitude, double julianDate) { | |
| 126 | 1 | const double localSidereal = glm::radians(GreenwichMeanSiderealTime(julianDate) + longitude); | |
| 127 | 1 | const double phi = glm::radians(latitude); | |
| 128 | // To hour-angle axes: x the equator's point on the meridian, y | ||
| 129 | // east, z the pole. | ||
| 130 | 1 | const double c = std::cos(localSidereal); | |
| 131 | 1 | const double s = std::sin(localSidereal); | |
| 132 | 1 | const glm::dmat3 hourAngle(glm::dvec3(c, -s, 0.0), glm::dvec3(s, c, 0.0), | |
| 133 | 1 | glm::dvec3(0.0, 0.0, 1.0)); | |
| 134 | // Those axes in east, north, up: the meridian's equator point due | ||
| 135 | // south at altitude 90 - latitude, and the pole due north at | ||
| 136 | // altitude latitude. | ||
| 137 | 1 | const glm::dmat3 site(glm::dvec3(0.0, -std::sin(phi), std::cos(phi)), glm::dvec3(1.0, 0.0, 0.0), | |
| 138 | 1 | glm::dvec3(0.0, std::cos(phi), std::sin(phi))); | |
| 139 | 1 | return glm::dmat4(site * hourAngle); | |
| 140 | } | ||
| 141 | |||
| 142 | } // namespace ernie | ||
| 143 |