Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Exclude outdated solar system objects #1880

Merged
merged 11 commits into from
Sep 23, 2021
Merged

Exclude outdated solar system objects #1880

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
27aa4ff
First attempt to exclude outdated comet data in ephemeris computations
alex-w Sep 11, 2021
0ba419a
Exclude outdated comet data in ephemeris and positions computations
alex-w Sep 11, 2021
122fe98
A small optimization
alex-w Sep 11, 2021
8889f52
Rename the method...
alex-w Sep 22, 2021
077df31
Update dev. docs
alex-w Sep 22, 2021
df76faf
Improve efficiency of checks
alex-w Sep 22, 2021
ddf7790
Notes for developers and change logic of filter
alex-w Sep 22, 2021
ffe8b3f
FIXME comments
alex-w Sep 22, 2021
927a249
Add a test for valid range of dates for Kepler orbits
gzotti Sep 22, 2021
c4e522a
Better interval query
gzotti Sep 22, 2021
7ea23d6
Fixed filtering in computing of the eclipses and shadows transits in …
alex-w Sep 23, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 12 additions & 0 deletions src/core/modules/Orbit.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -263,6 +263,18 @@ double KeplerOrbit::calculateSiderealPeriod() const
return calculateSiderealPeriod(a, centralMass);
}

Vec2d KeplerOrbit::objectDateValidRange() const
{
double min=std::numeric_limits<double>::min();
double max=std::numeric_limits<double>::max();
if (orbitGood>0)
{
min=t0-orbitGood;
max=t0+orbitGood;
}
return Vec2d(min, max);
}

GimbalOrbit::GimbalOrbit(double distance, double longitude, double latitude):
distance(distance),
longitude(longitude),
Expand Down
2 changes: 2 additions & 0 deletions src/core/modules/Orbit.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,8 @@ class KeplerOrbit : public Orbit {
virtual double getSemimajorAxis() const Q_DECL_OVERRIDE { return (e==1. ? 0. : q / (1.-e)); }
virtual double getEccentricity() const Q_DECL_OVERRIDE { return e; }
bool objectDateValid(const double JDE) const { return ((orbitGood<=0) || (fabs(t0-JDE)<orbitGood)); }
//! Return minimal and maximal JDE values where this orbit should be used. (if orbitGood is configured)
Vec2d objectDateValidRange() const;
//! Calculate sidereal period in days from semi-major axis and central mass. If SMA<=0 (hyperbolic orbit), return 0.
double calculateSiderealPeriod() const;
//! @param semiMajorAxis in AU. If SMA<=0 (hyperbolic orbit), return 0.
Expand Down
24 changes: 23 additions & 1 deletion src/core/modules/Planet.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4274,7 +4274,7 @@ void Planet::drawOrbit(const StelCore* core)
if (hidden || (pType==isObserver)) return;
if (orbitPtr && pType>=isArtificial)
{
if (!static_cast<KeplerOrbit*>(orbitPtr)->objectDateValid(lastJDE))
if (!hasValidPositionalData(lastJDE))
return;
}

Expand Down Expand Up @@ -4328,6 +4328,28 @@ void Planet::drawOrbit(const StelCore* core)
sPainter.setLineWidth(1);
}

bool Planet::hasValidPositionalData(const double JDE)
{
if (pType<isObserver)
return true;
else if (orbitPtr && pType>=isArtificial)
return static_cast<KeplerOrbit*>(orbitPtr)->objectDateValid(JDE);
else
return false;
}

Vec2d Planet::getValidPositionalDataRange()
{
double min=std::numeric_limits<double>::min();
double max=std::numeric_limits<double>::max();

if (orbitPtr && pType>=isArtificial)
{
return static_cast<KeplerOrbit*>(orbitPtr)->objectDateValidRange();
}
return Vec2d(min, max);
}

