From ed17d0f81d046b7cf8bf8692ef3e187f6940cb49 Mon Sep 17 00:00:00 2001 From: AstroChara Date: Thu, 24 Jul 2025 23:24:58 +0700 Subject: [PATCH] "Update astro constants with actual standards" 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. --- src/celastro/astro.cpp | 88 ++++++++++++++++++++++++++++++++++++++---- src/celastro/astro.h | 68 +++++++++++++++++++++----------- 2 files changed, 126 insertions(+), 30 deletions(-) diff --git a/src/celastro/astro.cpp b/src/celastro/astro.cpp index 62a421b657..5c669fa997 100644 --- a/src/celastro/astro.cpp +++ b/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 // // 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) { diff --git a/src/celastro/astro.h b/src/celastro/astro.h index e504b473ba..fae920d5ab 100644 --- a/src/celastro/astro.h +++ b/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 // // 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 -constexpr inline auto EARTH_RADIUS = detail::enable_if_fp(6378.14L); +constexpr inline auto EARTH_RADIUS = detail::enable_if_fp(6378.1L); // IAU 2015 Resolution B3 template -constexpr inline auto JUPITER_RADIUS = detail::enable_if_fp(71492.0L); +constexpr inline auto JUPITER_RADIUS = detail::enable_if_fp(71492.0L); // IAU 2015 Resolution B3 template -constexpr inline auto SOLAR_RADIUS = detail::enable_if_fp(696000.0L); +constexpr inline auto SOLAR_RADIUS = detail::enable_if_fp(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 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; }