Fix stellar appearance of lunar occultation

src/core/modules/Planet.cpp

@@ -590,7 +590,7 @@
 					cometType = qc_("short-period", "type of comet");
 				else
 					cometType = qc_("periodic", "type of comet");
 
 			}
 			oss << QString("%1: <b>%2</b> (%3)<br/>").arg(q_("Type"), getObjectTypeI18n(), cometType);
 		}
@@ -1577,156 +1577,156 @@
 	return str;
 }
 
 QVariantMap Planet::getInfoMap(const StelCore *core) const
 {
 	static SolarSystem *ssystem=GETSTELMODULE(SolarSystem);
 	PlanetP earth = ssystem->getEarth();
 	const bool onEarth = (core->getCurrentPlanet()==earth);
 	QVariantMap map = StelObject::getInfoMap(core);
 
 	if (getEnglishName()!="Sun")
 	{
 		const Vec3d& observerHelioPos = core->getObserverHeliocentricEclipticPos();
 		const double hdistanceAu = getHeliocentricEclipticPos().norm();
 		const double hdistanceKm = AU * hdistanceAu;
 		map.insert("heliocentric-distance", hdistanceAu);
 		map.insert("heliocentric-distance-km", hdistanceKm);
 		float phase=getPhase(observerHelioPos);
 		map.insert("phase", phase);
 		map.insert("illumination", 100.f*phase);
 		double phaseAngle = getPhaseAngle(observerHelioPos);
 		map.insert("phase-angle", phaseAngle);
 		map.insert("phase-angle-dms", StelUtils::radToDmsStr(phaseAngle));
 		map.insert("phase-angle-deg", StelUtils::radToDecDegStr(phaseAngle));
 		const bool waning = isWaning(observerHelioPos, core->getObserverHeliocentricEclipticVelocity());
 		map.insert("is-waning", waning);
 		double elongation = getElongation(observerHelioPos);
 		map.insert("elongation", elongation);
 		map.insert("elongation-dms", StelUtils::radToDmsStr(elongation));
 		map.insert("elongation-deg", StelUtils::radToDecDegStr(elongation));
 		map.insert("velocity", getEclipticVelocity().toString());
 		map.insert("velocity-kms", QString::number(getEclipticVelocity().norm()* AU/86400., 'f', 5));
 		map.insert("heliocentric-velocity", getHeliocentricEclipticVelocity().toString());
 		map.insert("heliocentric-velocity-kms", QString::number(getHeliocentricEclipticVelocity().norm()* AU/86400., 'f', 5));
 		map.insert("scale", sphereScale);
 		map.insert("albedo", getAlbedo());
 	}
 	else
 	{
 		QPair<double, PlanetP> eclObj = ssystem->getSolarEclipseFactor(core);
 		const double eclipseObscuration = 100.*(1.-eclObj.first);
 		if (eclipseObscuration>1.e-7)
 		{
 			map.insert("eclipse-obscuration", eclipseObscuration);
 			PlanetP obj = eclObj.second;
 			if (core->getCurrentPlanet()==ssystem->getEarth() && obj==ssystem->getMoon())
 			{
 				double angularSize = 2.*getAngularRadius(core)*M_PI_180;
 				const double eclipseMagnitude = (0.5*angularSize + (obj->getAngularRadius(core)*M_PI_180)/obj->getSphereScale() - getJ2000EquatorialPos(core).angle(obj->getJ2000EquatorialPos(core)))/angularSize;
 				map.insert("eclipse-magnitude", eclipseMagnitude);
 			}
 			else
 				map.insert("eclipse-magnitude", 0.0);
 		}
 		else
 		{
 			map.insert("eclipse-obscuration", 0.0);
 			map.insert("eclipse-magnitude", 0.0);
 		}
 	}
 	map.insert("type", getType());
 	map.insert("object-type", getObjectType());
 	const double distanceAu = getJ2000EquatorialPos(core).norm();
 	const double distanceKm = AU * distanceAu;
 	map.insert("distance", distanceAu);
 	map.insert("distance-km", distanceKm);
 
 	if (onEarth)
 	{
 		if (getEnglishName()!="Sun")
 		{
 			QPair<Vec4d, Vec3d>phys=getSubSolarObserverPoints(core);
 			map.insert("central_l", phys.first[2]*M_180_PI);
 			map.insert("central_b", phys.first[1]*M_180_PI);
 			map.insert("pa_axis", phys.first[3]*M_180_PI);
 			map.insert("subsolar_l", phys.second[2]*M_180_PI);
 			map.insert("subsolar_b", phys.second[1]*M_180_PI);
 			// some users require not "modern elongation" but just the DeltaLambda (GH:#1786)
 			double raSun, deSun, ra, de, lSun, ecLong, bSun, ecLat;
 			double obl=earth->getRotObliquity(core->getJDE());
 			if (core->getUseNutation())
 			{
 				double dEps, dPsi;
 				getNutationAngles(core->getJDE(), &dPsi, &dEps);
 				obl+=dEps;
 			}
 			StelUtils::rectToSphe(&raSun, &deSun, ssystem->getSun()->getEquinoxEquatorialPos(core));
 			StelUtils::rectToSphe(&ra, &de, getEquinoxEquatorialPos(core));
 			StelUtils::equToEcl(raSun, deSun, obl, &lSun, &bSun);
 			StelUtils::equToEcl(ra, de, obl, &ecLong, &ecLat);
 			double elongAlongEcliptic = StelUtils::fmodpos(ecLong-lSun, M_PI*2.);
 			if (elongAlongEcliptic > M_PI) elongAlongEcliptic-=2.*M_PI;
 			map.insert("ecl-elongation", elongAlongEcliptic);
 			map.insert("ecl-elongation-dms", StelUtils::radToDmsStr(elongAlongEcliptic));
 			map.insert("ecl-elongation-deg", StelUtils::radToDecDegStr(elongAlongEcliptic));
 
 			if (getEnglishName()=="Moon")
 			{
 				map.insert("libration_l", -phys.first[2]*M_180_PI); // longitude counted the other way!
 				map.insert("libration_b", phys.first[1]*M_180_PI);
 				map.insert("colongitude", StelUtils::fmodpos(450.0+phys.second[2]*M_PI_180, 360.));
 
 				QPair<double,double> magnitudes = getLunarEclipseMagnitudes();
 				map.insert("penumbral-eclipse-magnitude", magnitudes.first);
 				map.insert("umbral-eclipse-magnitude", magnitudes.second);
 
 				// For computing the Moon age we use geocentric coordinates
 				StelCore* core1 = StelApp::getInstance().getCore(); // we need non-const reference here.
 				const bool useTopocentric = core1->getUseTopocentricCoordinates();
 				core1->setUseTopocentricCoordinates(false);
 				core1->update(0); // enforce update cache!
 				const double eclJDE = earth->getRotObliquity(core1->getJDE());
 				double ra_equ, dec_equ, lambdaMoon, lambdaSun, betaMoon, betaSun, raSun, deSun;
 				StelUtils::rectToSphe(&ra_equ,&dec_equ, getEquinoxEquatorialPos(core1));
 				StelUtils::equToEcl(ra_equ, dec_equ, eclJDE, &lambdaMoon, &betaMoon);
 				StelUtils::rectToSphe(&raSun,&deSun, ssystem->getSun()->getEquinoxEquatorialPos(core1));
 				StelUtils::equToEcl(raSun, deSun, eclJDE, &lambdaSun, &betaSun);
 				core1->setUseTopocentricCoordinates(useTopocentric);
 				core1->update(0); // enforce update cache to avoid odd selection of Moon details!
 				const double deltaLong = StelUtils::fmodpos((lambdaMoon-lambdaSun)*M_180_PI, 360.);
 				QString moonPhase = "";
 				if ((deltaLong<0.5) || (deltaLong>359.5))
 					moonPhase = qc_("New Moon", "Moon phase");
 				else if (deltaLong<89.5)
 					moonPhase = qc_("Waxing Crescent", "Moon phase");
 				else if (deltaLong<90.5)
 					moonPhase = qc_("First Quarter", "Moon phase");
 				else if (deltaLong<179.5)
 					moonPhase = qc_("Waxing Gibbous", "Moon phase");
 				else if (deltaLong<180.5)
 					moonPhase = qc_("Full Moon", "Moon phase");
 				else if (deltaLong<269.5)
 					moonPhase = qc_("Waning Gibbous", "Moon phase");
 				else if (deltaLong<270.5)
 					moonPhase = qc_("Third Quarter", "Moon phase");
 				else if (deltaLong<359.5)
 					moonPhase = qc_("Waning Crescent", "Moon phase");
 				else
 				{
 					qWarning() << "ERROR IN PHASE STRING PROGRAMMING!";
 					Q_ASSERT(0);
 				}
 				map.insert("phase-name", moonPhase);
 
 				const double age = deltaLong*29.530588853/360.;
 				map.insert("age", QString::number(age, 'f', 2));
 			}
 		}
 	}
 
 	return map;
 }
 
 QPair<double,double> Planet::getLunarEclipseMagnitudes() const
 {
 	QPair<double,double> magnitudes;
@@ -3019,7 +3019,7 @@
 double Planet::getAngularRadius(const StelCore* core) const
 {
 	const double rad = (rings ? rings->getSize() : equatorialRadius);
-	return std::atan2(rad*sphereScale,getJ2000EquatorialPos(core).norm()) * M_180_PI;
+	return std::asin(rad*sphereScale/getJ2000EquatorialPos(core).norm()) * M_180_PI;
 }
 
 
@@ -4937,53 +4937,53 @@
 	return orbColor;
 }
 
 void Planet::computeOrbit()
 {
 	double dateJDE = lastJDE;
 	// For open Kepler orbits, compute only the static segment around epoch which was defined by orbit_good!
 	KeplerOrbit *keplerOrbit=static_cast<KeplerOrbit*>(orbitPtr);
 	if (keplerOrbit && !closeOrbit)
 		dateJDE = keplerOrbit->getEpochJDE();
 	double calc_date;
 	Vec3d parentPos;
 	if (parent)
 		parentPos = parent->getHeliocentricEclipticPos(dateJDE)+ parent->getAberrationPush(); // aberrationPush is not strictly correct, but helps a lot...
 
 	if (keplerOrbit && keplerOrbit->getEccentricity()>0.3)
 	{
 		// TODO: compute Mstart, Mend, Estart, Eend (or nuStart, nuEnd?), compute vertices from loop over M, E or nu.
 		double Estart = keplerOrbit->eccentricAnomaly(keplerOrbit->meanAnomaly(dateJDE + (-ORBIT_SEGMENTS/2)*deltaOrbitJDE));
 		double Eend   = keplerOrbit->eccentricAnomaly(keplerOrbit->meanAnomaly(dateJDE + ( ORBIT_SEGMENTS/2)*deltaOrbitJDE));
 		if (closeOrbit && Estart<0. && Eend<0.)
 			Eend+=2.*M_PI;
 		if (closeOrbit && Estart>0. && Eend>0.)
 			Estart-=2.*M_PI;
 		if (Estart>Eend)
 			Estart-=2.*M_PI;
 		const double Eincr=(Eend-Estart)/ORBIT_SEGMENTS; // TODO: probably MOD(2pi)?
 		double E=Estart;
 		for (int d=0; d<ORBIT_SEGMENTS; d++)
 		{
 			Vec3d v;
 			keplerOrbit->positionAtEccentricAnomalyInVSOP87Coordinates(E, v.v);
 			orbit[d]=v+parentPos;
 			E += Eincr;
 		}
 	}
 	else
 		for(int d = 0; d < ORBIT_SEGMENTS; d++)
 		{
 			calc_date = dateJDE + (d-ORBIT_SEGMENTS/2)*deltaOrbitJDE;
 			// Round to a number of deltaOrbitJDE to improve caching.
 			if (d != ORBIT_SEGMENTS / 2)
 			{
 				calc_date = nearbyint(calc_date / deltaOrbitJDE) * deltaOrbitJDE;
 			}
 			orbit[d] = getEclipticPos(calc_date) + parentPos;
 		}
 }
 
 // draw orbital path of Planet. For objects on open Kepler orbits, draw orbital segment of orbit_good days around epoch.
 void Planet::drawOrbit(const StelCore* core)
 {
 	if (!static_cast<bool>(orbitFader.getInterstate()))

src/core/modules/StarMgr.cpp

@@ -44,6 +44,7 @@
 #include "ConstellationMgr.hpp"
 #include "Planet.hpp"
 #include "StelUtils.hpp"
+#include "SolarSystem.hpp"
 
 #include <QTextStream>
 #include <QFile>
@@ -1320,6 +1321,20 @@ void StarMgr::draw(StelCore* core)
 
 	// Prepare a table for storing precomputed RCMag for all ZoneArrays
 	RCMag rcmag_table[RCMAG_TABLE_SIZE];
+
+	// Try to filter away stars hidden by the Moon!
+	bool filterMoon=false;
+	PlanetP moon = GETSTELMODULE(SolarSystem)->getMoon();
+	Vec3d moonPos=moon->getJ2000EquatorialPos(core);
+	const double moonRadius = moon->getEquatorialRadius() * moon->getSphereScale();
+	double angularSize = asin(moonRadius / moonPos.norm());
+	moonPos.normalize();
+	SphericalCap moonCap(moonPos, cos(angularSize));
+	for (auto  &cap : viewportCaps)
+	{
+	    if (cap.intersects(moonCap))
+		filterMoon=true;
+	}
 
 	// Draw all the stars of all the selected zones
 	for (const auto* z : std::as_const(gridLevels))
@@ -1363,9 +1378,9 @@ void StarMgr::draw(StelCore* core)
 		int zone;
 
 		for (GeodesicSearchInsideIterator it1(*geodesic_search_result,z->level);(zone = it1.next()) >= 0;)
-			z->draw(&sPainter, zone, true, rcmag_table, limitMagIndex, core, maxMagStarName, names_brightness, viewportCaps, withAberration, velf);
+			z->draw(&sPainter, zone, true, rcmag_table, limitMagIndex, core, maxMagStarName, names_brightness, viewportCaps, withAberration, velf, filterMoon, moonCap);
 		for (GeodesicSearchBorderIterator it1(*geodesic_search_result,z->level);(zone = it1.next()) >= 0;)
-			z->draw(&sPainter, zone, false, rcmag_table, limitMagIndex, core, maxMagStarName,names_brightness, viewportCaps, withAberration, velf);
+			z->draw(&sPainter, zone, false, rcmag_table, limitMagIndex, core, maxMagStarName,names_brightness, viewportCaps, withAberration, velf, filterMoon, moonCap);
 	}
 	exit_loop:
 