void Planet::update(int deltaTime)
{
hintFader.update(deltaTime);
Expand Down
9 changes: 9 additions & 0 deletions src/core/modules/Planet.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -222,6 +222,15 @@ class Planet : public StelObject
virtual double getAngularSize(const StelCore* core) const Q_DECL_OVERRIDE;
virtual bool hasAtmosphere(void) {return atmosphere;}
virtual bool hasHalo(void) {return halo;}
//! Returns whether planet positions are valid and useful for the current simulation time.
//! E.g. outdated orbital elements for Kepler orbits (beyond their orbit_good .ini file entries)
//! may lead to invalid positions which should better not be used.
//! @note for major planets and moons this method will always return true
bool hasValidPositionalData(const double JDE);
//! Returns JDE dates of presumably valid data for positional calculation.
//! For the major planets and moons, this is always (std::numeric_limits<double>::min(), std::numeric_limits<double>::max())
//! For planets with Keplerian orbits, this is (epoch-orbit_good, epoch+orbit_good)
Vec2d getValidPositionalDataRange();
float getAxisRotation(void) { return axisRotation;} //! return axisRotation last computed in computeTransMatrix().

///////////////////////////////////////////////////////////////////////////
Expand Down
46 changes: 33 additions & 13 deletions src/gui/AstroCalcDialog.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1195,6 +1195,9 @@ void AstroCalcDialog::currentCelestialPositions()
break;
}

if (!planet->hasValidPositionalData(JD))
passByType = false;

