From 33116aad4a7ee264438a84f8afbac52aac1df1a0 Mon Sep 17 00:00:00 2001 From: Andrew Tribick Date: Fri, 31 Dec 2021 14:04:47 +0100 Subject: [PATCH] Backport C++20 numbers header --- src/celcompat/CMakeLists.txt | 4 +- src/celcompat/charconv.h | 2 +- src/celcompat/{cc.cpp => charconv_impl.cpp} | 5 +- src/celcompat/{cc.h => charconv_impl.h} | 2 +- src/celcompat/numbers.h | 37 ++++++++++++ src/celcompat/numbers_impl.h | 41 ++++++++++++++ src/celengine/astro.cpp | 21 ++++--- src/celengine/axisarrow.cpp | 5 +- src/celengine/body.cpp | 12 ++-- src/celengine/glmarker.cpp | 10 ++-- src/celengine/lodspheremesh.cpp | 5 +- src/celengine/observer.cpp | 11 ++-- src/celengine/observer.h | 3 +- src/celengine/parseobject.cpp | 7 ++- src/celengine/planetgrid.cpp | 5 +- src/celengine/render.cpp | 8 +-- src/celengine/renderglsl.cpp | 6 +- src/celengine/skygrid.cpp | 48 ++++++++-------- src/celengine/spheremesh.cpp | 9 +-- src/celengine/universe.cpp | 3 +- src/celengine/visibleregion.cpp | 3 +- src/celephem/customorbit.cpp | 63 +++++++++++---------- src/celephem/customrotation.cpp | 11 ++-- src/celephem/orbit.cpp | 9 +-- src/celephem/precession.cpp | 5 +- src/celephem/rotation.cpp | 9 +-- src/celephem/samporient.cpp | 3 +- src/celephem/spicerotation.cpp | 5 +- src/celephem/vsop87.cpp | 5 +- src/celestia/celestiacore.cpp | 9 +-- src/celestia/gtk/dialog-eclipse.cpp | 3 +- src/celestia/qt/qtinfopanel.cpp | 6 +- src/celestia/win32/wineclipses.cpp | 3 +- src/celmath/mathlib.h | 15 +++-- src/celmath/randutils.h | 3 +- src/tools/atmosphere/scattersim.cpp | 46 ++++++++------- src/tools/cmod/cmodsphere/cmodsphere.cpp | 8 +-- src/tools/stardb/startextdump.cpp | 5 +- 38 files changed, 288 insertions(+), 167 deletions(-) rename src/celcompat/{cc.cpp => charconv_impl.cpp} (99%) rename src/celcompat/{cc.h => charconv_impl.h} (99%) create mode 100644 src/celcompat/numbers.h create mode 100644 src/celcompat/numbers_impl.h diff --git a/src/celcompat/CMakeLists.txt b/src/celcompat/CMakeLists.txt index 7611596b7..d159c04a7 100644 --- a/src/celcompat/CMakeLists.txt +++ b/src/celcompat/CMakeLists.txt @@ -1,7 +1,7 @@ if (NOT HAVE_FLOAT_CHARCONV) set(CELCOMPAT_SOURCES - cc.cpp - cc.h + charconv_impl.cpp + charconv_impl.h ) endif() diff --git a/src/celcompat/charconv.h b/src/celcompat/charconv.h index 27cd25471..552faed0a 100644 --- a/src/celcompat/charconv.h +++ b/src/celcompat/charconv.h @@ -11,5 +11,5 @@ namespace celestia::compat using std::from_chars; } #else -#include "cc.h" +#include "charconv_impl.h" #endif diff --git a/src/celcompat/cc.cpp b/src/celcompat/charconv_impl.cpp similarity index 99% rename from src/celcompat/cc.cpp rename to src/celcompat/charconv_impl.cpp index a4b0fd0d4..8b383a999 100644 --- a/src/celcompat/cc.cpp +++ b/src/celcompat/charconv_impl.cpp @@ -1,4 +1,4 @@ -// cc.cpp +// charconv_impl.cpp // // Copyright (C) 2021-present, Celestia Development Team. // @@ -17,7 +17,8 @@ #include #include -#include "cc.h" + +#include "charconv_impl.h" namespace celestia::compat { diff --git a/src/celcompat/cc.h b/src/celcompat/charconv_impl.h similarity index 99% rename from src/celcompat/cc.h rename to src/celcompat/charconv_impl.h index d1aaf3d1c..9725f857b 100644 --- a/src/celcompat/cc.h +++ b/src/celcompat/charconv_impl.h @@ -1,4 +1,4 @@ -// cc.h +// charconv_impl.h // // Copyright (C) 2021-present, Celestia Development Team. // diff --git a/src/celcompat/numbers.h b/src/celcompat/numbers.h new file mode 100644 index 000000000..df05d87a7 --- /dev/null +++ b/src/celcompat/numbers.h @@ -0,0 +1,37 @@ +#pragma once + +#if __cpp_lib_math_constants >= 201907L +#include +namespace celestia::numbers +{ +using std::numbers::e_v; +using std::numbers::log2e_v; +using std::numbers::log10e_v; +using std::numbers::pi_v; +using std::numbers::inv_pi_v; +using std::numbers::inv_sqrtpi_v; +using std::numbers::ln2_v; +using std::numbers::ln10_v; +using std::numbers::sqrt2_v; +using std::numbers::sqrt3_v; +using std::numbers::inv_sqrt3_v; +using std::numbers::egamma_v; +using std::numbers::phi_v; + +using std::numbers::e; +using std::numbers::log2e; +using std::numbers::log10e; +using std::numbers::pi; +using std::numbers::inv_pi; +using std::numbers::inv_sqrtpi; +using std::numbers::ln2; +using std::numbers::ln10; +using std::numbers::sqrt2; +using std::numbers::sqrt3; +using std::numbers::inv_sqrt3; +using std::numbers::egamma; +using std::numbers::phi; +} +#else +#include "numbers_impl.h" +#endif diff --git a/src/celcompat/numbers_impl.h b/src/celcompat/numbers_impl.h new file mode 100644 index 000000000..a55fee060 --- /dev/null +++ b/src/celcompat/numbers_impl.h @@ -0,0 +1,41 @@ +#pragma once + +#include + +namespace celestia::numbers +{ +namespace detail +{ +template +using enable_if_fp = std::enable_if_t, T>; +} + +// constants are given to 128 bits of precision, computed with sollya +template constexpr inline T e_v = detail::enable_if_fp{0x1.5bf0a8b1457695355fb8ac404e7a79e4p1}; +template constexpr inline T log2e_v = detail::enable_if_fp{0x1.71547652b82fe1777d0ffda0d23a7d12p0}; +template constexpr inline T log10e_v = detail::enable_if_fp{0x1.bcb7b1526e50e32a6ab7555f5a67b864p-2}; +template constexpr inline T pi_v = detail::enable_if_fp{0x1.921fb54442d18469898cc51701b839a2p1}; +template constexpr inline T inv_pi_v = detail::enable_if_fp{0x1.45f306dc9c882a53f84eafa3ea69bb82p-2}; +template constexpr inline T inv_sqrtpi_v = detail::enable_if_fp{0x1.20dd750429b6d11ae3a914fed7fd8688p-1}; +template constexpr inline T ln2_v = detail::enable_if_fp{0x1.62e42fefa39ef35793c7673007e5ed5ep-1}; +template constexpr inline T ln10_v = detail::enable_if_fp{0x1.26bb1bbb5551582dd4adac5705a61452p1}; +template constexpr inline T sqrt2_v = detail::enable_if_fp{0x1.6a09e667f3bcc908b2fb1366ea957d3ep0}; +template constexpr inline T sqrt3_v = detail::enable_if_fp{0x1.bb67ae8584caa73b25742d7078b83b8ap0}; +template constexpr inline T inv_sqrt3_v = detail::enable_if_fp{0x1.279a74590331c4d218f81e4afb257d06p-1}; +template constexpr inline T egamma_v = detail::enable_if_fp{0x1.2788cfc6fb618f49a37c7f0202a596aep-1}; +template constexpr inline T phi_v = detail::enable_if_fp{0x1.9e3779b97f4a7c15f39cc0605cedc834p0}; + +constexpr inline double e = e_v; +constexpr inline double log2e = log2e_v; +constexpr inline double log10e = log10e_v; +constexpr inline double pi = pi_v; +constexpr inline double inv_pi = inv_pi_v; +constexpr inline double inv_sqrtpi = inv_sqrtpi_v; +constexpr inline double ln2 = ln2_v; +constexpr inline double ln10 = ln10_v; +constexpr inline double sqrt2 = sqrt2_v; +constexpr inline double sqrt3 = sqrt3_v; +constexpr inline double inv_sqrt3 = inv_sqrt3_v; +constexpr inline double egamma = egamma_v; +constexpr inline double phi = phi_v; +} diff --git a/src/celengine/astro.cpp b/src/celengine/astro.cpp index 19ab65003..eadefa422 100644 --- a/src/celengine/astro.cpp +++ b/src/celengine/astro.cpp @@ -14,6 +14,7 @@ #include #include "astro.h" #include "univcoord.h" +#include #include #include @@ -140,7 +141,7 @@ static const UnitDefinition angleUnits[] = { "arcmin", 1.0 / MINUTES_PER_DEG }, { "deg", 1.0 }, { "hRA", DEG_PER_HRA }, - { "rad", 180.0 / PI }, + { "rad", 180.0 / celestia::numbers::pi }, }; @@ -208,8 +209,9 @@ void astro::decimalToHourMinSec(double angle, int& hours, int& minutes, double& Eigen::Vector3f astro::equatorialToCelestialCart(float ra, float dec, float distance) { - double theta = ra / 24.0 * PI * 2 + PI; - double phi = (dec / 90.0 - 1.0) * PI / 2; + using celestia::numbers::pi; + double theta = ra / 24.0 * pi * 2 + pi; + double phi = (dec / 90.0 - 1.0) * pi / 2; double x = cos(theta) * sin(phi) * distance; double y = cos(phi) * distance; double z = -sin(theta) * sin(phi) * distance; @@ -223,8 +225,9 @@ astro::equatorialToCelestialCart(float ra, float dec, float distance) Eigen::Vector3d astro::equatorialToCelestialCart(double ra, double dec, double distance) { - double theta = ra / 24.0 * PI * 2 + PI; - double phi = (dec / 90.0 - 1.0) * PI / 2; + using celestia::numbers::pi; + double theta = ra / 24.0 * pi * 2 + pi; + double phi = (dec / 90.0 - 1.0) * pi / 2; double x = cos(theta) * sin(phi) * distance; double y = cos(phi) * distance; double z = -sin(theta) * sin(phi) * distance; @@ -239,8 +242,9 @@ astro::equatorialToCelestialCart(double ra, double dec, double distance) Eigen::Vector3f astro::equatorialToEclipticCartesian(float ra, float dec, float distance) { - double theta = ra / 24.0 * PI * 2 + PI; - double phi = (dec / 90.0 - 1.0) * PI / 2; + using celestia::numbers::pi; + double theta = ra / 24.0 * pi * 2 + pi; + double phi = (dec / 90.0 - 1.0) * pi / 2; double x = cos(theta) * sin(phi) * distance; double y = cos(phi) * distance; double z = -sin(theta) * sin(phi) * distance; @@ -252,11 +256,12 @@ astro::equatorialToEclipticCartesian(float ra, float dec, float distance) void astro::anomaly(double meanAnomaly, double eccentricity, double& trueAnomaly, double& eccentricAnomaly) { + using celestia::numbers::pi; double e, delta, err; double tol = 0.00000001745; int iterations = 20; // limit while() to maximum of 20 iterations. - e = meanAnomaly - 2*PI * (int) (meanAnomaly / (2*PI)); + e = meanAnomaly - 2*pi * (int) (meanAnomaly / (2*pi)); err = 1; while(abs(err) > tol && iterations > 0) { diff --git a/src/celengine/axisarrow.cpp b/src/celengine/axisarrow.cpp index eb803bced..c3715bbe3 100644 --- a/src/celengine/axisarrow.cpp +++ b/src/celengine/axisarrow.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include "axisarrow.h" #include "body.h" @@ -58,7 +59,7 @@ static size_t initArrow(VertexObject &vo) for (unsigned i = 0; i <= nSections; i++) { float c, s; - sincos((i * 2.0f * (float)PI) / nSections, c, s); + sincos((i * 2.0f * celestia::numbers::pi_v) / nSections, c, s); // circle at bottom Vector3f v0(shaftRadius * c, shaftRadius * s, 0.0f); @@ -581,7 +582,7 @@ BodyAxisArrows::BodyAxisArrows(const Body& _body) : Quaterniond BodyAxisArrows::getOrientation(double tdb) const { - return (Quaterniond(AngleAxis(PI, Vector3d::UnitY())) * body.getEclipticToBodyFixed(tdb)).conjugate(); + return (Quaterniond(AngleAxis(celestia::numbers::pi, Vector3d::UnitY())) * body.getEclipticToBodyFixed(tdb)).conjugate(); } diff --git a/src/celengine/body.cpp b/src/celengine/body.cpp index b13a34ef2..5646a3f8e 100644 --- a/src/celengine/body.cpp +++ b/src/celengine/body.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -261,7 +262,7 @@ float Body::getBoundingRadius() const if (geometry == InvalidResource) return radius; - return radius * 1.7320508f; // sqrt(3) + return radius * celestia::numbers::sqrt3_v; // sqrt(3) } @@ -300,7 +301,7 @@ float Body::getDensity() const // @astro::EarthMass unit is kg // @radius unit km // so we divide density by 1e9 to have kg/m^3 - double volume = 4.0 / 3.0 * PI * ::pow(radius, 3); + double volume = 4.0 / 3.0 * celestia::numbers::pi * ::pow(radius, 3); return (float) mass * astro::EarthMass / 1e9 / volume; } @@ -725,8 +726,9 @@ Matrix4d Body::getBodyFixedToAstrocentric(double tdb) const Vector3d Body::planetocentricToCartesian(double lon, double lat, double alt) const { - double phi = -degToRad(lat) + PI / 2; - double theta = degToRad(lon) - PI; + using celestia::numbers::pi; + double phi = -degToRad(lat) + pi / 2; + double theta = degToRad(lon) - pi; Vector3d pos(cos(theta) * sin(phi), cos(phi), @@ -749,7 +751,7 @@ Vector3d Body::cartesianToPlanetocentric(const Vector3d& v) const { Vector3d w = v.normalized(); - double lat = PI / 2.0 - acos(w.y()); + double lat = celestia::numbers::pi / 2.0 - acos(w.y()); double lon = atan2(w.z(), -w.x()); return Vector3d(lon, lat, v.norm() - getRadius()); diff --git a/src/celengine/glmarker.cpp b/src/celengine/glmarker.cpp index 3325ade13..ceaea41e7 100644 --- a/src/celengine/glmarker.cpp +++ b/src/celengine/glmarker.cpp @@ -9,6 +9,7 @@ // of the License, or (at your option) any later version. #include +#include #include #include #include "marker.h" @@ -148,7 +149,7 @@ static void fillCircleValue(GLfloat* data, int size, float scale) float s, c; for (int i = 0; i < size; i++) { - sincos((float) (2 * i) / (float) size * ((float) PI), s, c); + sincos((float) (2 * i) / (float) size * celestia::numbers::pi_v, s, c); data[i * 2] = c * scale; data[i * 2 + 1] = s * scale; } @@ -320,7 +321,7 @@ static void initEclipticVO(VertexObject& vo) float scale = 1000.0f; for (int i = 0; i < eclipticCount; i++) { - sincos((float) (2 * i) / (float) eclipticCount * ((float) PI), s, c); + sincos((float) (2 * i) / (float) eclipticCount * celestia::numbers::pi_v, s, c); ecliptic[i * 3] = c * scale; ecliptic[i * 3 + 1] = 0; ecliptic[i * 3 + 2] = s * scale; @@ -604,7 +605,8 @@ void Renderer::renderCrosshair(float selectionSizeInPixels, const float cursorPulsePeriod = 1.5f; float cursorRadius = selectionSizeInPixels + cursorMinRadius; - cursorRadius += cursorRadiusVariability * (float) (0.5 + 0.5 * sin(tsec * 2 * PI / cursorPulsePeriod)); + cursorRadius += cursorRadiusVariability * + (float) (0.5 + 0.5 * sin(tsec * 2 * celestia::numbers::pi / cursorPulsePeriod)); // Enlarge the size of the cross hair sligtly when the selection // has a large apparent size @@ -621,7 +623,7 @@ void Renderer::renderCrosshair(float selectionSizeInPixels, const unsigned int markCount = 4; for (unsigned int i = 0; i < markCount; i++) { - float theta = (float) (PI / 4.0) + (float) i / (float) markCount * (float) (2 * PI); + float theta = (float) (celestia::numbers::pi / 4.0) + (float) i / (float) markCount * (float) (2 * celestia::numbers::pi); prog->floatParam("angle") = theta; markerVO.draw(GL_TRIANGLES, CrosshairCount, CrosshairOffset); } diff --git a/src/celengine/lodspheremesh.cpp b/src/celengine/lodspheremesh.cpp index 274c69da3..856524b72 100644 --- a/src/celengine/lodspheremesh.cpp +++ b/src/celengine/lodspheremesh.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include "glsupport.h" #include "lodspheremesh.h" @@ -64,14 +65,14 @@ static void InitTrigArrays() int i; for (i = 0; i <= thetaDivisions; i++) { - double theta = (double) i / (double) thetaDivisions * 2.0 * PI; + double theta = (double) i / (double) thetaDivisions * 2.0 * celestia::numbers::pi; sinTheta[i] = (float) sin(theta); cosTheta[i] = (float) cos(theta); } for (i = 0; i <= phiDivisions; i++) { - double phi = ((double) i / (double) phiDivisions - 0.5) * PI; + double phi = ((double) i / (double) phiDivisions - 0.5) * celestia::numbers::pi; sinPhi[i] = (float) sin(phi); cosPhi[i] = (float) cos(phi); } diff --git a/src/celengine/observer.cpp b/src/celengine/observer.cpp index ce56fb685..a8c7ad9e2 100644 --- a/src/celengine/observer.cpp +++ b/src/celengine/observer.cpp @@ -386,7 +386,8 @@ void Observer::update(double dt, double timeScale) else { v = pow(sin((t - journey.startInterpolation) / - (journey.endInterpolation - journey.startInterpolation) * PI / 2), 2); + (journey.endInterpolation - journey.startInterpolation) * + celestia::numbers::pi / 2.0), 2); } q = journey.initialOrientation.slerp(v, journey.finalOrientation); @@ -493,7 +494,7 @@ void Observer::setLocationFilter(uint64_t _locationFilter) void Observer::reverseOrientation() { - setOrientation(getOrientation() * Quaterniond(AngleAxisd(PI, Vector3d::UnitY()))); + setOrientation(getOrientation() * Quaterniond(AngleAxisd(celestia::numbers::pi, Vector3d::UnitY()))); reverseFlag = !reverseFlag; } @@ -1178,7 +1179,7 @@ void Observer::gotoSelectionLongLat(const Selection& selection, { if (!selection.empty()) { - double phi = -latitude + PI / 2; + double phi = -latitude + celestia::numbers::pi / 2.0; double theta = longitude; double x = cos(theta) * sin(phi); double y = cos(phi); @@ -1240,7 +1241,7 @@ void Observer::getSelectionLongLat(const Selection& selection, distance = bfPos.norm(); longitude = radToDeg(atan2(y, x)); - latitude = radToDeg(PI/2 - acos(z / distance)); + latitude = radToDeg(celestia::numbers::pi/2.0 - acos(z / distance)); } } @@ -1358,7 +1359,7 @@ Vector3f Observer::getPickRay(float x, float y) const Vector3f Observer::getPickRayFisheye(float x, float y) const { float r = hypot(x, y); - float phi = float(PI) * r; + float phi = celestia::numbers::pi_v * r; float sin_phi = sin(phi); float theta = atan2(y, x); float newX = sin_phi * cos(theta); diff --git a/src/celengine/observer.h b/src/celengine/observer.h index 37d1d9cb5..cdd8755f5 100644 --- a/src/celengine/observer.h +++ b/src/celengine/observer.h @@ -18,6 +18,7 @@ #ifndef _CELENGINE_OBSERVER_H_ #define _CELENGINE_OBSERVER_H_ +#include #include #include #include @@ -312,7 +313,7 @@ public: Eigen::Quaterniond trackingOrientation{ Eigen::Quaternionf::Identity() }; // orientation prior to selecting tracking - float fov{ (float) (PI / 4.0) }; + float fov{ static_cast(celestia::numbers::pi / 4.0) }; bool reverseFlag{ false }; uint64_t locationFilter{ ~0ull }; diff --git a/src/celengine/parseobject.cpp b/src/celengine/parseobject.cpp index 9a4874c45..ce7cb31b7 100644 --- a/src/celengine/parseobject.cpp +++ b/src/celengine/parseobject.cpp @@ -16,6 +16,7 @@ #include "trajmanager.h" #include "rotationmanager.h" #include "universe.h" +#include #include #include #ifdef USE_SPICE @@ -833,7 +834,7 @@ CreateFixedRotationModel(double offset, double inclination, double ascendingNode) { - Quaterniond q = YRotation(-PI - offset) * + Quaterniond q = YRotation(-celestia::numbers::pi - offset) * XRotation(-inclination) * YRotation(-ascendingNode); @@ -910,7 +911,7 @@ CreateFixedRotationModel(Hash* rotationData) ascendingNode = degToRad(ascendingNode); } - Quaterniond q = YRotation(-PI - offset) * + Quaterniond q = YRotation(-celestia::numbers::pi - offset) * XRotation(-inclination) * YRotation(-ascendingNode); @@ -939,7 +940,7 @@ CreateFixedAttitudeRotationModel(Hash* rotationData) roll = degToRad(roll); } - Quaterniond q = YRotation(-PI - heading) * + Quaterniond q = YRotation(-celestia::numbers::pi - heading) * XRotation(-tilt) * ZRotation(-roll); diff --git a/src/celengine/planetgrid.cpp b/src/celengine/planetgrid.cpp index 97f58622f..1f593590f 100644 --- a/src/celengine/planetgrid.cpp +++ b/src/celengine/planetgrid.cpp @@ -12,6 +12,7 @@ #include #include +#include #include #include "body.h" #include "planetgrid.h" @@ -105,7 +106,7 @@ PlanetographicGrid::render(Renderer* renderer, return; // Compatibility - Quaterniond q = Quaterniond(AngleAxis(PI, Vector3d::UnitY())) * body.getEclipticToBodyFixed(tdb); + Quaterniond q = Quaterniond(AngleAxis(celestia::numbers::pi, Vector3d::UnitY())) * body.getEclipticToBodyFixed(tdb); Quaternionf qf = q.cast(); // The grid can't be rendered exactly on the planet sphere, or @@ -356,7 +357,7 @@ PlanetographicGrid::InitializeGeometry() xzCircle.reserve((circleSubdivisions + 2) * 2); for (unsigned int i = 0; i <= circleSubdivisions + 1; i++) { - float theta = (float) (2.0 * PI) * (float) i / (float) circleSubdivisions; + float theta = (float) (2.0 * celestia::numbers::pi) * (float) i / (float) circleSubdivisions; float s, c; sincos(theta, s, c); Vector3f thisPointXY(c, s, 0.0f); diff --git a/src/celengine/render.cpp b/src/celengine/render.cpp index 7a479cfe7..1eaa89650 100644 --- a/src/celengine/render.cpp +++ b/src/celengine/render.cpp @@ -341,7 +341,7 @@ static void BuildGaussianDiscMipLevel(unsigned char* mipPixels, unsigned int size = 1 << log2size; float sigma = fwhm / 2.3548f; float isig2 = 1.0f / (2.0f * sigma * sigma); - float s = 1.0f / (sigma * (float) sqrt(2.0 * PI)); + float s = 1.0f / (sigma * (float) sqrt(2.0 * celestia::numbers::pi)); for (unsigned int i = 0; i < size; i++) { @@ -1999,7 +1999,7 @@ void Renderer::renderEllipsoidAtmosphere(const Atmosphere& atmosphere, { // We want rays with an origin at the eye point and tangent to the the // ellipsoid. - float theta = (float) i / (float) nSlices * 2 * (float) PI; + float theta = (float) i / (float) nSlices * 2 * celestia::numbers::pi_v; Vector3f w = (float) cos(theta) * uAxis + (float) sin(theta) * vAxis; w = w * (float) centerDist; @@ -2705,7 +2705,7 @@ void Renderer::renderObject(const Vector3f& pos, cloudNormalMap = atmosphere->cloudNormalMap.find(textureResolution); } if (atmosphere->cloudSpeed != 0.0f) - cloudTexOffset = (float) (-pfmod(now * atmosphere->cloudSpeed / (2 * PI), 1.0)); + cloudTexOffset = (float) (-pfmod(now * atmosphere->cloudSpeed / (2 * celestia::numbers::pi), 1.0)); } if (obj.geometry == InvalidResource) @@ -3546,7 +3546,7 @@ void Renderer::renderCometTail(const Body& body, float radius = (float) i / (float) nTailPoints * dustTailRadius; for (int j = 0; j < nTailSlices; j++) { - float theta = (float) (2 * PI * (float) j / nTailSlices); + float theta = (float) (2 * celestia::numbers::pi * (float) j / nTailSlices); float s, c; sincos(theta, s, c); CometTailVertex& vtx = cometTailVertices[i * nTailSlices + j]; diff --git a/src/celengine/renderglsl.cpp b/src/celengine/renderglsl.cpp index aaf53f37d..46649ba64 100644 --- a/src/celengine/renderglsl.cpp +++ b/src/celengine/renderglsl.cpp @@ -16,6 +16,8 @@ #include #include +#include + #include #include #include @@ -722,7 +724,7 @@ static void renderRingSystem(GLuint *vboId, GLshort tex[2]; }; - constexpr const float angle = 2.0f * static_cast(PI); + constexpr const float angle = 2.0f * celestia::numbers::pi_v; if (*vboId == 0) { @@ -920,7 +922,7 @@ void renderRings_GLSL(RingSystem& rings, std::size_t i = 0; for (i = 0; i < data->vboId.size() - 1; i++) { - float s = segmentSizeInPixels * tan(PI / nSections); + float s = segmentSizeInPixels * tan(celestia::numbers::pi / nSections); if (s < 30.0f) // TODO: make configurable break; nSections <<= 1; diff --git a/src/celengine/skygrid.cpp b/src/celengine/skygrid.cpp index 4a6000a7b..c4fe05db0 100644 --- a/src/celengine/skygrid.cpp +++ b/src/celengine/skygrid.cpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include "render.h" #include "vecgl.h" @@ -147,9 +148,10 @@ toStandardCoords(const Vector3d& v) template static T angleDiff(T a, T b) { + using celestia::numbers::pi_v; T diff = std::fabs(a - b); - if (diff > PI) - return (T) (2.0 * PI - diff); + if (diff > pi_v) + return (T) (2.0 * pi_v - diff); else return diff; } @@ -346,7 +348,7 @@ SkyGrid::parallelSpacing(double idealSpacing) const unsigned int tableSize = sizeof(DEG_MIN_SEC_SPACING) / sizeof(DEG_MIN_SEC_SPACING[0]); for (unsigned int i = 0; i < tableSize; i++) { - if (PI * (double) DEG_MIN_SEC_SPACING[i] / (double) DEG_MIN_SEC_TOTAL < idealSpacing) + if (celestia::numbers::pi * (double) DEG_MIN_SEC_SPACING[i] / (double) DEG_MIN_SEC_TOTAL < idealSpacing) break; spacing = DEG_MIN_SEC_SPACING[i]; } @@ -375,7 +377,7 @@ SkyGrid::meridianSpacing(double idealSpacing) const for (unsigned int i = 0; i < tableSize; i++) { - if (2 * PI * (double) spacingTable[i] / (double) totalUnits < idealSpacing) + if (2 * celestia::numbers::pi * (double) spacingTable[i] / (double) totalUnits < idealSpacing) break; spacing = spacingTable[i]; } @@ -404,7 +406,7 @@ SkyGrid::render(Renderer& renderer, // 90 degree rotation about the x-axis used to transform coordinates // to Celestia's system. - Quaterniond xrot90 = XRotation(-PI / 2.0); + Quaterniond xrot90 = XRotation(-celestia::numbers::pi / 2.0); double vfov = observer.getFOV(); double viewAspectRatio = (double) windowWidth / (double) windowHeight; @@ -462,7 +464,7 @@ SkyGrid::render(Renderer& renderer, updateAngleRange(thetaC1, thetaC3, &maxDiff, &minTheta, &maxTheta); updateAngleRange(thetaC2, thetaC3, &maxDiff, &minTheta, &maxTheta); - if (std::fabs(maxTheta - minTheta) < PI) + if (std::fabs(maxTheta - minTheta) < celestia::numbers::pi) { if (minTheta > maxTheta) std::swap(minTheta, maxTheta); @@ -495,26 +497,26 @@ SkyGrid::render(Renderer& renderer, if (fabs(viewCenter.z()) < 1.0) centerDec = std::asin(viewCenter.z()); else if (viewCenter.z() < 0.0) - centerDec = -PI / 2.0; + centerDec = -celestia::numbers::pi / 2.0; else - centerDec = PI / 2.0; + centerDec = celestia::numbers::pi / 2.0; double minDec = centerDec - halfFov; double maxDec = centerDec + halfFov; - if (maxDec >= PI / 2.0) + if (maxDec >= celestia::numbers::pi / 2.0) { // view cone contains north pole - maxDec = PI / 2.0; - minTheta = -PI; - maxTheta = PI; + maxDec = celestia::numbers::pi / 2.0; + minTheta = -celestia::numbers::pi; + maxTheta = celestia::numbers::pi; } - else if (minDec <= -PI / 2.0) + else if (minDec <= -celestia::numbers::pi / 2.0) { // view cone contains south pole - minDec = -PI / 2.0; - minTheta = -PI; - maxTheta = PI; + minDec = -celestia::numbers::pi / 2.0; + minTheta = -celestia::numbers::pi; + maxTheta = celestia::numbers::pi; } double idealParallelSpacing = 2.0 * halfFov / MAX_VISIBLE_ARCS; @@ -542,10 +544,10 @@ SkyGrid::render(Renderer& renderer, int raIncrement = meridianSpacing(idealMeridianSpacing); int decIncrement = parallelSpacing(idealParallelSpacing); - int startRa = (int) std::ceil (totalLongitudeUnits * (minTheta / (PI * 2.0f)) / (float) raIncrement) * raIncrement; - int endRa = (int) std::floor(totalLongitudeUnits * (maxTheta / (PI * 2.0f)) / (float) raIncrement) * raIncrement; - int startDec = (int) std::ceil (DEG_MIN_SEC_TOTAL * (minDec / PI) / (float) decIncrement) * decIncrement; - int endDec = (int) std::floor(DEG_MIN_SEC_TOTAL * (maxDec / PI) / (float) decIncrement) * decIncrement; + int startRa = (int) std::ceil (totalLongitudeUnits * (minTheta / (celestia::numbers::pi * 2.0)) / (float) raIncrement) * raIncrement; + int endRa = (int) std::floor(totalLongitudeUnits * (maxTheta / (celestia::numbers::pi * 2.0)) / (float) raIncrement) * raIncrement; + int startDec = (int) std::ceil (DEG_MIN_SEC_TOTAL * (minDec / celestia::numbers::pi) / (float) decIncrement) * decIncrement; + int endDec = (int) std::floor(DEG_MIN_SEC_TOTAL * (maxDec / celestia::numbers::pi) / (float) decIncrement) * decIncrement; // Get the orientation at single precision Quaterniond q = xrot90 * m_orientation * xrot90.conjugate(); @@ -590,7 +592,7 @@ SkyGrid::render(Renderer& renderer, for (int dec = startDec; dec <= endDec; dec += decIncrement) { - double phi = PI * (double) dec / (double) DEG_MIN_SEC_TOTAL; + double phi = celestia::numbers::pi * (double) dec / (double) DEG_MIN_SEC_TOTAL; double cosPhi = cos(phi); double sinPhi = sin(phi); @@ -660,7 +662,7 @@ SkyGrid::render(Renderer& renderer, // Render meridians only to the last latitude circle; this looks better // than spokes radiating from the pole. - double maxMeridianAngle = PI / 2.0 * (1.0 - 2.0 * (double) decIncrement / (double) DEG_MIN_SEC_TOTAL); + double maxMeridianAngle = celestia::numbers::pi / 2.0 * (1.0 - 2.0 * (double) decIncrement / (double) DEG_MIN_SEC_TOTAL); minDec = std::max(minDec, -maxMeridianAngle); maxDec = std::min(maxDec, maxMeridianAngle); arcStep = (maxDec - minDec) / (double) ARC_SUBDIVISIONS; @@ -670,7 +672,7 @@ SkyGrid::render(Renderer& renderer, for (int ra = startRa; ra <= endRa; ra += raIncrement) { - double theta = 2.0 * PI * (double) ra / (double) totalLongitudeUnits; + double theta = 2.0 * celestia::numbers::pi * (double) ra / (double) totalLongitudeUnits; double cosTheta = cos(theta); double sinTheta = sin(theta); diff --git a/src/celengine/spheremesh.cpp b/src/celengine/spheremesh.cpp index b93c40db8..2d0eec852 100644 --- a/src/celengine/spheremesh.cpp +++ b/src/celengine/spheremesh.cpp @@ -14,6 +14,7 @@ #include #include +#include #include #include #include @@ -30,8 +31,8 @@ constexpr const int MinSlices = 3; float SphereMeshParameters::value(float u, float v) const { - float theta = u * static_cast(PI) * 2; - float phi = (v - 0.5f) * static_cast(PI); + float theta = u * celestia::numbers::pi_v * 2; + float phi = (v - 0.5f) * celestia::numbers::pi_v; Eigen::Vector3f p = Eigen::Vector3f(std::cos(phi) * std::cos(theta), std::sin(phi), @@ -65,10 +66,10 @@ void SphereMesh::createSphere() for (int i = 0; i < nRings; i++) { - float phi = (static_cast(i) / static_cast(nRings - 1) - 0.5f) * static_cast(PI); + float phi = (static_cast(i) / static_cast(nRings - 1) - 0.5f) * celestia::numbers::pi_v; for (int j = 0; j <= nSlices; j++) { - float theta = static_cast(j) / static_cast(nSlices) * static_cast(PI * 2.0); + float theta = static_cast(j) / static_cast(nSlices) * static_cast(celestia::numbers::pi * 2.0); auto x = std::cos(phi) * std::cos(theta); auto y = std::sin(phi); auto z = std::cos(phi) * std::sin(theta); diff --git a/src/celengine/universe.cpp b/src/celengine/universe.cpp index d12ae848f..1c4ffc500 100644 --- a/src/celengine/universe.cpp +++ b/src/celengine/universe.cpp @@ -18,6 +18,7 @@ #include "universe.h" #include "timelinephase.h" #include "frametree.h" +#include #include #include #include @@ -400,7 +401,7 @@ static bool ExactPlanetPickTraversal(Body* body, void* info) double sinAngle2 = bodyMiss.norm() / 2.0; - if (sinAngle2 < sin(PI/4.0) && distance > 0.0 && + if (sinAngle2 < sin(celestia::numbers::pi/4.0) && distance > 0.0 && distance <= pickInfo->closestDistance) { pickInfo->closestDistance = distance; diff --git a/src/celengine/visibleregion.cpp b/src/celengine/visibleregion.cpp index 456fc71e4..ed87198fd 100644 --- a/src/celengine/visibleregion.cpp +++ b/src/celengine/visibleregion.cpp @@ -12,6 +12,7 @@ #include #include +#include #include #include "body.h" #include "render.h" @@ -219,7 +220,7 @@ VisibleRegion::render(Renderer* renderer, for (unsigned i = 0; i <= nSections + 1; i++) { - double theta = (double) i / (double) (nSections) * 2.0 * PI; + double theta = (double) i / (double) (nSections) * 2.0 * celestia::numbers::pi; Vector3d w = cos(theta) * uAxis + sin(theta) * vAxis; Vector3d toCenter = ellipsoidTangent(recipSemiAxes, w, e, e_, ee); diff --git a/src/celephem/customorbit.cpp b/src/celephem/customorbit.cpp index 0371e96a3..e0fda836c 100644 --- a/src/celephem/customorbit.cpp +++ b/src/celephem/customorbit.cpp @@ -10,6 +10,7 @@ #include "customorbit.h" #include "vsop87.h" #include "jpleph.h" +#include #include #include #include @@ -236,7 +237,7 @@ static void EclipticToEquatorial(double fEclLat, double fEclLon, dec = asin((sy*ceps)+(cy*seps*sx)); RA = atan(((sx*ceps)-(ty*seps))/cx); if (cx<0) - RA += PI; // account for atan quad ambiguity + RA += celestia::numbers::pi; // account for atan quad ambiguity RA = pfmod(RA, TWOPI); } @@ -345,7 +346,7 @@ void computePlanetCoords(int p, double map, double da, double dhl, double dl, spsi = sin(eclLat); eclLong = atan(y/clo) + om + degToRad(dl); if (clo < 0) - eclLong += PI; + eclLong += celestia::numbers::pi; eclLong = pfmod(eclLong, TWOPI); distance *= KM_PER_AU; } @@ -430,8 +431,8 @@ class MercuryOrbit : public CachingOrbit eclLong, eclLat, distance); // Corrections for internal coordinate system - eclLat -= (PI/2); - eclLong += PI; + eclLat -= (celestia::numbers::pi/2); + eclLong += celestia::numbers::pi; return Vector3d(cos(eclLong) * sin(eclLat) * distance, cos(eclLat) * distance, @@ -503,8 +504,8 @@ class VenusOrbit : public CachingOrbit eclLong, eclLat, distance); //Corrections for internal coordinate system - eclLat -= (PI/2); - eclLong += PI; + eclLat -= (celestia::numbers::pi/2); + eclLong += celestia::numbers::pi; return Vector3d(cos(eclLong) * sin(eclLat) * distance, cos(eclLat) * distance, @@ -566,12 +567,12 @@ class EarthOrbit : public CachingOrbit dr = 5.43e-06*sin(a1)+1.575e-05*sin(b1)+1.627e-05*sin(c1)+ 3.076e-05*cos(d1)+9.27e-06*sin(h1); - eclLong = nu+degToRad(ls-ms+dl) + PI; + eclLong = nu+degToRad(ls-ms+dl) + celestia::numbers::pi; eclLong = pfmod(eclLong, TWOPI); distance = KM_PER_AU * (1.0000002*(1-s*cos(ea))+dr); // Correction for internal coordinate system - eclLong += PI; + eclLong += celestia::numbers::pi; return Vector3d(-cos(eclLong) * distance, 0, @@ -733,8 +734,8 @@ class LunarOrbit : public CachingOrbit EpochConvert(jd, astro::J2000, RA, dec, RA, dec); // Corrections for internal coordinate system - dec -= (PI/2); - RA += PI; + dec -= (celestia::numbers::pi/2); + RA += celestia::numbers::pi; return Vector3d(cos(RA) * sin(dec) * distance, cos(dec) * distance, @@ -826,8 +827,8 @@ class MarsOrbit : public CachingOrbit eclLong, eclLat, distance); //Corrections for internal coordinate system - eclLat -= (PI/2); - eclLong += PI; + eclLat -= (celestia::numbers::pi/2); + eclLong += celestia::numbers::pi; return Vector3d(cos(eclLong) * sin(eclLat) * distance, cos(eclLat) * distance, @@ -937,8 +938,8 @@ class JupiterOrbit : public CachingOrbit eclLong, eclLat, distance); //Corrections for internal coordinate system - eclLat -= (PI/2); - eclLong += PI; + eclLat -= (celestia::numbers::pi/2); + eclLong += celestia::numbers::pi; return Vector3d(cos(eclLong) * sin(eclLat) * distance, cos(eclLat) * distance, @@ -1075,8 +1076,8 @@ class SaturnOrbit : public CachingOrbit eclLong, eclLat, distance); //Corrections for internal coordinate system - eclLat -= (PI/2); - eclLong += PI; + eclLat -= (celestia::numbers::pi/2); + eclLong += celestia::numbers::pi; return Vector3d(cos(eclLong) * sin(eclLat) * distance, cos(eclLat) * distance, @@ -1168,8 +1169,8 @@ class UranusOrbit : public CachingOrbit eclLong, eclLat, distance); //Corrections for internal coordinate system - eclLat -= (PI/2); - eclLong += PI; + eclLat -= (celestia::numbers::pi/2); + eclLong += celestia::numbers::pi; return Vector3d(cos(eclLong) * sin(eclLat) * distance, cos(eclLat) * distance, @@ -1251,8 +1252,8 @@ class NeptuneOrbit : public CachingOrbit eclLong, eclLat, distance); //Corrections for internal coordinate system - eclLat -= (PI/2); - eclLong += PI; + eclLat -= (celestia::numbers::pi/2); + eclLong += celestia::numbers::pi; return Vector3d(cos(eclLong) * sin(eclLat) * distance, cos(eclLat) * distance, @@ -1296,8 +1297,8 @@ class PlutoOrbit : public CachingOrbit eclLong, eclLat, distance); //Corrections for internal coordinate system - eclLat -= PI / 2; - eclLong += PI; + eclLat -= celestia::numbers::pi / 2; + eclLong += celestia::numbers::pi; return Vector3d(cos(eclLong) * sin(eclLat) * distance, cos(eclLat) * distance, @@ -1607,8 +1608,8 @@ class IoOrbit : public CachingOrbit L += JupAscendingNode; // Corrections for internal coordinate system - B -= (PI/2); - L += PI; + B -= (celestia::numbers::pi/2); + L += celestia::numbers::pi; return Vector3d(cos(L) * sin(B) * R, cos(B) * R, @@ -1699,8 +1700,8 @@ class EuropaOrbit : public CachingOrbit L += JupAscendingNode; // Corrections for internal coordinate system - B -= (PI/2); - L += PI; + B -= (celestia::numbers::pi/2); + L += celestia::numbers::pi; return Vector3d(cos(L) * sin(B) * R, cos(B) * R, @@ -1794,8 +1795,8 @@ class GanymedeOrbit : public CachingOrbit L += JupAscendingNode; //Corrections for internal coordinate system - B -= (PI/2); - L += PI; + B -= (celestia::numbers::pi/2); + L += celestia::numbers::pi; return Vector3d(cos(L) * sin(B) * R, cos(B) * R, @@ -1934,8 +1935,8 @@ class CallistoOrbit : public CachingOrbit L += JupAscendingNode; //Corrections for internal coordinate system - B -= (PI/2); - L += PI; + B -= (celestia::numbers::pi/2); + L += celestia::numbers::pi; return Vector3d(cos(L) * sin(B) * R, cos(B) * R, @@ -2575,7 +2576,7 @@ class UranianSatelliteOrbit : public CachingOrbit double getPeriod() const override { - return 2 * PI / n; + return 2 * celestia::numbers::pi / n; } double getBoundingRadius() const override diff --git a/src/celephem/customrotation.cpp b/src/celephem/customrotation.cpp index b0f743d28..166cb49eb 100644 --- a/src/celephem/customrotation.cpp +++ b/src/celephem/customrotation.cpp @@ -13,6 +13,7 @@ #include "customrotation.h" #include "rotation.h" #include "precession.h" +#include #include #include #include @@ -78,7 +79,7 @@ public: double inclination = 90.0 - poleDec; if (flipped) - return XRotation(PI) * XRotation(degToRad(-inclination)) * YRotation(degToRad(-node)); + return XRotation(celestia::numbers::pi) * XRotation(degToRad(-inclination)) * YRotation(degToRad(-node)); else return XRotation(degToRad(-inclination)) * YRotation(degToRad(-node)); } @@ -122,7 +123,7 @@ public: { // TODO: Use a more accurate model for sidereal time double t = tjd - astro::J2000; - double theta = 2 * PI * (t * 24.0 / 23.9344694 - 259.853 / 360.0); + double theta = 2 * celestia::numbers::pi * (t * 24.0 / 23.9344694 - 259.853 / 360.0); return YRotation(-theta); } @@ -147,8 +148,8 @@ public: // P and Q: // P = sin(pi)*sin(Pi) // Q = sin(pi)*cos(Pi) - double P = pole.PA * 2.0 * PI / 1296000; - double Q = pole.QA * 2.0 * PI / 1296000; + double P = pole.PA * 2.0 * celestia::numbers::pi / 1296000; + double Q = pole.QA * 2.0 * celestia::numbers::pi / 1296000; double piA = asin(sqrt(P * P + Q * Q)); double PiA = atan2(P, Q); @@ -161,7 +162,7 @@ public: Quaterniond q = XRotation(obliquity) * ZRotation(-precession) * eclRotation.conjugate(); // convert to Celestia's coordinate system - return XRotation(PI / 2.0) * q * XRotation(-PI / 2.0); + return XRotation(celestia::numbers::pi / 2.0) * q * XRotation(-celestia::numbers::pi / 2.0); } double getPeriod() const override diff --git a/src/celephem/orbit.cpp b/src/celephem/orbit.cpp index 0af68897d..ab877ce9f 100644 --- a/src/celephem/orbit.cpp +++ b/src/celephem/orbit.cpp @@ -9,6 +9,7 @@ // of the License, or (at your option) any later version. #include "orbit.h" +#include #include #include #include @@ -373,7 +374,7 @@ Vector3d EllipticalOrbit::velocityAtE(double E) const double sinE = sin(E); double cosE = cos(E); - double meanMotion = 2.0 * PI / period; + double meanMotion = 2.0 * celestia::numbers::pi / period; double edot = meanMotion / (1 - eccentricity * cosE); x = -a * sinE * edot; @@ -403,7 +404,7 @@ Vector3d EllipticalOrbit::velocityAtE(double E) const Vector3d EllipticalOrbit::positionAtTime(double t) const { t = t - epoch; - double meanMotion = 2.0 * PI / period; + double meanMotion = 2.0 * celestia::numbers::pi / period; double meanAnomaly = meanAnomalyAtEpoch + t * meanMotion; double E = eccentricAnomaly(meanAnomaly); @@ -414,7 +415,7 @@ Vector3d EllipticalOrbit::positionAtTime(double t) const Vector3d EllipticalOrbit::velocityAtTime(double t) const { t = t - epoch; - double meanMotion = 2.0 * PI / period; + double meanMotion = 2.0 * celestia::numbers::pi / period; double meanAnomaly = meanAnomalyAtEpoch + t * meanMotion; double E = eccentricAnomaly(meanAnomaly); @@ -543,7 +544,7 @@ static EllipticalOrbit* StateVectorToOrbit(const Vector3d& position, double om = atan2(P.y(), Q.y()); // Compute the period - double T = 2 * PI * sqrt(cube(a) / GM); + double T = 2 * celestia::numbers::pi * sqrt(cube(a) / GM); T = T / 86400.0; // Convert from seconds to days return new EllipticalOrbit(a * (1 - e), e, i, Om, om, M, T, t); diff --git a/src/celephem/precession.cpp b/src/celephem/precession.cpp index e30f19860..51a6ee3d9 100644 --- a/src/celephem/precession.cpp +++ b/src/celephem/precession.cpp @@ -12,6 +12,7 @@ #include #include +#include #include #include "precession.h" @@ -108,7 +109,7 @@ astro::EclipticPrecession_P03LP(double T) for (unsigned int i = 0; i < nTerms; i++) { const EclipticPrecessionTerm& p = EclipticPrecessionTerms[i]; - double theta = 2.0 * PI * T / p.period; + double theta = 2.0 * celestia::numbers::pi * T / p.period; double s = sin(theta); double c = cos(theta); pole.PA += p.Pc * c + p.Ps * s; @@ -146,7 +147,7 @@ astro::PrecObliquity_P03LP(double T) for (unsigned int i = 0; i < nTerms; i++) { const PrecessionTerm& p = PrecessionTerms[i]; - double theta = 2.0 * PI * T / p.period; + double theta = 2.0 * celestia::numbers::pi * T / p.period; double s = sin(theta); double c = cos(theta); angles.pA += p.pc * c + p.ps * s; diff --git a/src/celephem/rotation.cpp b/src/celephem/rotation.cpp index 12fdfd9b9..ed0a6a96b 100644 --- a/src/celephem/rotation.cpp +++ b/src/celephem/rotation.cpp @@ -12,6 +12,7 @@ // of the License, or (at your option) any later version. #include "rotation.h" +#include #include #include #include @@ -212,7 +213,7 @@ UniformRotationModel::spin(double tjd) const // the texture. remainder += 0.5; - return YRotation(-remainder * 2 * PI - offset); + return YRotation(-remainder * 2 * celestia::numbers::pi - offset); } @@ -227,7 +228,7 @@ Vector3d UniformRotationModel::angularVelocityAtTime(double tdb) const { Vector3d v = equatorOrientationAtTime(tdb).conjugate() * Vector3d::UnitY();; - return v * (2.0 * PI / period); + return v * (2.0 * celestia::numbers::pi / period); } @@ -276,7 +277,7 @@ PrecessingRotationModel::spin(double tjd) const // the texture. remainder += 0.5; - return YRotation(-remainder * 2 * PI - offset); + return YRotation(-remainder * 2 * celestia::numbers::pi - offset); } @@ -293,7 +294,7 @@ PrecessingRotationModel::equatorOrientationAtTime(double tjd) const else { nodeOfDate = (double) ascendingNode - - (2.0 * PI / precessionPeriod) * (tjd - epoch); + (2.0 * celestia::numbers::pi / precessionPeriod) * (tjd - epoch); } return XRotation((double) -inclination) * YRotation(-nodeOfDate); diff --git a/src/celephem/samporient.cpp b/src/celephem/samporient.cpp index 9aa3b2b36..b4a9053a7 100644 --- a/src/celephem/samporient.cpp +++ b/src/celephem/samporient.cpp @@ -12,6 +12,7 @@ // of the License, or (at your option) any later version. #include "samporient.h" +#include #include #include #include @@ -59,7 +60,7 @@ typedef vector OrientationSampleVector; // 90 degree rotation about x-axis to convert orientation to Celestia's // coordinate system. -static Quaternionf coordSysCorrection = XRotation((float) (PI / 2.0)); +static Quaternionf coordSysCorrection = XRotation((float) (celestia::numbers::pi / 2.0)); bool operator<(const OrientationSample& a, const OrientationSample& b) diff --git a/src/celephem/spicerotation.cpp b/src/celephem/spicerotation.cpp index 3d5f9b1f3..cd69fe77d 100644 --- a/src/celephem/spicerotation.cpp +++ b/src/celephem/spicerotation.cpp @@ -12,6 +12,7 @@ #include "spicerotation.h" #include "spiceinterface.h" +#include #include #include #include @@ -27,8 +28,8 @@ using namespace celmath; using celestia::util::GetLogger; static const double MILLISEC = astro::secsToDays(0.001); -static const Quaterniond Rx90 = XRotation(PI / 2.0); -static const Quaterniond Ry180 = YRotation(PI); +static const Quaterniond Rx90 = XRotation(celestia::numbers::pi / 2.0); +static const Quaterniond Ry180 = YRotation(celestia::numbers::pi); /*! Create a new rotation model based on a SPICE frame. The * orientation of the rotation model is the orientation of the diff --git a/src/celephem/vsop87.cpp b/src/celephem/vsop87.cpp index 2796b3a5a..60a8f07c2 100644 --- a/src/celephem/vsop87.cpp +++ b/src/celephem/vsop87.cpp @@ -13,6 +13,7 @@ // of the License, or (at your option) any later version. #include +#include #include #include #include "vsop87.h" @@ -11088,8 +11089,8 @@ class VSOP87Orbit : public CachingOrbit r *= KM_PER_AU; // Corrections for internal coordinate system - b -= PI / 2; - l += PI; + b -= celestia::numbers::pi / 2; + l += celestia::numbers::pi; return Vector3d(cos(l) * sin(b) * r, cos(b) * r, diff --git a/src/celestia/celestiacore.cpp b/src/celestia/celestiacore.cpp index 24db2beab..48e8c46ce 100644 --- a/src/celestia/celestiacore.cpp +++ b/src/celestia/celestiacore.cpp @@ -16,6 +16,7 @@ #include "celestiacore.h" #include "favorites.h" #include "url.h" +#include #include #include #include @@ -2708,15 +2709,15 @@ static void displayApparentMagnitude(Overlay& overlay, static void displayRADec(Overlay& overlay, const Vector3d& v) { - double phi = atan2(v.x(), v.z()) - PI / 2; + double phi = atan2(v.x(), v.z()) - celestia::numbers::pi / 2; if (phi < 0) - phi = phi + 2 * PI; + phi = phi + 2 * celestia::numbers::pi; double theta = atan2(sqrt(v.x() * v.x() + v.z() * v.z()), v.y()); if (theta > 0) - theta = PI / 2 - theta; + theta = celestia::numbers::pi / 2 - theta; else - theta = -PI / 2 - theta; + theta = -celestia::numbers::pi / 2 - theta; displayRightAscension(overlay, radToDeg(phi)); diff --git a/src/celestia/gtk/dialog-eclipse.cpp b/src/celestia/gtk/dialog-eclipse.cpp index 0ed13c6d2..5d147be4d 100644 --- a/src/celestia/gtk/dialog-eclipse.cpp +++ b/src/celestia/gtk/dialog-eclipse.cpp @@ -12,6 +12,7 @@ #include +#include #include #include #include @@ -328,7 +329,7 @@ static gint eclipseGoto(GtkButton*, EclipseData* ed) double distance = target.radius() * 4.0; sim->gotoLocation(UniversalCoord::Zero().offsetKm(Vector3d::UnitX() * distance), - (YRotation(-PI / 2) * XRotation(-PI / 2)), 2.5); + (YRotation(-celestia::numbers::pi / 2) * XRotation(-celestia::numbers::pi / 2)), 2.5); return TRUE; } diff --git a/src/celestia/qt/qtinfopanel.cpp b/src/celestia/qt/qtinfopanel.cpp index 201e53ca3..b4a661397 100644 --- a/src/celestia/qt/qtinfopanel.cpp +++ b/src/celestia/qt/qtinfopanel.cpp @@ -10,7 +10,7 @@ // as published by the Free Software Foundation; either version 2 // of the License, or (at your option) any later version. - +#include #include #include #include @@ -276,7 +276,7 @@ Vector3d rectToSpherical(const Vector3d& v) double r = v.norm(); double theta = atan2(v.y(), v.x()); if (theta < 0) - theta = theta + 2 * PI; + theta = theta + 2 * celestia::numbers::pi; double phi = asin(v.z() / r); return Vector3d(theta, phi, r); @@ -407,7 +407,7 @@ static void StateVectorToElements(const Vector3d& position, double om = atan2(P.y(), Q.y()); // Compute the period - double T = 2 * PI * sqrt(cube(a) / GM); + double T = 2 * celestia::numbers::pi * sqrt(cube(a) / GM); elements->semimajorAxis = a; elements->eccentricity = e; diff --git a/src/celestia/win32/wineclipses.cpp b/src/celestia/win32/wineclipses.cpp index 7c125b1af..9accfb7d1 100644 --- a/src/celestia/win32/wineclipses.cpp +++ b/src/celestia/win32/wineclipses.cpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include "wineclipses.h" #include "res/resource.h" @@ -383,7 +384,7 @@ BOOL APIENTRY EclipseFinderProc(HWND hDlg, double distance = target.radius() * 4.0; sim->gotoLocation(UniversalCoord::Zero().offsetKm(Vector3d::UnitX() * distance), - YRotation(-PI / 2) * XRotation(-PI / 2), + YRotation(-celestia::numbers::pi / 2) * XRotation(-celestia::numbers::pi / 2), 2.5); } break; diff --git a/src/celmath/mathlib.h b/src/celmath/mathlib.h index 70c212095..dba70c395 100644 --- a/src/celmath/mathlib.h +++ b/src/celmath/mathlib.h @@ -10,9 +10,10 @@ #pragma once #include + #include -#define PI 3.14159265358979323846 +#include namespace celmath { @@ -32,12 +33,14 @@ template constexpr T lerp(T t, T a, T b) template inline constexpr T degToRad(T d) { - return d / static_cast(180) * static_cast(PI); + using celestia::numbers::pi_v; + return d / static_cast(180) * pi_v; } template inline constexpr T radToDeg(T r) { - return r * static_cast(180) / static_cast(PI); + using celestia::numbers::inv_pi_v; + return r * static_cast(180) * inv_pi_v; } template inline constexpr T square(T x) @@ -72,12 +75,14 @@ template T pfmod(T x, T y) template inline constexpr T circleArea(T r) { - return static_cast(PI) * r * r; + using celestia::numbers::pi_v; + return pi_v * r * r; } template inline constexpr T sphereArea(T r) { - return static_cast(4) * static_cast(PI) * r * r; + using celestia::numbers::pi_v; + return static_cast(4) * pi_v * r * r; } template static Eigen::Matrix diff --git a/src/celmath/randutils.h b/src/celmath/randutils.h index d3f222a5b..08e7a24eb 100644 --- a/src/celmath/randutils.h +++ b/src/celmath/randutils.h @@ -5,6 +5,7 @@ #include +#include #include "mathlib.h" namespace celmath @@ -36,7 +37,7 @@ RealDists::SignedUnit{static_cast(-1), static_cast(1)}; template std::uniform_real_distribution -RealDists::SignedFullAngle{static_cast(-PI), static_cast(PI)}; +RealDists::SignedFullAngle{-celestia::numbers::pi_v, celestia::numbers::pi_v}; template Eigen::Matrix randomOnCircle(RNG&& rng) diff --git a/src/tools/atmosphere/scattersim.cpp b/src/tools/atmosphere/scattersim.cpp index 081ae473e..00b7e1649 100644 --- a/src/tools/atmosphere/scattersim.cpp +++ b/src/tools/atmosphere/scattersim.cpp @@ -19,6 +19,7 @@ #include #include +#include #include #include #include @@ -31,6 +32,7 @@ using namespace std; using namespace Eigen; using namespace celmath; +using celestia::numbers::pi; // Extinction lookup table dimensions constexpr const unsigned int ExtinctionLUTHeightSteps = 256; @@ -141,8 +143,8 @@ public: } else { - double phi = -viewportY * fov / 2.0 + PI / 2.0; - double theta = viewportX * fov / 2.0 + PI / 2.0; + double phi = -viewportY * fov / 2.0 + pi / 2.0; + double theta = viewportX * fov / 2.0 + pi / 2.0; viewDir.x() = sin(phi) * cos(theta); viewDir.y() = cos(phi); viewDir.z() = sin(phi) * sin(theta); @@ -153,7 +155,7 @@ public: return transformRay(viewRay, transform); } - double fov{PI / 2.0}; + double fov{pi / 2.0}; double front{1.0}; Matrix4d transform{Matrix4d::Identity()}; CameraType type{Planar}; @@ -496,9 +498,9 @@ void DumpLUT(const LUT3& lut, const string& filename) { unsigned int x = l; Vector4d v = lut.getValue(x, y, z); - Color c(mapColor(v.x() * 0.000617f * 4.0 * PI), - mapColor(v.y() * 0.00109f * 4.0 * PI), - mapColor(v.z() * 0.00195f * 4.0 * PI)); + Color c(mapColor(v.x() * 0.000617f * 4.0 * pi), + mapColor(v.y() * 0.00109f * 4.0 * pi), + mapColor(v.z() * 0.00195f * 4.0 * pi)); img.setPixel(j * tileWidth + k, i * tileHeight + l, c); } } @@ -930,8 +932,8 @@ Vector3d integrateInscattering(const Scene& scene, // Sum the optical depths to get the depth on the complete path from sun // to sample point to eye. OpticalDepths totalDepth = sumOpticalDepths(sunDepth, eyeDepth); - totalDepth.rayleigh *= 4.0 * PI; - totalDepth.mie *= 4.0 * PI; + totalDepth.rayleigh *= 4.0 * pi; + totalDepth.mie *= 4.0 * pi; Vector3d extinction = scene.atmosphere.computeExtinction(totalDepth); @@ -988,8 +990,8 @@ Vector4d integrateInscatteringFactors(const Scene& scene, // Sum the optical depths to get the depth on the complete path from sun // to sample point to eye. OpticalDepths totalDepth = sumOpticalDepths(sunDepth, eyeDepth); - totalDepth.rayleigh *= 4.0 * PI; - totalDepth.mie *= 4.0 * PI; + totalDepth.rayleigh *= 4.0 * pi; + totalDepth.mie *= 4.0 * pi; Vector3d extinction = scene.atmosphere.computeExtinction(totalDepth); @@ -1056,8 +1058,8 @@ buildExtinctionLUT(const Scene& scene) OpticalDepths depth = integrateOpticalDepth(scene, atmStart, ray.pointAt(dist)); - depth.rayleigh *= 4.0 * PI; - depth.mie *= 4.0 * PI; + depth.rayleigh *= 4.0 * pi; + depth.mie *= 4.0 * pi; Vector3d ext = scene.atmosphere.computeExtinction(depth); lut->setValue(i, j, ext.cwiseMax(1.0e-18)); @@ -1116,8 +1118,8 @@ buildOpticalDepthLUT(const Scene& scene) OpticalDepths depth = integrateOpticalDepth(scene, atmStart, ray.pointAt(dist)); - depth.rayleigh *= 4.0 * PI; - depth.mie *= 4.0 * PI; + depth.rayleigh *= 4.0 * pi; + depth.mie *= 4.0 * pi; lut->setValue(i, j, Vector3d(depth.rayleigh, depth.mie, depth.absorption)); } @@ -1390,8 +1392,8 @@ Color getPlanetColor(const Scene& scene, const Vector3d& p) // Give the planet a checkerboard texture double phi = atan2(n.z(), n.x()); double theta = asin(n.y()); - int tx = (int) (8 + 8 * phi / PI); - int ty = (int) (8 + 8 * theta / PI); + int tx = (int) (8 + 8 * phi / pi); + int ty = (int) (8 + 8 * theta / pi); return ((tx ^ ty) & 0x1) != 0 ? scene.planetColor : scene.planetColor2; } @@ -1442,8 +1444,8 @@ Color Scene::raytrace(const Eigen::ParametrizedLine& ray) const OpticalDepths sunDepth = integrateOpticalDepth(*this, surfacePt, sunRay.pointAt(sunDist)); OpticalDepths eyeDepth = integrateOpticalDepth(*this, atmStart, surfacePt); OpticalDepths totalDepth = sumOpticalDepths(sunDepth, eyeDepth); - totalDepth.rayleigh *= 4.0 * PI; - totalDepth.mie *= 4.0 * PI; + totalDepth.rayleigh *= 4.0 * pi; + totalDepth.mie *= 4.0 * pi; Vector3d extinction = atmosphere.computeExtinction(totalDepth); // Reflected color of planet surface is: @@ -1453,7 +1455,7 @@ Color Scene::raytrace(const Eigen::ParametrizedLine& ray) const atmEnd = ray.origin() + dist * ray.direction(); } - Vector3d inscatter = integrateInscattering(*this, atmStart, atmEnd) * 4.0 * PI; + Vector3d inscatter = integrateInscattering(*this, atmStart, atmEnd) * 4.0 * pi; return Color((float) inscatter.x(), (float) inscatter.y(), (float) inscatter.z()) + baseColor; @@ -1537,7 +1539,7 @@ Scene::raytrace_LUT(const Eigen::ParametrizedLine& ray) const if (LUTUsage == UseExtinctionLUT) { bool hitPlanet = hit && planetEnter > 0.0; - inscatter = integrateInscattering_LUT(*this, atmStart, atmEnd, eyePt, hitPlanet) * 4.0 * PI; + inscatter = integrateInscattering_LUT(*this, atmStart, atmEnd, eyePt, hitPlanet) * 4.0 * pi; //if (!hit) //inscatter = Vector3d::Zero(); } @@ -1579,7 +1581,7 @@ Scene::raytrace_LUT(const Eigen::ParametrizedLine& ray) const const Vector3d& rayleigh = atmosphere.rayleighCoeff; double cosSunAngle = mray.direction().dot(-light.direction); inscatter = phaseRayleigh(cosSunAngle) * rayleighScatter.cwiseProduct(rayleigh); - inscatter = inscatter * 4.0 * PI; + inscatter = inscatter * 4.0 * pi; } return Color((float) inscatter.x(), (float) inscatter.y(), (float) inscatter.z()) + @@ -1952,7 +1954,7 @@ int main(int argc, char* argv[]) scene.setParameters(sceneParams); cout << "atmosphere height: " << scene.atmosphereShellHeight << '\n'; - cout << "attenuation coeffs: " << scene.atmosphere.rayleighCoeff.transpose() * 4 * PI << '\n'; + cout << "attenuation coeffs: " << scene.atmosphere.rayleighCoeff.transpose() * 4 * pi << '\n'; if (LUTUsage != NoLUT) diff --git a/src/tools/cmod/cmodsphere/cmodsphere.cpp b/src/tools/cmod/cmodsphere/cmodsphere.cpp index de22cb173..5207ba63c 100644 --- a/src/tools/cmod/cmodsphere/cmodsphere.cpp +++ b/src/tools/cmod/cmodsphere/cmodsphere.cpp @@ -10,6 +10,8 @@ #include +#include + using namespace std; @@ -18,8 +20,6 @@ using namespace std; static int latSamples = 1440; static int longSamples = 2880; -constexpr float pi = 3.14159265f; - static float* samples = nullptr; // Read a big-endian 32-bit unsigned integer @@ -130,8 +130,8 @@ void triangleSection(unsigned int subdiv, { float theta = (float) acos(w.y()); float phi = (float) atan2(-w.z(), w.x()); - float s = phi / (2.0f * pi) + 0.5f; - float t = theta / pi; + float s = phi / (2.0f * celestia::numbers::pi_v) + 0.5f; + float t = theta / celestia::numbers::pi_v; float r = sampleBilinear(samples, longSamples, latSamples, s, t); diff --git a/src/tools/stardb/startextdump.cpp b/src/tools/stardb/startextdump.cpp index 86cc6ecc2..b10937674 100644 --- a/src/tools/stardb/startextdump.cpp +++ b/src/tools/stardb/startextdump.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -257,8 +258,8 @@ bool DumpStarDatabase(istream& in, ostream& out, bool spherical) Eigen::Vector3d pos = astro::eclipticToEquatorial(eclipticpos); double distance = sqrt(x * x + y * y + z * z); // acos outputs angles in interval [0, pi], use negative sign for interval [-pi, 0] - double phi = -acos(pos.y() / distance) * 180 / PI; - double theta = atan2(pos.z(), -pos.x()) * 180 / PI; + double phi = -acos(pos.y() / distance) * 180 / celestia::numbers::pi; + double theta = atan2(pos.z(), -pos.x()) * 180 / celestia::numbers::pi; // atan2 outputs angles in interval [-pi, pi], so we add 360 to fix this double ra = theta - 180 + 360; double dec = phi + 90;