"Update astro constants with up-to-date standards" (#2371)

Split from #1952, as suggested by Askaniy themself, allowing this component of that PR to be updated even if the PSF itself remains work in progress.

Author: AstroChara

src/celastro/astro.cpp

@@ -1,6 +1,6 @@
 // astro.cpp
 //
-// Copyright (C) 2001-2009, the Celestia Development Team
+// Copyright (C) 2001-present, the Celestia Development Team
 // Original version by Chris Laurel <[email protected]>
 //
 // This program is free software; you can redistribute it and/or
@@ -54,20 +54,51 @@ negateIf(double& d, bool condition)
 
 } // end unnamed namespace
 
+
+// Computes the luminosity of a perfectly reflective disc.
+// It is also used as an upper bound for the irradiance of an object when culling invisible objects.
+// The function translates luminosity into luminosity of the same units
+// (distanceFromSun and objRadius must also be in the same units).
+float reflectedLuminosity(float sunLuminosity,
+                          float distanceFromSun,
+                          float objRadius)
+{
+    float lengthRatio = objRadius / distanceFromSun;
+    return sunLuminosity * 0.25 * lengthRatio * lengthRatio;
+}
+
+
+// The following notation rules apply in the functions below and in the code in general:
+// - Luminosity is implied in solar units (in SI, flux is measured in W)
+// - Irradiance is implied in vegan units (in SI, it is measured in W/m^2)
+
+// Absolute magnitude is the logarithmic inverse of luminosity.
+// Apparent magnitude is the logarithmic inverse of irradiance.
+
+
+// Luminosity conversions:
+
 float
 lumToAbsMag(float lum)
 {
     return SOLAR_ABSMAG - std::log(lum) * LN_MAG;
 }
 
-// Return the apparent magnitude of a star with lum times solar
-// luminosity viewed at lyrs light years
 float
 lumToAppMag(float lum, float lyrs)
 {
     return absToAppMag(lumToAbsMag(lum), lyrs);
 }
 
+float
+lumToIrradiance(float lum, float km)
+{
+    return lum * SOLAR_POWER / (math::sphereArea(km * 1000) * VEGAN_IRRADIANCE);
+}
+
+
+// Magnitude conversions:
+
 float
 absMagToLum(float mag)
 {
@@ -80,6 +111,51 @@ appMagToLum(float mag, float lyrs)
     return absMagToLum(appToAbsMag(mag, lyrs));
 }
 
+float
+absMagToIrradiance(float mag, float km)
+{
+    return lumToIrradiance(absMagToLum(mag), km);
+}
+
+
+// Logarithmic magnitude system <-> linear irradiance system in Vega units:
+
+float
+magToIrradiance(float mag)
+{
+    return std::exp(- mag / LN_MAG);
+    // slower solution:
+    // return std::pow(10.0f, -0.4f * mag);
+}
+
+float
+irradianceToMag(float irradiance)
+{
+    return - std::log(irradiance) * LN_MAG;
+    // equivalent solution:
+    // return -2.5f * std::log10(irradiance);
+}
+
+
+// Faintest star magnitude system <-> exposure time:
+
+float
+faintestMagToExposure(float faintestMag)
+{
+    return std::exp(faintestMag / LN_MAG) * LOWEST_IRRADIATION;
+    // slower solution:
+    // return std::pow(10.0f, 0.4f * faintestMag) * LOWEST_IRRADIATION;
+}
+
+float
+exposureToFaintestMag(float exposure)
+{
+    return std::log(exposure / LOWEST_IRRADIATION) * LN_MAG;
+    // equivalent solution:
+    // return 2.5f * std::log10(exposure / LOWEST_IRRADIATION);
+}
+
+
 void
 decimalToDegMinSec(double angle, int& degrees, int& minutes, double& seconds)
 {
@@ -108,8 +184,7 @@ decimalToHourMinSec(double angle, int& hours, int& minutes, double& seconds)
     seconds = (B - (double) minutes) * 60.0;
 }
 
-// Convert equatorial coordinates to Cartesian celestial (or ecliptical)
-// coordinates.
+// Convert equatorial coordinates to Cartesian celestial (or ecliptical) coordinates.
 Eigen::Vector3f
 equatorialToCelestialCart(float ra, float dec, float distance)
 {
@@ -129,8 +204,7 @@ equatorialToCelestialCart(float ra, float dec, float distance)
     return EQUATORIAL_TO_ECLIPTIC_MATRIX_F * Eigen::Vector3f(x, y, z);
 }
 