if (passByType && planet != core->getCurrentPlanet() && planet->getVMagnitudeWithExtinction(core) <= mag && planet->isAboveRealHorizon(core))
{
pos = planet->getJ2000EquatorialPos(core);
Expand Down Expand Up @@ -1717,6 +1720,10 @@ void AstroCalcDialog::generateEphemeris()
double JD = firstJD + i * currentStep;
core->setJD(JD);
core->update(0); // force update to get new coordinates

if (!obj->hasValidPositionalData(JD))
continue;

if (useHorizontalCoords)
{
pos = obj->getAltAzPosAuto(core);
Expand Down Expand Up @@ -3887,7 +3894,7 @@ void AstroCalcDialog::calculatePhenomena()
if (selectedObject!=planet && selectedObject->getType() != "Satellite")
{
// conjunction
fillPhenomenaTable(findClosestApproach(planet, selectedObject, startJD, stopJD, separation, PhenomenaTypeIndex::Conjuction), planet, selectedObject, PhenomenaTypeIndex::Conjuction);
fillPhenomenaTable(findClosestApproach(planet, selectedObject, startJD, stopJD, separation, PhenomenaTypeIndex::Conjunction), planet, selectedObject, PhenomenaTypeIndex::Conjunction);
// opposition
if (opposition)
fillPhenomenaTable(findClosestApproach(planet, selectedObject, startJD, stopJD, separation, PhenomenaTypeIndex::Opposition), planet, selectedObject, PhenomenaTypeIndex::Opposition);
Expand All @@ -3901,7 +3908,7 @@ void AstroCalcDialog::calculatePhenomena()
{
// conjunction
StelObjectP mObj = qSharedPointerCast<StelObject>(obj);
fillPhenomenaTable(findClosestApproach(planet, mObj, startJD, stopJD, separation, PhenomenaTypeIndex::Conjuction), planet, obj, PhenomenaTypeIndex::Conjuction);
fillPhenomenaTable(findClosestApproach(planet, mObj, startJD, stopJD, separation, PhenomenaTypeIndex::Conjunction), planet, obj, PhenomenaTypeIndex::Conjunction);
// opposition
if (opposition)
fillPhenomenaTable(findClosestApproach(planet, mObj, startJD, stopJD, separation, PhenomenaTypeIndex::Opposition), planet, obj, PhenomenaTypeIndex::Opposition);
Expand All @@ -3917,7 +3924,7 @@ void AstroCalcDialog::calculatePhenomena()
{
// conjunction
StelObjectP mObj = qSharedPointerCast<StelObject>(obj);
fillPhenomenaTable(findClosestApproach(planet, mObj, startJD, stopJD, separation, PhenomenaTypeIndex::Conjuction), planet, obj, PhenomenaTypeIndex::Conjuction);
fillPhenomenaTable(findClosestApproach(planet, mObj, startJD, stopJD, separation, PhenomenaTypeIndex::Conjunction), planet, obj, PhenomenaTypeIndex::Conjunction);
}
}
else
Expand All @@ -3927,7 +3934,7 @@ void AstroCalcDialog::calculatePhenomena()
{
// conjunction
StelObjectP mObj = qSharedPointerCast<StelObject>(obj);
fillPhenomenaTable(findClosestApproach(planet, mObj, startJD, stopJD, separation, PhenomenaTypeIndex::Conjuction), planet, obj);
fillPhenomenaTable(findClosestApproach(planet, mObj, startJD, stopJD, separation, PhenomenaTypeIndex::Conjunction), planet, obj);
}
}

Expand Down Expand Up @@ -4162,7 +4169,7 @@ void AstroCalcDialog::fillPhenomenaTable(const QMap<double, double> list, const
}

QString elongStr = "";
if (((object1 == sun || object2 == sun) && mode==PhenomenaTypeIndex::Conjuction) || (object2 == sun && mode==PhenomenaTypeIndex::Opposition))
if (((object1 == sun || object2 == sun) && mode==PhenomenaTypeIndex::Conjunction) || (object2 == sun && mode==PhenomenaTypeIndex::Opposition))
elongStr = dash;
else
{
Expand Down Expand Up @@ -4466,13 +4473,26 @@ double AstroCalcDialog::findInitialStep(double startJD, double stopJD, QStringLi
return step;
}

QMap<double, double> AstroCalcDialog::findClosestApproach(PlanetP& object1, StelObjectP& object2, double startJD, double stopJD, double maxSeparation, int mode)
QMap<double, double> AstroCalcDialog::findClosestApproach(PlanetP& object1, StelObjectP& object2, const double startJD, const double stopJD, const double maxSeparation, const int mode)
{
double dist, prevDist, step, step0;
int sgn, prevSgn = 0;
QMap<double, double> separations;
QPair<double, double> extremum;

if (object2->getType()=="Planet")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, do you also see the inefficiency?
For a simple date validity test you set the whole program state twice.

Please do the following: extend
bool Planet::hasValidPositionalData(const double jde)
{
if (pType<isObserver)
return true;
else if (orbitPtr && pType>=isArtificial)
return static_cast<KeplerOrbit*>(orbitPtr)->objectDateValid(jde);
else
return false;
}

And then you can simply test here:

if ( (!planet->hasValidPositionalData(stopJD)) || (!planet->hasValidPositionalData(startJD)) )
return separations; // empty result.

And probably this is not even what we need here? Imagine some object with orbit_good=100 (days)

You make an ephemeris for 250 days which enclose the orbit's validity range. The object would be excluded with the current solution! Maybe add a test for (startJD+stopJD)/2? Or sample every 50 days between startJD and stopJD? Or is there a limit or typical value for (stopJD-startJD)? If this is always something like 10 days, you could exclude the object if both dates are invalid.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've added JDE as parameter to reduce number of calls, but for second proposal I think it should be postponed to the future releases, because this subsystem should be refactored after pencil and paper works.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK. The inefficiency is solved.
The next question is just what this function should do. Is this logic correct? Whenever at least one limit of the search interval is outside the validity range for the planet, the planet is discarded. I see a danger of another error. If you need to merge this now, leave a FIXME with this question for later, and maybe somebody reports an error. Else please don't merge, and let's solve this with better tests for the overlapping cases. (We may need to add some really simple date validity queries to the Planet.)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please add the QWarning() logging, or some other warning mechanism. the QWarning() is not much asked, and will help if a user opens an issue; the log will explain what happened.

thank you!

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The next question is just what this function should do. Is this logic correct?

Yes, but I'll change it as first iteration of excluding the data.

{
PlanetP planet = qSharedPointerCast<Planet>(object2);
// We shouldn't compute eclipses and shadow transits not for planets and their moons
if (mode==PhenomenaTypeIndex::Shadows && object2->getEnglishName()!="Sun" && planet->getParent()!=object1)
return separations;

// If we don't have at least partial overlap between planet valid dates and our interval, skip by returning an empty map.
const Vec2d planetValidityLimits=planet->getValidPositionalDataRange();
if ( (planetValidityLimits[0] > stopJD) || (planetValidityLimits[1] < startJD) )
return separations;
}

QStringList objects;
objects.clear();
objects.append(object1->getEnglishName());
Expand Down Expand Up @@ -4590,11 +4610,11 @@ QMap<double, double> AstroCalcDialog::findGreatestElongationApproach(PlanetP& ob
step0 = findInitialStep(startJD, stopJD, objects);
step = step0;
double jd = startJD;
prevDist = findDistance(jd, object1, object2, PhenomenaTypeIndex::Conjuction);
prevDist = findDistance(jd, object1, object2, PhenomenaTypeIndex::Conjunction);
jd += step;
while (jd <= stopJD)
{
dist = findDistance(jd, object1, object2, PhenomenaTypeIndex::Conjuction);
dist = findDistance(jd, object1, object2, PhenomenaTypeIndex::Conjunction);
double factor = qAbs((dist - prevDist) / dist);
if (factor > 10.)
step = step0 * factor / 10.;
Expand All @@ -4609,7 +4629,7 @@ QMap<double, double> AstroCalcDialog::findGreatestElongationApproach(PlanetP& ob
step = step0;
while (jd <= stopJD)
{
dist = findDistance(jd, object1, object2, PhenomenaTypeIndex::Conjuction);
dist = findDistance(jd, object1, object2, PhenomenaTypeIndex::Conjunction);
if (dist<prevDist)
break;

Expand Down Expand Up @@ -4637,19 +4657,19 @@ bool AstroCalcDialog::findPreciseGreatestElongation(QPair<double, double>* out,
if (out == Q_NULLPTR)
return false;

prevDist = findDistance(JD, object1, object2, PhenomenaTypeIndex::Conjuction);
prevDist = findDistance(JD, object1, object2, PhenomenaTypeIndex::Conjunction);
step = -step / 2.;

while (true)
{
JD += step;
dist = findDistance(JD, object1, object2, PhenomenaTypeIndex::Conjuction);
dist = findDistance(JD, object1, object2, PhenomenaTypeIndex::Conjunction);

if (qAbs(step) < 1. / 1440.)
{
out->first = JD - step / 2.0;
out->second = findDistance(JD - step / 2.0, object1, object2, PhenomenaTypeIndex::Conjuction);
if (out->second > findDistance(JD - 5.0, object1, object2, PhenomenaTypeIndex::Conjuction))
out->second = findDistance(JD - step / 2.0, object1, object2, PhenomenaTypeIndex::Conjunction);
if (out->second > findDistance(JD - 5.0, object1, object2, PhenomenaTypeIndex::Conjunction))
{
if (object1->getJ2000EquatorialPos(core).longitude()>object2->getJ2000EquatorialPos(core).longitude())
out->second *= -1.0; // let's use negative value for eastern elongations
Expand Down
4 changes: 2 additions & 2 deletions src/gui/AstroCalcDialog.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,7 @@ class AstroCalcDialog : public StelDialog
};

enum PhenomenaTypeIndex {
Conjuction = 0,
Conjunction = 0,
Opposition = 1,
GreatestElongation = 2,
StationaryPoint = 3,
Expand Down Expand Up @@ -386,7 +386,7 @@ private slots:
//! angular separation ("conjunction" defined as equality of right ascension
//! of two bodies), and current solution is not accurate and slow.
//! @note modes: 0 - conjuction, 1 - opposition, 2 - greatest elongation
QMap<double, double> findClosestApproach(PlanetP& object1, StelObjectP& object2, double startJD, double stopJD, double maxSeparation, int mode);
QMap<double, double> findClosestApproach(PlanetP& object1, StelObjectP& object2, const double startJD, const double stopJD, const double maxSeparation, const int mode);
// TODO: Doc?
double findDistance(double JD, PlanetP object1, StelObjectP object2, int mode);
// TODO: Doc?
Expand Down