src/core/modules/ZoneArray.cpp

@@ -408,7 +408,8 @@ template<class Star>
 void SpecialZoneArray<Star>::draw(StelPainter* sPainter, int index, bool isInsideViewport, const RCMag* rcmag_table,
 				  int limitMagIndex, StelCore* core, int maxMagStarName, float names_brightness,
 				  const QVector<SphericalCap> &boundingCaps,
-				  const bool withAberration, const Vec3f vel) const
+				  const bool withAberration, const Vec3f vel,
+				  const bool filterMoon, SphericalCap moonCap) const
 {
 	StelSkyDrawer* drawer = core->getSkyDrawer();
 	Vec3f vf;
@@ -474,6 +475,14 @@ void SpecialZoneArray<Star>::draw(StelPainter* sPainter, int index, bool isInsid
 				continue;
 		}
 
+		// If the star is covered by the Moon, don't draw (avoid displaying part of halo of occulted star!)
+		// TODO: This is a bad test. We should draw the moon into a stencil buffer and check that.
+		if (filterMoon)
+		{
+			if (moonCap.contains(vf))
+				continue;
+		}
+
 		int extinctedMagIndex = s->getMag();
 		float twinkleFactor=1.0f; // allow height-dependent twinkle.
 		if (withExtinction)

src/core/modules/ZoneArray.hpp

@@ -102,7 +102,8 @@ class ZoneArray
 					  const RCMag* rcmag_table, int limitMagIndex, StelCore* core,
 					  int maxMagStarName, float names_brightness,
 					  const QVector<SphericalCap>& boundingCaps,
-					  const bool withAberration, const Vec3f vel) const = 0;
+					  const bool withAberration, const Vec3f vel,
+					  const bool filterMoon, SphericalCap moonCap) const = 0;
 
 	//! Get whether or not the catalog was successfully loaded.
 	//! @return @c true if at least one zone was loaded, otherwise @c false
@@ -186,7 +187,8 @@ class SpecialZoneArray : public ZoneArray
 	          const RCMag *rcmag_table, int limitMagIndex, StelCore* core,
 	          int maxMagStarName, float names_brightness,
 	          const QVector<SphericalCap>& boundingCaps,
-	          const bool withAberration, const Vec3f vel) const override;
+	          const bool withAberration, const Vec3f vel,
+            const bool filterMoon, SphericalCap moonCap) const override;
 
 	void scaleAxis() override;
 	void searchAround(const StelCore* core, int index,const Vec3d &v,double cosLimFov,