Add parallax effect to stars (#4023)

guide/ch_astronomical_concepts.tex

@@ -1123,15 +1123,15 @@ \subsection{Stellar Parallax}
 Even the nearest stars exhibit very small movement due to
 parallax. The closest star to the Earth other than the Sun is Proxima
 Centauri. It has an annual parallax of $0.77199''$, corresponding to a
-distance of $1.295\pc$ (4.22 light years).
-
-Even with the most sensitive instruments for measuring the positions of
-the stars it is only possible to use parallax to determine the distance
-of stars up to about 1,600 light years from the Earth, after which the
-annual parallax is so small it cannot be measured accurately enough.
+distance of $1.295\pc$ (4.22 light years). Currently, Gaia satellite is the most
+accurate instrument for measuring stellar parallax. It can measure parallax
+at the precision level of sub milli-arcseconds for bright stars ($G<15$), measuring distance
+with $\approx 10\%$ to stars up to 10,000 parsecs.
 
 In Stellarium, the annual parallax can be listed in the object information for stars
-when available. It is not accounted for in the positional calculations.
+when available and is accounted for by default in the positional calculations on any 
+Solar System objects including the Earth. You can also exaggerate the parallax effect up to $10,000\times$ 
+to make it more apparent.
 
 \section{Aberration of Light}
 \label{sec:Concepts:Aberration}
@@ -1165,7 +1165,7 @@ \section{Aberration of Light}
 In Stellarium, you can observe a slight annual wobble in all objects'
 positions against the fixed celestial grids. Just enable the
 equatorial grid, zoom in and run through the months. You can also
-exaggerate the aberration effect upto $5\times$ to make it more
+exaggerate the aberration effect up to $5\times$ to make it more
 apparent. \footnote{Exaggerating even more causes graphics defects and
   has been disabled for now.}
 
@@ -1174,7 +1174,7 @@ \section{Proper Motion}
 
 \indexterm[proper motion]{Proper motion} is the change in the position of a star over time as a
 result of its motion through space relative to the Sun. It does not
-include the apparent shift in position of star due to annular parallax.
+include the apparent shift in position of star due to annual parallax.
 The star exhibiting the greatest proper motion is \indexterm{Barnard's Star} which
 moves more than ten seconds of arc per year.
 
@@ -1186,14 +1186,14 @@ \section{Proper Motion}
 
 Note however some limitations:
 \begin{enumerate}
-\item Stellarium will stop at $\pm 100.000$ years. This limit may be
+\item Stellarium will stop at $\pm 100,000$ years. This limit may be
   still suitable for most stellar locations. The planetary locations
   are not trustworthy outside of a much narrower temporal window (see
   section~\ref{ch:Accuracy}). You cannot simulate the sky over the
   dinosaurs or such things.
-\item Proper motion is only modeled by linear components. True 3D
-  motion in space requires more computation, which would slow down the
-  program.
+\item For most bright stars ($V<10.5$), their full 6D astrometry (2D sky positions, 
+  2D proper motion, parallax, radial velocity) computed by assuming straight-line motion at constant speed. 
+  Otherwise proper motion are only modelled by linear components.
 \item Double stars are listed in catalogs as two individual stars with
   their current proper motion which consists of a common spatial motion and the 
   orbital motion around their common center of gravity. However, their orbital

plugins/Exoplanets/src/Exoplanet.cpp

@@ -670,8 +670,7 @@ Vec3d Exoplanet::getJ2000EquatorialPos(const StelCore* core) const
 {
 	if ((core) && (core->getUseAberration()) && (core->getCurrentPlanet()))
 	{
-		Vec3d vel=core->getCurrentPlanet()->getHeliocentricEclipticVelocity();
-		vel=StelCore::matVsop87ToJ2000*vel*core->getAberrationFactor()*(AU/(86400.0*SPEED_OF_LIGHT));
+		const Vec3d vel = core->getAberrationVec(core->getJDE());
 		Vec3d pos=XYZ+vel;
 		pos.normalize();
 		return pos;

plugins/Novae/src/Nova.cpp

@@ -343,8 +343,7 @@ Vec3d Nova::getJ2000EquatorialPos(const StelCore* core) const
 {
 	if ((core) && (core->getUseAberration()) && (core->getCurrentPlanet()))
 	{
-		Vec3d vel=core->getCurrentPlanet()->getHeliocentricEclipticVelocity();
-		vel=StelCore::matVsop87ToJ2000*vel*core->getAberrationFactor()*(AU/(86400.0*SPEED_OF_LIGHT));
+		const Vec3d vel = core->getAberrationVec(core->getJDE());
 		Vec3d pos=XYZ+vel;
 		pos.normalize();
 		return pos;

plugins/Pulsars/src/Pulsar.cpp

@@ -465,8 +465,7 @@ Vec3d Pulsar::getJ2000EquatorialPos(const StelCore* core) const
 
 	if ((core) && (core->getUseAberration()) && (core->getCurrentPlanet()))
 	{
-		Vec3d vel=core->getCurrentPlanet()->getHeliocentricEclipticVelocity();
-		vel=StelCore::matVsop87ToJ2000*vel*core->getAberrationFactor()*(AU/(86400.0*SPEED_OF_LIGHT));
+		const Vec3d vel = core->getAberrationVec(core->getJDE());
 		Vec3d pos=v+vel;
 		pos.normalize();
 		return pos;

plugins/Quasars/src/Quasar.cpp

@@ -252,8 +252,7 @@ Vec3d Quasar::getJ2000EquatorialPos(const StelCore* core) const
 {
 	if ((core) && (core->getUseAberration()) && (core->getCurrentPlanet()))
 	{
-		Vec3d vel=core->getCurrentPlanet()->getHeliocentricEclipticVelocity();
-		vel=StelCore::matVsop87ToJ2000*vel*core->getAberrationFactor()*(AU/(86400.0*SPEED_OF_LIGHT));
+		const Vec3d vel = core->getAberrationVec(core->getJDE());
 		Vec3d pos=XYZ+vel;
 		pos.normalize();
 		return pos;

plugins/Supernovae/src/Supernova.cpp

@@ -272,8 +272,7 @@ Vec3d Supernova::getJ2000EquatorialPos(const StelCore* core) const
 {
 	if ((core) && (core->getUseAberration()) && (core->getCurrentPlanet()))
 	{
-		Vec3d vel=core->getCurrentPlanet()->getHeliocentricEclipticVelocity();
-		vel=StelCore::matVsop87ToJ2000*vel*core->getAberrationFactor()*(AU/(86400.0*SPEED_OF_LIGHT));
+		const Vec3d vel = core->getAberrationVec(core->getJDE());
 		Vec3d pos=XYZ+vel;
 		pos.normalize();
 		return pos;

src/core/StelCore.cpp

@@ -41,6 +41,7 @@
 #include "EphemWrapper.hpp"
 #include "NomenclatureItem.hpp"
 #include "precession.h"
+#include "Star.hpp"
 
 #include <QSettings>
 #include <QDebug>
@@ -71,6 +72,13 @@ const double StelCore::JD_DAY    = 1.;
 const double StelCore::ONE_OVER_JD_SECOND = 86400;		// 86400
 const double StelCore::TZ_ERA_BEGINNING = 2395996.5;		// December 1, 1847
 
+Vec3d StelCore::cachedParallaxDiff = Vec3d(0.,0.,0.);
+double StelCore::cachedParallaxJD = 0.0;
+PlanetP StelCore::cachedParallaxPlanet;
+Vec3d StelCore::cachedAberrationVec = Vec3d(0.,0.,0.);
+double StelCore::cachedAberrationJD = 0.0;
+PlanetP StelCore::cachedAberrationPlanet = Q_NULLPTR;
+
 StelCore::StelCore()
 	: skyDrawer(Q_NULLPTR)
 	, movementMgr(Q_NULLPTR)
@@ -82,6 +90,8 @@ StelCore::StelCore()
 	, flagUseNutation(true)
 	, flagUseAberration(true)
 	, aberrationFactor(1.0)
+	, flagUseParallax(true)
+	, parallaxFactor(1.0)
 	, flagUseTopocentricCoordinates(true)
 	, timeSpeed(JD_SECOND)
 	, JD(0.,0.)
@@ -141,6 +151,8 @@ StelCore::StelCore()
 	flagUseNutation=conf->value("astro/flag_nutation", true).toBool();
 	flagUseAberration=conf->value("astro/flag_aberration", true).toBool();
 	aberrationFactor=conf->value("astro/aberration_factor", 1.0).toDouble();
+	flagUseParallax=conf->value("astro/flag_parallax", true).toBool();
+	parallaxFactor=conf->value("astro/parallax_factor", 1.0).toDouble();
 	flagUseTopocentricCoordinates=conf->value("astro/flag_topocentric_coordinates", true).toBool();
 	flagUseDST=conf->value("localization/flag_dst", true).toBool();
 
@@ -159,6 +171,8 @@ StelCore::~StelCore()
 	delete geodesicGrid; geodesicGrid=Q_NULLPTR;
 	delete skyDrawer; skyDrawer=Q_NULLPTR;
 	delete position; position=Q_NULLPTR;
+	cachedParallaxPlanet=Q_NULLPTR;
+	cachedAberrationPlanet=Q_NULLPTR;
 }
 
 const QMap<QString, DitheringMode>StelCore::ditheringMap={
@@ -3091,17 +3105,58 @@ Vec3d StelCore::getMouseJ2000Pos() const
 	return mousePosition;
 }
 
+Vec3d StelCore::calculateParallaxDiff(double JD) const {
+	// ICRS coordinates are barycentric (Gaia gives barycentric RA/DEC coordinates)
+	// diff between solar system bayrcentric location at STAR_CATALOG_JDEPOCH and current solar system bayrcentric location
+	Vec3d PosNow = getCurrentPlanet()->getBarycentricEclipticPos(JD);
+	// Transform from heliocentric ecliptic to equatorial coordinates
+	PosNow = matVsop87ToJ2000.upper3x3() * -1. * PosNow;  // need to times -1 because technically it is doing (0,0,0) - PosNow
+	return PosNow;
+}
+
+Vec3d StelCore::getParallaxDiff(double JD) const {
+	// if isArtificial meaning transitioning between planets, use cache and don't recalculate because it will crash
+	if ((fuzzyEquals(JD, cachedParallaxJD, JD_SECOND) && (getCurrentPlanet() == cachedParallaxPlanet)) || (getCurrentPlanet()->getPlanetType() == Planet::isArtificial))
+	{
+		return cachedParallaxDiff;
+	}
+	cachedParallaxDiff = calculateParallaxDiff(JD); // set the cache to the new value
+	cachedParallaxJD = JD; // set the cache to the new value
+	cachedParallaxPlanet = getCurrentPlanet(); 
+    return cachedParallaxDiff;
+}
+
+Vec3d StelCore::calculateAberrationVec(double JD) const {
+	// Solar system barycentric velocity
+	Q_UNUSED(JD);
+	Vec3d vel = getCurrentPlanet()->getBarycentricEclipticVelocity();
+	vel = StelCore::matVsop87ToJ2000 * vel * (AU/(86400.0*SPEED_OF_LIGHT));
+	return vel;
+}
+
+Vec3d StelCore::getAberrationVec(double JD) const {
+	// need to recompute the aberration vector if the JD has changed or the planet has changed
+	if (fuzzyEquals(JD, cachedAberrationJD, JD_SECOND) && (getCurrentPlanet() == cachedAberrationPlanet))
+	{
+		return getAberrationFactor() * cachedAberrationVec;
+	}
+	cachedAberrationVec = StelCore::calculateAberrationVec(JD); // set the cache to the new value
+	cachedAberrationJD = JD; // set the cache to the new value
+	cachedAberrationPlanet = getCurrentPlanet();
+	return getAberrationFactor() * cachedAberrationVec;
+}
+
 QByteArray StelCore::getAberrationShader() const
 {
 	return 1+R"(
-uniform vec3 STELCORE_currentPlanetHeliocentricEclipticVelocity;
+uniform vec3 STELCORE_currentPlanetBarycentricEclipticVelocity;
 // objectDir points to the object as viewed from its comoving frame.
 // Return value represents the apparent direction to this object from a frame
 // that moves with respect to the object at slightly relativistic speeds (v<0.1c).
 // Relative error in aberration angle is about 0.5v/c.
 vec3 applyAberrationToObject(vec3 objectDir)
 {
-	vec3 velocity = STELCORE_currentPlanetHeliocentricEclipticVelocity;
+	vec3 velocity = STELCORE_currentPlanetBarycentricEclipticVelocity;
 	return normalize(objectDir + velocity);
 }
 // viewDir is the direction where the object appears to be when viewed from a
@@ -3110,7 +3165,7 @@ vec3 applyAberrationToObject(vec3 objectDir)
 // Relative error in aberration angle is about 0.5v/c.
 vec3 applyAberrationToViewDir(vec3 viewDir)
 {
-	vec3 velocity = STELCORE_currentPlanetHeliocentricEclipticVelocity;
+	vec3 velocity = STELCORE_currentPlanetBarycentricEclipticVelocity;
 	return normalize(viewDir - velocity);
 }
 )";
@@ -3121,14 +3176,11 @@ void StelCore::setAberrationUniforms(QOpenGLShaderProgram& program) const
 	Vec3d velocity;
 	if(getUseAberration())
 	{
-		const auto p = getCurrentPlanet();
-		const auto hev = p->getHeliocentricEclipticVelocity();
-		velocity = StelCore::matVsop87ToJ2000 * hev;
-		velocity *= getAberrationFactor() * (AU/(86400.0*SPEED_OF_LIGHT));
+		velocity = cachedAberrationVec;
 	}
 	else
 	{
 		velocity = Vec3d(0,0,0);
 	}
-	program.setUniformValue("STELCORE_currentPlanetHeliocentricEclipticVelocity", velocity.toQVector());
+	program.setUniformValue("STELCORE_currentPlanetBarycentricEclipticVelocity", velocity.toQVector());
 }

src/core/StelCore.hpp

@@ -26,6 +26,7 @@
 #include "StelLocation.hpp"
 #include "StelSkyDrawer.hpp"
 #include "StelPropertyMgr.hpp"
+#include "SolarSystem.hpp"
 #include "Dithering.hpp"
 #include <QString>
 #include <QStringList>
@@ -55,6 +56,8 @@ class StelCore : public QObject
 	Q_PROPERTY(bool flagUseNutation READ getUseNutation WRITE setUseNutation NOTIFY flagUseNutationChanged)
 	Q_PROPERTY(bool flagUseAberration READ getUseAberration WRITE setUseAberration NOTIFY flagUseAberrationChanged)
 	Q_PROPERTY(double aberrationFactor READ getAberrationFactor WRITE setAberrationFactor NOTIFY aberrationFactorChanged)
+	Q_PROPERTY(bool flagUseParallax READ getUseParallax WRITE setUseParallax NOTIFY flagUseParallaxChanged)
+	Q_PROPERTY(double parallaxFactor READ getParallaxFactor WRITE setParallaxFactor NOTIFY parallaxFactorChanged)
 	Q_PROPERTY(bool flagUseTopocentricCoordinates READ getUseTopocentricCoordinates WRITE setUseTopocentricCoordinates NOTIFY flagUseTopocentricCoordinatesChanged)
 	Q_PROPERTY(ProjectionType currentProjectionType READ getCurrentProjectionType WRITE setCurrentProjectionType NOTIFY currentProjectionTypeChanged)
 	//! This is just another way to access the projection type, by string instead of enum
@@ -374,6 +377,11 @@ class StelCore : public QObject
 
 	Vec3d getMouseJ2000Pos(void) const;
 
+	//! get vector used to compute parallax effect
+	Vec3d getParallaxDiff(double JD) const;
+	//! get vector used to compute aberration effect
+	Vec3d getAberrationVec(double JD) const;
+
 public slots:
 	//! Smoothly move the observer to the given location
 	//! @param target the target location
@@ -551,6 +559,16 @@ public slots:
 	QByteArray getAberrationShader() const;
 	void setAberrationUniforms(QOpenGLShaderProgram& program) const;
 
+	//! @return whether parallax effect is currently used.
+	bool getUseParallax() const {return flagUseParallax;}
+	//! Set whether you want computation and simulation of parallax effect.
+	void setUseParallax(bool use) { if (flagUseParallax != use) { flagUseParallax=use; StelApp::immediateSave("astro/flag_parallax", use); emit flagUseParallaxChanged(use); }}
+
+	//! @return parallax factor. 1 is realistic simulation, but higher values may be useful for didactic purposes.
+	double getParallaxFactor() const {return parallaxFactor;}
+	//! Set aberration factor. Values are clamped to 0...5. (Values above 5 cause graphical problems.)
+	void setParallaxFactor(double factor) { if (!fuzzyEquals(parallaxFactor, factor)) { parallaxFactor=qBound(0.,factor, 10000.); StelApp::immediateSave("astro/parallax_factor", parallaxFactor); emit parallaxFactorChanged(factor); }}
+
 	//! @return whether topocentric coordinates are currently used.
 	bool getUseTopocentricCoordinates() const {return flagUseTopocentricCoordinates;}
 	//! Set whether you want topocentric or planetocentric data
@@ -862,6 +880,10 @@ public slots:
 	void flagUseAberrationChanged(bool b);
 	//! This signal indicates a change in aberration exaggeration factor
 	void aberrationFactorChanged(double val);
+	//! This signal indicates a switch in use of parallax
+	void flagUseParallaxChanged(bool b);
+	//! This signal indicates a change in parallax exaggeration factor
+	void parallaxFactorChanged(double val);
 	//! This signal indicates a switch in use of topocentric coordinates
 	void flagUseTopocentricCoordinatesChanged(bool b);
 	//! Emitted whenever the projection type changes
@@ -943,6 +965,10 @@ private slots:
 	bool flagUseAberration;
 	// value to allow exaggerating aberration effects. 1 is natural value, stretching to e.g. 1000 may be useful for explanations.
 	double aberrationFactor;
+	// flag to indicate we want to include parallax effect
+	bool flagUseParallax;
+	// value to allow exaggerating parallax effects. 1 is natural value, stretching to e.g. 1000 may be useful for explanations.
+	double parallaxFactor;
 	// flag to indicate that we show topocentrically corrected coordinates. (Switching to false for planetocentric coordinates is new for 0.14)
 	bool flagUseTopocentricCoordinates;
 
@@ -984,6 +1010,17 @@ private slots:
 	bool de440Active;    // available and user-activated.
 	bool de441Active;    // available and user-activated.
 	QPair<int, int> minMaxEphemRange;
+
+	// Variables for caching the observer position relative to the star catalog reference frame
+	static Vec3d cachedParallaxDiff;
+    static double cachedParallaxJD; // Cached Julian Date
+	static PlanetP cachedParallaxPlanet;
+	Vec3d calculateParallaxDiff(double JD) const; // Actual calculation
+
+	// Variables for caching the aberration effect
+	static Vec3d cachedAberrationVec;
+	static double cachedAberrationJD;
+	static PlanetP cachedAberrationPlanet;
+	Vec3d calculateAberrationVec(double JD) const; // Actual calculation
 };
-
 #endif // STELCORE_HPP

src/core/StelSkyImageTile.cpp

@@ -90,8 +90,7 @@ void StelSkyImageTile::draw(StelCore* core, StelPainter& sPainter, float opacity
 	Vec3d vel(0.0);
 	if ((core) && (core->getUseAberration()) && (core->getCurrentPlanet()) && (withAberration))
 	{
-		vel=core->getCurrentPlanet()->getHeliocentricEclipticVelocity();
-		vel=StelCore::matVsop87ToJ2000*vel*core->getAberrationFactor()*(AU/(86400.0*SPEED_OF_LIGHT));
+		vel = core->getAberrationVec(core->getJDE());
 	}
 
 	const float limitLuminance = core->getSkyDrawer()->getLimitLuminance();

src/core/StelToast.cpp

@@ -144,9 +144,7 @@ void ToastTile::prepareDraw(Vec3f color)
 
 	if (core->getUseAberration())
 	{
-		Vec3d vel=core->getCurrentPlanet()->getHeliocentricEclipticVelocity();
-		vel=StelCore::matVsop87ToJ2000*vel;
-		vel*=core->getAberrationFactor() * (AU/(86400.0*SPEED_OF_LIGHT));
+		const Vec3d vel = core->getAberrationVec(core->getJDE());
 		vertexArray=QVector<Vec3d>(originalVertexArray);
 		for (int i=0; i<originalVertexArray.size(); i++)
 		{

src/core/modules/AsterismMgr.cpp

@@ -357,9 +357,6 @@ void AsterismMgr::drawNames(StelPainter& sPainter) const
 		return;
 
 	StelCore *core=StelApp::getInstance().getCore();
-	Vec3d vel=core->getCurrentPlanet()->getHeliocentricEclipticVelocity();
-	vel=StelCore::matVsop87ToJ2000*vel;
-	vel*=core->getAberrationFactor() * (AU/(86400.0*SPEED_OF_LIGHT));
 
 	sPainter.setBlending(true);
 	for (auto* asterism : asterisms)
@@ -369,6 +366,7 @@ void AsterismMgr::drawNames(StelPainter& sPainter) const
 		Vec3d XYZname=asterism->XYZname;
 		if (core->getUseAberration())
 		{
+			const Vec3d vel = core->getAberrationVec(core->getJDE());
 			XYZname.normalize();
 			XYZname+=vel;
 			XYZname.normalize();

src/core/modules/ConstellationMgr.cpp

@@ -684,9 +684,7 @@ void ConstellationMgr::draw(StelCore* core)
 	Vec3d vel(0.);
 	if (core->getUseAberration())
 	{
-		vel=core->getCurrentPlanet()->getHeliocentricEclipticVelocity();
-		vel=StelCore::matVsop87ToJ2000*vel;
-		vel*=core->getAberrationFactor() * (AU/(86400.0*SPEED_OF_LIGHT));
+		vel = core->getAberrationVec(core->getJDE());
 	}
 	drawNames(sPainter, vel);
 	drawArt(sPainter, vel);

src/core/modules/CustomObject.cpp

@@ -99,8 +99,7 @@ Vec3d CustomObject::getJ2000EquatorialPos(const StelCore* core) const
 
 	if  (!isMarker && (core && core->getUseAberration() && core->getCurrentPlanet()))
 	{
-		Vec3d vel=core->getCurrentPlanet()->getHeliocentricEclipticVelocity();
-		vel=StelCore::matVsop87ToJ2000*vel*core->getAberrationFactor()*(AU/(86400.0*SPEED_OF_LIGHT));
+		const Vec3d vel = core->getAberrationVec(core->getJDE());
 		Vec3d pos=XYZ+vel;
 		pos.normalize();
 		return pos;