-// Convert equatorial coordinates to Cartesian celestial (or ecliptical)
-// coordinates.
+// Convert equatorial coordinates to Cartesian celestial (or ecliptical) coordinates.
 Eigen::Vector3d
 equatorialToCelestialCart(double ra, double dec, double distance)
 {

src/celastro/astro.h

@@ -1,6 +1,6 @@
 // astro.h
 //
-// Copyright (C) 2001-2009, the Celestia Development Team
+// Copyright (C) 2001-present, the Celestia Development Team
 // Original version by Chris Laurel <[email protected]>
 //
 // This program is free software; you can redistribute it and/or
@@ -25,8 +25,36 @@
 namespace celestia::astro
 {
 
-constexpr inline float SOLAR_ABSMAG = 4.83f;
-constexpr inline float LN_MAG = 1.0857362f; // 5/ln(100)
+// Angle between J2000 mean equator and the ecliptic plane: 23 deg 26' 21".448
+// Seidelmann, Explanatory Supplement to the Astronomical Almanac (1992), eqn 3.222-1
+constexpr inline double J2000Obliquity    = 23.4392911_deg;
+
+// CODATA 2022
+constexpr inline double speedOfLight      = 299792.458; // km/s
+constexpr inline double G                 = 6.67430e-11; // N m^2 / kg^2
+
+// IAU 2015 Resolution B3 + CODATA 2022
+constexpr inline double SolarMass         = 1.3271244e20 / G; // kg
+constexpr inline double EarthMass         = 3.986004e14 / G; // kg
+constexpr inline double LunarMass         = 7.346e22; // kg
+constexpr inline double JupiterMass       = 1.2668653e17 / G; // kg
+
+// https://mips.as.arizona.edu/~cnaw/sun.html for Johnson V filter
+constexpr inline float SOLAR_ABSMAG       = 4.81f;
+
+// IAU 2015 Resolution B3
+constexpr inline float SOLAR_IRRADIANCE   = 1361.0f; // W / m^2
+constexpr inline float SOLAR_POWER        = 3.828e26f; // W
+
+// Bessel (1979) for Johnson V filter
+constexpr inline float VEGAN_IRRADIANCE   = 3.640e-11f; // W / m^2
+
+// Auxiliary magnitude conversion factor
+constexpr inline float LN_MAG             = 1.0857362f; // 5/ln(100)
+
+// Lowest screen brightness of a point to render
+constexpr inline float LOWEST_IRRADIATION = 1.0f / 255.0f;
+//                                        = 1.0f / (255.0f * 12.92f); after implementing gamma correction
 
 namespace detail
 {
@@ -60,19 +88,27 @@ constexpr inline double SECONDS_PER_DEG = 3600.0;
 constexpr inline double DEG_PER_HRA     = 15.0;
 
 template<typename T>
-constexpr inline auto EARTH_RADIUS = detail::enable_if_fp<T>(6378.14L);
+constexpr inline auto EARTH_RADIUS = detail::enable_if_fp<T>(6378.1L); // IAU 2015 Resolution B3
 
 template<typename T>
-constexpr inline auto JUPITER_RADIUS = detail::enable_if_fp<T>(71492.0L);
+constexpr inline auto JUPITER_RADIUS = detail::enable_if_fp<T>(71492.0L); // IAU 2015 Resolution B3
 
 template<typename T>
-constexpr inline auto SOLAR_RADIUS = detail::enable_if_fp<T>(696000.0L);
+constexpr inline auto SOLAR_RADIUS = detail::enable_if_fp<T>(695700.0L); // IAU 2015 Resolution B3
+
+float reflectedLuminosity(float sunLuminosity, float distanceFromSun, float objRadius);
 
 // Magnitude conversions
 float lumToAbsMag(float lum);
 float lumToAppMag(float lum, float lyrs);
+float lumToIrradiance(float lum, float km);
 float absMagToLum(float mag);
 float appMagToLum(float mag, float lyrs);
+float absMagToIrradiance(float mag, float lyrs);
+float magToIrradiance(float mag);
+float irradianceToMag(float irradiance);
+float faintestMagToExposure(float faintestMag);
+float exposureToFaintestMag(float exposure);
 
 template<class T>
 CELESTIA_CMATH_CONSTEXPR T
@@ -175,33 +211,19 @@ void anomaly(double meanAnomaly, double eccentricity,
              double& trueAnomaly, double& eccentricAnomaly);
 double meanEclipticObliquity(double jd);
 
-constexpr inline double speedOfLight     = 299792.458; // km/s
-constexpr inline double G                = 6.672e-11; // N m^2 / kg^2; gravitational constant
-constexpr inline double SolarMass        = 1.989e30;
-constexpr inline double EarthMass        = 5.972e24;
-constexpr inline double LunarMass        = 7.346e22;
-constexpr inline double JupiterMass      = 1.898e27;
-
-// Angle between J2000 mean equator and the ecliptic plane.
-// 23 deg 26' 21".448 (Seidelmann, _Explanatory Supplement to the
-// Astronomical Almanac_ (1992), eqn 3.222-1.
-constexpr inline double J2000Obliquity   = 23.4392911_deg;
-
-constexpr inline double SOLAR_IRRADIANCE = 1367.6; // Watts / m^2
-constexpr inline double SOLAR_POWER      = 3.8462e26;  // in Watts
 
 namespace literals
 {
 
-constexpr long double operator ""_au(long double au)
+constexpr long double operator "" _au (long double au)
 {
     return AUtoKilometers(au);
 }
-constexpr long double operator ""_ly(long double ly)
+constexpr long double operator "" _ly (long double ly)
 {
     return lightYearsToKilometers(ly);
 }
-constexpr long double operator ""_c(long double n)
+constexpr long double operator "" _c (long double n)
 {
     return speedOfLight * n;
 }