Planet.cpp   Planet.cpp 
skipping to change at line 27 skipping to change at line 27
* Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA. * Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA.
*/ */
#include "StelApp.hpp" #include "StelApp.hpp"
#include "StelCore.hpp" #include "StelCore.hpp"
#include "StelFileMgr.hpp" #include "StelFileMgr.hpp"
#include "StelTexture.hpp" #include "StelTexture.hpp"
#include "StelSkyDrawer.hpp" #include "StelSkyDrawer.hpp"
#include "SolarSystem.hpp" #include "SolarSystem.hpp"
#include "Planet.hpp" #include "Planet.hpp"
#include "planetsephems/precession.h"
#include "StelObserver.hpp"
#include "StelProjector.hpp" #include "StelProjector.hpp"
#include "sidereal_time.h" #include "sidereal_time.h"
#include "StelTextureMgr.hpp" #include "StelTextureMgr.hpp"
#include "StelModuleMgr.hpp" #include "StelModuleMgr.hpp"
#include "StarMgr.hpp" #include "StarMgr.hpp"
#include "StelMovementMgr.hpp" #include "StelMovementMgr.hpp"
#include "StelPainter.hpp" #include "StelPainter.hpp"
#include "StelTranslator.hpp" #include "StelTranslator.hpp"
#include "StelUtils.hpp" #include "StelUtils.hpp"
#include "StelOpenGL.hpp" #include "StelOpenGL.hpp"
#include <iomanip> #include <iomanip>
#include <limits>
#include <QTextStream> #include <QTextStream>
#include <QString> #include <QString>
#include <QDebug> #include <QDebug>
#include <QVarLengthArray> #include <QVarLengthArray>
#include <QOpenGLContext> #include <QOpenGLContext>
#include <QOpenGLShader> #include <QOpenGLShader>
Vec3f Planet::labelColor = Vec3f(0.4,0.4,0.8); Vec3f Planet::labelColor = Vec3f(0.4,0.4,0.8);
Vec3f Planet::orbitColor = Vec3f(1,0.6,1); Vec3f Planet::orbitColor = Vec3f(1,0.6,1);
StelTextureSP Planet::hintCircleTex; StelTextureSP Planet::hintCircleTex;
skipping to change at line 66 skipping to change at line 68
QOpenGLShaderProgram* Planet::moonShaderProgram=NULL; QOpenGLShaderProgram* Planet::moonShaderProgram=NULL;
Planet::MoonShaderVars Planet::moonShaderVars; Planet::MoonShaderVars Planet::moonShaderVars;
QMap<Planet::PlanetType, QString> Planet::pTypeMap; QMap<Planet::PlanetType, QString> Planet::pTypeMap;
QMap<Planet::ApparentMagnitudeAlgorithm, QString> Planet::vMagAlgorithmMap; QMap<Planet::ApparentMagnitudeAlgorithm, QString> Planet::vMagAlgorithmMap;
Planet::Planet(const QString& englishName, Planet::Planet(const QString& englishName,
int flagLighting, int flagLighting,
double radius, double radius,
double oblateness, double oblateness,
Vec3f color, Vec3f halocolor,
float albedo, float albedo,
const QString& atexMapName, const QString& atexMapName,
const QString& anormalMapName, const QString& anormalMapName,
posFuncType coordFunc, posFuncType coordFunc,
void* auserDataPtr, void* auserDataPtr,
OsculatingFunctType *osculatingFunc, OsculatingFunctType *osculatingFunc,
bool acloseOrbit, bool acloseOrbit,
bool hidden, bool hidden,
bool hasAtmosphere, bool hasAtmosphere,
bool hasHalo, bool hasHalo,
const QString& pTypeStr) const QString& pTypeStr)
: englishName(englishName), : englishName(englishName),
flagLighting(flagLighting), flagLighting(flagLighting),
radius(radius), radius(radius),
oneMinusOblateness(1.0-oblateness), oneMinusOblateness(1.0-oblateness),
color(color), haloColor(halocolor),
albedo(albedo), albedo(albedo),
axisRotation(0.), axisRotation(0.),
rings(NULL), rings(NULL),
sphereScale(1.f), sphereScale(1.f),
lastJD(J2000), lastJDE(J2000),
coordFunc(coordFunc), coordFunc(coordFunc),
userDataPtr(auserDataPtr), userDataPtr(auserDataPtr),
osculatingFunc(osculatingFunc), osculatingFunc(osculatingFunc),
parent(NULL), parent(NULL),
hidden(hidden), hidden(hidden),
atmosphere(hasAtmosphere), atmosphere(hasAtmosphere),
halo(hasHalo) halo(hasHalo)
{ {
texMapName = atexMapName; texMapName = atexMapName;
normalMapName = anormalMapName; normalMapName = anormalMapName;
lastOrbitJD =0; lastOrbitJDE =0;
deltaJD = StelCore::JD_SECOND; deltaJDE = StelCore::JD_SECOND;
orbitCached = 0; orbitCached = 0;
closeOrbit = acloseOrbit; closeOrbit = acloseOrbit;
deltaOrbitJD = 0; deltaOrbitJDE = 0;
distance = 0; distance = 0;
// Initialize pType with the key found in pTypeMap, or mark planet t ype as undefined. // Initialize pType with the key found in pTypeMap, or mark planet t ype as undefined.
// The latter condition should obviously never happen. // The latter condition should obviously never happen.
pType = pTypeMap.key(pTypeStr, Planet::isUNDEFINED); pType = pTypeMap.key(pTypeStr, Planet::isUNDEFINED);
vMagAlgorithm = Planet::UndefinedAlgorithm; vMagAlgorithm = Planet::UndefinedAlgorithm;
eclipticPos = Vec3d(0.,0.,0.); eclipticPos = Vec3d(0.,0.,0.);
rotLocalToParent = Mat4d::identity(); rotLocalToParent = Mat4d::identity();
texMap = StelApp::getInstance().getTextureManager().createTextureThr ead(StelFileMgr::getInstallationDir()+"/textures/"+texMapName, StelTexture: :StelTextureParams(true, GL_LINEAR, GL_REPEAT)); texMap = StelApp::getInstance().getTextureManager().createTextureThr ead(StelFileMgr::getInstallationDir()+"/textures/"+texMapName, StelTexture: :StelTextureParams(true, GL_LINEAR, GL_REPEAT));
normalMap = StelApp::getInstance().getTextureManager().createTexture Thread(StelFileMgr::getInstallationDir()+"/textures/"+normalMapName, StelTe xture::StelTextureParams(true, GL_LINEAR, GL_REPEAT)); normalMap = StelApp::getInstance().getTextureManager().createTexture Thread(StelFileMgr::getInstallationDir()+"/textures/"+normalMapName, StelTe xture::StelTextureParams(true, GL_LINEAR, GL_REPEAT));
nameI18 = englishName; nameI18 = englishName;
nativeName = ""; nativeName = "";
if (englishName!="Pluto") if (englishName!="Pluto")
{ {
deltaJD = 0.001*StelCore::JD_SECOND; deltaJDE = 0.001*StelCore::JD_SECOND;
} }
flagLabels = true; flagLabels = true;
flagNativeName = true; flagNativeName = true;
flagTranslatedName = true; flagTranslatedName = true;
} }
// called in SolarSystem::init() before first planet is created. Loads pTyp eMap. // called in SolarSystem::init() before first planet is created. Loads pTyp eMap.
void Planet::init() void Planet::init()
{ {
if (pTypeMap.count() > 0 ) if (pTypeMap.count() > 0 )
skipping to change at line 198 skipping to change at line 200
QString Planet::getNameI18n() const QString Planet::getNameI18n() const
{ {
return nameI18; return nameI18;
} }
// Return the information string "ready to print" :) // Return the information string "ready to print" :)
QString Planet::getInfoString(const StelCore* core, const InfoStringGroup& flags) const QString Planet::getInfoString(const StelCore* core, const InfoStringGroup& flags) const
{ {
QString str; QString str;
QTextStream oss(&str); QTextStream oss(&str);
double az_app, alt_app;
StelUtils::rectToSphe(&az_app,&alt_app,getAltAzPosApparent(core));
Q_UNUSED(az_app);
if (flags&Name) if (flags&Name)
{ {
oss << "<h2>" << getNameI18n(); // UI translation can diffe r from sky translation oss << "<h2>" << getNameI18n(); // UI translation can diffe r from sky translation
oss.setRealNumberNotation(QTextStream::FixedNotation); oss.setRealNumberNotation(QTextStream::FixedNotation);
oss.setRealNumberPrecision(1); oss.setRealNumberPrecision(1);
if (sphereScale != 1.f) if (sphereScale != 1.f)
oss << QString::fromUtf8(" (\xC3\x97") << sphereScal e << ")"; oss << QString::fromUtf8(" (\xC3\x97") << sphereScal e << ")";
oss << "</h2>"; oss << "</h2>";
} }
if (flags&ObjectType && getPlanetType()!=isUNDEFINED) if (flags&ObjectType && getPlanetType()!=isUNDEFINED)
{ {
oss << q_("Type: <b>%1</b>").arg(q_(getPlanetTypeString())) << "<br />"; oss << q_("Type: <b>%1</b>").arg(q_(getPlanetTypeString())) << "<br />";
} }
if (flags&Magnitude) if (flags&Magnitude && getVMagnitude(core)!=std::numeric_limits<floa t>::infinity())
{ {
if (core->getSkyDrawer()->getFlagHasAtmosphere()) if (core->getSkyDrawer()->getFlagHasAtmosphere() && (alt_app
oss << q_("Magnitude: <b>%1</b> (extincted to: <b>%2</b> >-3.0*M_PI/180.0)) // Don't show extincted magnitude much below horizon whe
)").arg(QString::number(getVMagnitude(core), 'f', 2), re model is meaningless.
oss << q_("Magnitude: <b>%1</b> (extincted to: <b>%2
QString::number(getVMagnitudeWithExtinction(core), 'f', 2)) << "<br </b>)").arg(QString::number(getVMagnitude(core), 'f', 2),
>";
QString::number(getVMagnitudeWithExtinction(core), 'f', 2)) <<
"<br>";
else else
oss << q_("Magnitude: <b>%1</b>").arg(getVMagnitude(core ), 0, 'f', 2) << "<br>"; oss << q_("Magnitude: <b>%1</b>").arg(getVMagnitude( core), 0, 'f', 2) << "<br>";
} }
if (flags&AbsoluteMagnitude) if (flags&AbsoluteMagnitude && getVMagnitude(core)!=std::numeric_lim its<float>::infinity())
oss << q_("Absolute Magnitude: %1").arg(getVMagnitude(core)- 5.*(std::log10(getJ2000EquatorialPos(core).length()*AU/PARSEC)-1.), 0, 'f', 2) << "<br>"; oss << q_("Absolute Magnitude: %1").arg(getVMagnitude(core)- 5.*(std::log10(getJ2000EquatorialPos(core).length()*AU/PARSEC)-1.), 0, 'f', 2) << "<br>";
oss << getPositionInfoString(core, flags); oss << getPositionInfoString(core, flags);
if (flags&Extra) // GZ This is mostly for debugging. Maybe also useful for letting pe
{ ople use our results to cross-check theirs, but we should not act as refere
static SolarSystem *ssystem=GETSTELMODULE(SolarSystem); nce, currently...
double ecl= ssystem->getEarth()->getRotObliquity(core->getJD // TODO: maybe separate this out into:
ay()); //if (flags&EclipticCoordXYZ)
if (core->getCurrentLocation().planetName=="Earth") // For now: add to EclipticCoord
oss << q_("Obliquity (of date, for Earth): %1").arg( //if (flags&EclipticCoord)
StelUtils::radToDmsStr(ecl, true)) << "<br>"; // oss << q_("Ecliptical XYZ (VSOP87A): %1/%2/%3").arg(QString:
//if (englishName!="Sun") :number(eclipticPos[0], 'f', 3), QString::number(eclipticPos[1], 'f', 3), Q
// oss << q_("Obliquity (of date): %1").arg(StelUtils:: String::number(eclipticPos[2], 'f', 3)) << "<br>";
radToDmsStr(getRotObliquity(core->getJDay()), true)) << "<br>"; // if (flags&Extra)
} // {
// static SolarSystem *ssystem=GETSTELMODULE(SolarSystem);
// double ecl= ssystem->getEarth()->getRotObliquity(core->getJD
ay());
// if (core->getCurrentLocation().planetName=="Earth")
// oss << q_("Ecliptic obliquity (of date): %1").arg(St
elUtils::radToDmsStr(ecl, true)) << "<br>";
// //if (englishName!="Sun")
// // oss << q_("Obliquity (of date): %1").arg(StelUtils::
radToDmsStr(getRotObliquity(core->getJDay()), true)) << "<br>";
// }
if (flags&Distance) if (flags&Distance)
{ {
double distanceAu = getJ2000EquatorialPos(core).length(); double distanceAu = getJ2000EquatorialPos(core).length();
double distanceKm = AU * distanceAu; double distanceKm = AU * distanceAu;
if (distanceAu < 0.1) if (distanceAu < 0.1)
{ {
// xgettext:no-c-format // xgettext:no-c-format
oss << QString(q_("Distance: %1AU (%2 km)")) oss << QString(q_("Distance: %1AU (%2 km)"))
.arg(distanceAu, 0, 'f', 6) .arg(distanceAu, 0, 'f', 6)
skipping to change at line 279 skipping to change at line 290
{ {
oss << q_("Apparent diameter: %1").arg(StelUtils::ra dToDmsStr(angularSize, true)); oss << q_("Apparent diameter: %1").arg(StelUtils::ra dToDmsStr(angularSize, true));
} }
oss << "<br>"; oss << "<br>";
} }
double siderealPeriod = getSiderealPeriod(); double siderealPeriod = getSiderealPeriod();
double siderealDay = getSiderealDay(); double siderealDay = getSiderealDay();
if (flags&Extra) if (flags&Extra)
{ {
// This is a string you can activate for debugging. It shows
the distance between observer and center of the body you are standing on.
// May be helpful for debugging critical parallax correction
s for eclipses.
// For general use, find a better location first.
// oss << q_("Planetocentric distance &rho;: %1 (km)").arg(c
ore->getCurrentObserver()->getDistanceFromCenter() * AU) <<"<br>";
if (siderealPeriod>0) if (siderealPeriod>0)
{ {
// TRANSLATORS: Sidereal (orbital) period for solar system bodies in days and in Julian years (symbol: a) // TRANSLATORS: Sidereal (orbital) period for solar system bodies in days and in Julian years (symbol: a)
oss << q_("Sidereal period: %1 days (%2 a)").arg(QSt ring::number(siderealPeriod, 'f', 2)).arg(QString::number(siderealPeriod/36 5.25, 'f', 3)) << "<br>"; oss << q_("Sidereal period: %1 days (%2 a)").arg(QSt ring::number(siderealPeriod, 'f', 2)).arg(QString::number(siderealPeriod/36 5.25, 'f', 3)) << "<br>";
if (qAbs(siderealDay)>0) if (qAbs(siderealDay)>0)
{ {
oss << q_("Sidereal day: %1").arg(StelUtils: :hoursToHmsStr(qAbs(siderealDay*24))) << "<br>"; oss << q_("Sidereal day: %1").arg(StelUtils: :hoursToHmsStr(qAbs(siderealDay*24))) << "<br>";
oss << q_("Mean solar day: %1").arg(StelUtil s::hoursToHmsStr(qAbs(getMeanSolarDay()*24))) << "<br>"; oss << q_("Mean solar day: %1").arg(StelUtil s::hoursToHmsStr(qAbs(getMeanSolarDay()*24))) << "<br>";
} }
else if (re.period==0.) else if (re.period==0.)
{ {
oss << q_("The period of rotation is chaotic ") << "<br>"; oss << q_("The period of rotation is chaotic ") << "<br>";
} }
} }
if (englishName.compare("Sun")!=0) if (englishName != "Sun")
{ {
const Vec3d& observerHelioPos = core->getObserverHel iocentricEclipticPos(); const Vec3d& observerHelioPos = core->getObserverHel iocentricEclipticPos();
const double elongation = getElongation(observerHeli
oPos);
QString phase = "";
double deg;
bool sign;
StelUtils::radToDecDeg(elongation, sign, deg);
if (deg>=357.5 && deg<=2.5)
phase = q_("New Moon");
if (deg>=87.5 && deg<=92.5)
phase = q_("Last Quarter");
if (deg>=177.5 && deg<=182.5)
phase = q_("Full Moon");
if (deg>=267.5 && deg<=272.5)
phase = q_("First Quarter");
oss << QString(q_("Phase Angle: %1")).arg(StelUtils: :radToDmsStr(getPhaseAngle(observerHelioPos))) << "<br>"; oss << QString(q_("Phase Angle: %1")).arg(StelUtils: :radToDmsStr(getPhaseAngle(observerHelioPos))) << "<br>";
oss << QString(q_("Elongation: %1")).arg(StelUtils:: oss << QString(q_("Elongation: %1")).arg(StelUtils::
radToDmsStr(getElongation(observerHelioPos))) << "<br>"; radToDmsStr(elongation)) << "<br>";
oss << QString(q_("Phase: %1")).arg(getPhase(observe
rHelioPos), 0, 'f', 2) << "<br>"; if (englishName=="Moon" && !phase.isEmpty())
{
if (qRound(deg)==180.f || qRound(deg)==90.f
|| qRound(deg)==270.f || qRound(deg)==0.f)
phase = QString("%1 (<b>%2</b>)").ar
g(getPhase(observerHelioPos), 0, 'f', 2).arg(phase);
else
phase = QString("%1 (%2)").arg(getPh
ase(observerHelioPos), 0, 'f', 2).arg(phase);
}
else
phase = QString("%1").arg(getPhase(observerH
elioPos), 0, 'f', 2);
oss << QString(q_("Phase: %1")).arg(phase) << "<br>"
;
oss << QString(q_("Illuminated: %1%")).arg(getPhase( observerHelioPos) * 100, 0, 'f', 1) << "<br>"; oss << QString(q_("Illuminated: %1%")).arg(getPhase( observerHelioPos) * 100, 0, 'f', 1) << "<br>";
} }
if (englishName=="Sun")
{
// Only show during eclipse, show percent?
static SolarSystem *ssystem=GETSTELMODULE(SolarSyste
m);
// Debug solution:
// float eclipseFactor = ssystem->getEclipseFactor(core
);
// oss << QString(q_("Eclipse Factor: %1 alpha: %2")).a
rg(eclipseFactor).arg(-0.1f*qMax(-10.0f, (float) std::log10(eclipseFactor))
) << "<br>";
// Release version:
float eclipseFactor = 100.f*(1.f-ssystem->getEclipse
Factor(core));
if (eclipseFactor>0.f)
oss << QString(q_("Eclipse Factor: %1%")).ar
g(eclipseFactor) << "<br>";
}
} }
postProcessInfoString(str, flags); postProcessInfoString(str, flags);
return str; return str;
} }
//! Get sky label (sky translation) //! Get sky label (sky translation)
QString Planet::getSkyLabel(const StelCore*) const QString Planet::getSkyLabel(const StelCore*) const
{ {
skipping to change at line 373 skipping to change at line 426
void Planet::setRotationElements(float _period, float _offset, double _epoc h, float _obliquity, float _ascendingNode, float _precessionRate, double _s iderealPeriod ) void Planet::setRotationElements(float _period, float _offset, double _epoc h, float _obliquity, float _ascendingNode, float _precessionRate, double _s iderealPeriod )
{ {
re.period = _period; re.period = _period;
re.offset = _offset; re.offset = _offset;
re.epoch = _epoch; re.epoch = _epoch;
re.obliquity = _obliquity; re.obliquity = _obliquity;
re.ascendingNode = _ascendingNode; re.ascendingNode = _ascendingNode;
re.precessionRate = _precessionRate; re.precessionRate = _precessionRate;
re.siderealPeriod = _siderealPeriod; // used for drawing orbit line s re.siderealPeriod = _siderealPeriod; // used for drawing orbit line s
deltaOrbitJD = re.siderealPeriod/ORBIT_SEGMENTS; deltaOrbitJDE = re.siderealPeriod/ORBIT_SEGMENTS;
} }
Vec3d Planet::getJ2000EquatorialPos(const StelCore *core) const Vec3d Planet::getJ2000EquatorialPos(const StelCore *core) const
{ {
return StelCore::matVsop87ToJ2000.multiplyWithoutTranslation(getHeli ocentricEclipticPos() - core->getObserverHeliocentricEclipticPos()); return StelCore::matVsop87ToJ2000.multiplyWithoutTranslation(getHeli ocentricEclipticPos() - core->getObserverHeliocentricEclipticPos());
} }
// Compute the position in the parent Planet coordinate system // Compute the position in the parent Planet coordinate system
// Actually call the provided function to compute the ecliptical position // Actually call the provided function to compute the ecliptical position
void Planet::computePositionWithoutOrbits(const double dateJD) void Planet::computePositionWithoutOrbits(const double dateJDE)
{ {
if (fabs(lastJD-dateJD)>deltaJD) if (fabs(lastJDE-dateJDE)>deltaJDE)
{ {
coordFunc(dateJD, eclipticPos, userDataPtr); coordFunc(dateJDE, eclipticPos, userDataPtr);
lastJD = dateJD; lastJDE = dateJDE;
} }
} }
double Planet::getRotObliquity(double JDay) const // return value in radians!
// For Earth, this is epsilon_A, the angle between earth's rotational axis
and mean ecliptic of date.
// Details: e.g. Hilton etal, Report on Precession and the Ecliptic, Cel.Me
ch.Dyn.Astr.94:351-67 (2006), Fig1.
double Planet::getRotObliquity(double JDE) const
{ {
// JDay=2451545.0 for J2000.0 // JDay=2451545.0 for J2000.0
if (englishName=="Earth") if (englishName=="Earth")
return get_mean_ecliptical_obliquity(JDay) *M_PI/180.0; return getPrecessionAngleVondrakEpsilon(JDE);
else else
return re.obliquity; return re.obliquity;
} }
bool willCastShadow(const Planet* thisPlanet, const Planet* p) bool willCastShadow(const Planet* thisPlanet, const Planet* p)
{ {
Vec3d thisPos = thisPlanet->getHeliocentricEclipticPos(); Vec3d thisPos = thisPlanet->getHeliocentricEclipticPos();
Vec3d planetPos = p->getHeliocentricEclipticPos(); Vec3d planetPos = p->getHeliocentricEclipticPos();
// If the planet p is farther from the sun than this planet, it can' t cast shadow on it. // If the planet p is farther from the sun than this planet, it can' t cast shadow on it.
skipping to change at line 444 skipping to change at line 500
{ {
if (willCastShadow(this, planet.data())) if (willCastShadow(this, planet.data()))
res.append(planet.data()); res.append(planet.data());
} }
if (willCastShadow(this, parent.data())) if (willCastShadow(this, parent.data()))
res.append(parent.data()); res.append(parent.data());
return res; return res;
} }
void Planet::computePosition(const double dateJD) void Planet::computePosition(const double dateJDE)
{ {
// Make sure the parent position is computed for the dateJD, otherwi se // Make sure the parent position is computed for the dateJD, otherwi se
// getHeliocentricPos() would return incorect values. // getHeliocentricPos() would return incorect values.
if (parent) if (parent)
parent->computePositionWithoutOrbits(dateJD); parent->computePositionWithoutOrbits(dateJDE);
if (orbitFader.getInterstate()>0.000001 && deltaOrbitJD > 0 && (fabs (lastOrbitJD-dateJD)>deltaOrbitJD || !orbitCached)) if (orbitFader.getInterstate()>0.000001 && deltaOrbitJDE > 0 && (fab s(lastOrbitJDE-dateJDE)>deltaOrbitJDE || !orbitCached))
{ {
StelCore *core=StelApp::getInstance().getCore();
double calc_date; double calc_date;
// int delta_points = (int)(0.5 + (date - lastOrbitJD)/date_ increment); // int delta_points = (int)(0.5 + (date - lastOrbitJD)/date_ increment);
int delta_points; int delta_points;
if( dateJD > lastOrbitJD ) if( dateJDE > lastOrbitJDE )
{ {
delta_points = (int)(0.5 + (dateJD - lastOrbitJD)/de ltaOrbitJD); delta_points = (int)(0.5 + (dateJDE - lastOrbitJDE)/ deltaOrbitJDE);
} }
else else
{ {
delta_points = (int)(-0.5 + (dateJD - lastOrbitJD)/d eltaOrbitJD); delta_points = (int)(-0.5 + (dateJDE - lastOrbitJDE) /deltaOrbitJDE);
} }
double new_date = lastOrbitJD + delta_points*deltaOrbitJD; double new_date = lastOrbitJDE + delta_points*deltaOrbitJDE;
// qDebug( "Updating orbit coordinates for %s (delta %f) (%d points)\n", getEnglishName().toUtf8().data(), deltaOrbitJD, delta_points); // qDebug( "Updating orbit coordinates for %s (delta %f) (%d points)\n", getEnglishName().toUtf8().data(), deltaOrbitJD, delta_points);
if( delta_points > 0 && delta_points < ORBIT_SEGMENTS && orb itCached) if( delta_points > 0 && delta_points < ORBIT_SEGMENTS && orb itCached)
{ {
for( int d=0; d<ORBIT_SEGMENTS; d++ ) for( int d=0; d<ORBIT_SEGMENTS; d++ )
{ {
if(d + delta_points >= ORBIT_SEGMENTS ) if(d + delta_points >= ORBIT_SEGMENTS )
{ {
// calculate new points // calculate new points
calc_date = new_date + (d-ORBIT_SEGM ENTS/2)*deltaOrbitJD; calc_date = new_date + (d-ORBIT_SEGM ENTS/2)*deltaOrbitJDE;
// date increments between points wi ll not be completely constant though // date increments between points wi ll not be completely constant though
computeTransMatrix(calc_date); computeTransMatrix(calc_date-core->c omputeDeltaT(calc_date)/86400.0, calc_date);
if (osculatingFunc) if (osculatingFunc)
{ {
(*osculatingFunc)(dateJD,cal c_date,eclipticPos); (*osculatingFunc)(dateJDE,ca lc_date,eclipticPos);
} }
else else
{ {
coordFunc(calc_date, eclipti cPos, userDataPtr); coordFunc(calc_date, eclipti cPos, userDataPtr);
} }
orbitP[d] = eclipticPos; orbitP[d] = eclipticPos;
orbit[d] = getHeliocentricEclipticPo s(); orbit[d] = getHeliocentricEclipticPo s();
} }
else else
{ {
orbitP[d] = orbitP[d+delta_points]; orbitP[d] = orbitP[d+delta_points];
orbit[d] = getHeliocentricPos(orbitP [d]); orbit[d] = getHeliocentricPos(orbitP [d]);
} }
} }
lastOrbitJD = new_date; lastOrbitJDE = new_date;
} }
else if( delta_points < 0 && abs(delta_points) < ORBIT_SEGME NTS && orbitCached) else if( delta_points < 0 && abs(delta_points) < ORBIT_SEGME NTS && orbitCached)
{ {
for( int d=ORBIT_SEGMENTS-1; d>=0; d-- ) for( int d=ORBIT_SEGMENTS-1; d>=0; d-- )
{ {
if(d + delta_points < 0 ) if(d + delta_points < 0 )
{ {
// calculate new points // calculate new points
calc_date = new_date + (d-ORBIT_SEGM ENTS/2)*deltaOrbitJD; calc_date = new_date + (d-ORBIT_SEGM ENTS/2)*deltaOrbitJDE;
computeTransMatrix(calc_date); computeTransMatrix(calc_date-core->c omputeDeltaT(calc_date)/86400.0, calc_date);
if (osculatingFunc) { if (osculatingFunc) {
(*osculatingFunc)(dateJD,cal c_date,eclipticPos); (*osculatingFunc)(dateJDE,ca lc_date,eclipticPos);
} }
else else
{ {
coordFunc(calc_date, eclipti cPos, userDataPtr); coordFunc(calc_date, eclipti cPos, userDataPtr);
} }
orbitP[d] = eclipticPos; orbitP[d] = eclipticPos;
orbit[d] = getHeliocentricEclipticPo s(); orbit[d] = getHeliocentricEclipticPo s();
} }
else else
{ {
orbitP[d] = orbitP[d+delta_points]; orbitP[d] = orbitP[d+delta_points];
orbit[d] = getHeliocentricPos(orbitP [d]); orbit[d] = getHeliocentricPos(orbitP [d]);
} }
} }
lastOrbitJD = new_date; lastOrbitJDE = new_date;
} }
else if( delta_points || !orbitCached) else if( delta_points || !orbitCached)
{ {
// update all points (less efficient) // update all points (less efficient)
for( int d=0; d<ORBIT_SEGMENTS; d++ ) for( int d=0; d<ORBIT_SEGMENTS; d++ )
{ {
calc_date = dateJD + (d-ORBIT_SEGMENTS/2)*de calc_date = dateJDE + (d-ORBIT_SEGMENTS/2)*d
ltaOrbitJD; eltaOrbitJDE;
computeTransMatrix(calc_date); computeTransMatrix(calc_date-core->computeDe
ltaT(calc_date)/86400.0, calc_date);
if (osculatingFunc) if (osculatingFunc)
{ {
(*osculatingFunc)(dateJD,calc_date,e clipticPos); (*osculatingFunc)(dateJDE,calc_date, eclipticPos);
} }
else else
{ {
coordFunc(calc_date, eclipticPos, us erDataPtr); coordFunc(calc_date, eclipticPos, us erDataPtr);
} }
orbitP[d] = eclipticPos; orbitP[d] = eclipticPos;
orbit[d] = getHeliocentricEclipticPos(); orbit[d] = getHeliocentricEclipticPos();
} }
lastOrbitJD = dateJD; lastOrbitJDE = dateJDE;
if (!osculatingFunc) orbitCached = 1; if (!osculatingFunc) orbitCached = 1;
} }
// calculate actual Planet position // calculate actual Planet position
coordFunc(dateJD, eclipticPos, userDataPtr); coordFunc(dateJDE, eclipticPos, userDataPtr);
lastJD = dateJD; lastJDE = dateJDE;
} }
else if (fabs(lastJD-dateJD)>deltaJD) else if (fabs(lastJDE-dateJDE)>deltaJDE)
{ {
// calculate actual Planet position // calculate actual Planet position
coordFunc(dateJD, eclipticPos, userDataPtr); coordFunc(dateJDE, eclipticPos, userDataPtr);
// XXX: do we need to do that even when the orbit is not vis ible? // XXX: do we need to do that even when the orbit is not vis ible?
for( int d=0; d<ORBIT_SEGMENTS; d++ ) for( int d=0; d<ORBIT_SEGMENTS; d++ )
orbit[d]=getHeliocentricPos(orbitP[d]); orbit[d]=getHeliocentricPos(orbitP[d]);
lastJD = dateJD; lastJDE = dateJDE;
} }
} }
// Compute the transformation matrix from the local Planet coordinate to th // Compute the transformation matrix from the local Planet coordinate syste
e parent Planet coordinate m to the parent Planet coordinate system.
void Planet::computeTransMatrix(double jd) // In case of the planets, this makes the axis point to their respective ce
lestial poles.
// TODO: Verify for the other planets if their axes are relative to J2000 e
cliptic (VSOP87A XY plane) or relative to (precessed) ecliptic of date?
void Planet::computeTransMatrix(double JD, double JDE)
{ {
axisRotation = getSiderealTime(jd); // We have to call with both to correct this for earth with the new
model.
axisRotation = getSiderealTime(JD, JDE);
// Special case - heliocentric coordinates are on ecliptic, // Special case - heliocentric coordinates are relative to eclipticJ 2000 (VSOP87A XY plane),
// not solar equator... // not solar equator...
if (parent) if (parent)
{ {
rotLocalToParent = Mat4d::zrotation(re.ascendingNode - re.pr // We can inject a proper precession plus even nutation matr
ecessionRate*(jd-re.epoch)) * Mat4d::xrotation(re.obliquity); ix in this stage, if available.
if (englishName=="Earth")
{
// rotLocalToParent = Mat4d::zrotation(re.ascendingN
ode - re.precessionRate*(jd-re.epoch)) * Mat4d::xrotation(-getRotObliquity(
jd));
// We follow Capitaine's (2003) formulation P=Rz(Chi
_A)*Rx(-omega_A)*Rz(-psi_A)*Rx(eps_o).
// ADS: 2011A&A...534A..22V = A&A 534, A22 (2011): V
ondrak, Capitane, Wallace: New Precession Expressions, valid for long time
intervals:
// See also Hilton et al., Report on Precession and
the Ecliptic. Cel.Mech.Dyn.Astr. 94:351-367 (2006), eqn (6) and (21).
double eps_A, chi_A, omega_A, psi_A;
getPrecessionAnglesVondrak(JDE, &eps_A, &chi_A, &ome
ga_A, &psi_A);
// Canonical precession rotations: Nodal rotation ps
i_A,
// then rotation by omega_A, the angle between EclPo
leJ2000 and EarthPoleOfDate.
// The final rotation by chi_A rotates the equinox (
zero degree).
// To achieve ecliptical coords of date, you just ha
ve now to add a rotX by epsilon_A (obliquity of date).
rotLocalToParent= Mat4d::zrotation(-psi_A) * Mat4d::
xrotation(-omega_A) * Mat4d::zrotation(chi_A);
// Plus nutation IAU-2000B:
if (StelApp::getInstance().getCore()->getUseNutation
())
{
double deltaEps, deltaPsi;
getNutationAngles(JDE, &deltaPsi, &deltaEps)
;
//qDebug() << "deltaEps, arcsec" << deltaEps
*180./M_PI*3600. << "deltaPsi" << deltaPsi*180./M_PI*3600.;
Mat4d nut2000B=Mat4d::xrotation(eps_A) * Mat
4d::zrotation(deltaPsi)* Mat4d::xrotation(-eps_A-deltaEps);
rotLocalToParent=rotLocalToParent*nut2000B;
}
}
else
rotLocalToParent = Mat4d::zrotation(re.ascendingNode
- re.precessionRate*(JDE-re.epoch)) * Mat4d::xrotation(re.obliquity);
} }
} }
Mat4d Planet::getRotEquatorialToVsop87(void) const Mat4d Planet::getRotEquatorialToVsop87(void) const
{ {
Mat4d rval = rotLocalToParent; Mat4d rval = rotLocalToParent;
if (parent) if (parent)
{ {
for (PlanetP p=parent;p->parent;p=p->parent) for (PlanetP p=parent;p->parent;p=p->parent)
rval = p->rotLocalToParent * rval; rval = p->rotLocalToParent * rval;
skipping to change at line 609 skipping to change at line 697
{ {
Mat4d a = Mat4d::identity(); Mat4d a = Mat4d::identity();
if (parent) if (parent)
{ {
for (PlanetP p=parent;p->parent;p=p->parent) for (PlanetP p=parent;p->parent;p=p->parent)
a = p->rotLocalToParent * a; a = p->rotLocalToParent * a;
} }
rotLocalToParent = a.transpose() * m; rotLocalToParent = a.transpose() * m;
} }
// Compute the z rotation to use from equatorial to geographic coordinates // Compute the z rotation to use from equatorial to geographic coordinates.
double Planet::getSiderealTime(double jd) const // We need both JD and JDE here for Earth. (For other planets only JDE.)
double Planet::getSiderealTime(double JD, double JDE) const
{ {
if (englishName=="Earth") if (englishName=="Earth")
{ { // Check to make sure that nutation is just those few arcsec
return get_apparent_sidereal_time(jd); onds.
if (StelApp::getInstance().getCore()->getUseNutation())
return get_apparent_sidereal_time(JD, JDE);
else
return get_mean_sidereal_time(JD, JDE);
} }
double t = jd - re.epoch; double t = JDE - re.epoch;
// oops... avoid division by zero (typical case for moons with chaot ic period of rotation) // oops... avoid division by zero (typical case for moons with chaot ic period of rotation)
double rotations = 1.f; // NOTE: Maybe 1e-3 will be better? double rotations = 1.f; // NOTE: Maybe 1e-3 will be better?
if (re.period!=0.) // OK, it's not a moon with chaotic period of rot ation :) if (re.period!=0.) // OK, it's not a moon with chaotic period of rot ation :)
{ {
rotations = t / (double) re.period; rotations = t / (double) re.period;
} }
double wholeRotations = floor(rotations); double wholeRotations = floor(rotations);
double remainder = rotations - wholeRotations; double remainder = rotations - wholeRotations;
if (englishName=="Jupiter") if (englishName=="Jupiter")
{ {
if( re.offset >= 0.0 ) if( re.offset >= 0.0 )
{ {
// use semi-empirical coefficient for GRS drift // use semi-empirical coefficient for GRS drift
// qDebug() << "Jupiter: offset = " << re.offset << " --> rotation = " << (remainder * 360. + re.offset - 0.2483 * qAbs(jd - 24 56172)); // qDebug() << "Jupiter: offset = " << re.offset << " --> rotation = " << (remainder * 360. + re.offset - 0.2483 * qAbs(jd - 24 56172));
return remainder * 360. + re.offset - 0.2483 * qAbs( jd - 2456172); return remainder * 360. + re.offset - 0.2483 * qAbs( JDE - 2456172);
} }
else else
{ {
// http://www.projectpluto.com/grs_form.htm // http://www.projectpluto.com/grs_form.htm
// CM( System II) = 181.62 + 870.1869147 * jd + cor rection [870d rotation every day] // CM( System II) = 181.62 + 870.1869147 * jd + cor rection [870d rotation every day]
const double rad = M_PI/180.; const double rad = M_PI/180.;
double jup_mean = (jd - 2455636.938) * 360. / 4332.8 9709; double jup_mean = (JDE - 2455636.938) * 360. / 4332. 89709;
double eqn_center = 5.55 * sin( rad*jup_mean); double eqn_center = 5.55 * sin( rad*jup_mean);
double angle = (jd - 2451870.628) * 360. / 398.884 - eqn_center; double angle = (JDE - 2451870.628) * 360. / 398.884 - eqn_center;
//double correction = 11 * sin( rad*angle) + 5 * cos ( rad*angle)- 1.25 * cos( rad*jup_mean) - eqn_center; // original correctio n //double correction = 11 * sin( rad*angle) + 5 * cos ( rad*angle)- 1.25 * cos( rad*jup_mean) - eqn_center; // original correctio n
double correction = 25.8 + 11 * sin( rad*angle) - 2. 5 * cos( rad*jup_mean) - eqn_center; // light speed correction not used bec ause in stellarium the jd is manipulated for that double correction = 25.8 + 11 * sin( rad*angle) - 2. 5 * cos( rad*jup_mean) - eqn_center; // light speed correction not used bec ause in stellarium the jd is manipulated for that
double cm2=181.62 + 870.1869147 * jd + correction; double cm2=181.62 + 870.1869147 * JDE + correction; // Central Meridian II
cm2=cm2 - 360.0*(int)(cm2/360.); cm2=cm2 - 360.0*(int)(cm2/360.);
// http://www.skyandtelescope.com/observing/transit- times-of-jupiters-great-red-spot/ writes: // http://www.skyandtelescope.com/observing/transit- times-of-jupiters-great-red-spot/ writes:
// The predictions assume the Red Spot was at Jovian System II longitude 216° in September 2014 and continues to drift 1.25° pe r month, based on historical trends noted by JUPOS. // The predictions assume the Red Spot was at Jovian System II longitude 216° in September 2014 and continues to drift 1.25° pe r month, based on historical trends noted by JUPOS.
// GRS longitude was at 2014-09-08 216d with a drift of 1.25d every month // GRS longitude was at 2014-09-08 216d with a drift of 1.25d every month
double longitudeGRS=216+1.25*( jd - 2456908)/30; double longitudeGRS=216+1.25*( JDE - 2456908)/30;
// qDebug() << "Jupiter: CM2 = " << cm2 << " longitu deGRS = " << longitudeGRS << " --> rotation = " << (cm2 - longitudeGRS); // qDebug() << "Jupiter: CM2 = " << cm2 << " longitu deGRS = " << longitudeGRS << " --> rotation = " << (cm2 - longitudeGRS);
return cm2 - longitudeGRS + 25.; // + 25 = Jupiter T exture not 0d return cm2 - longitudeGRS + 25.; // + 25 = Jupiter T exture not 0d
// To verify: // To verify:
// GRS at 2015-02-26 23:07 UT on picture at https:// maximusphotography.files.wordpress.com/2015/03/jupiter-febr-26-2015.jpg // GRS at 2015-02-26 23:07 UT on picture at https:// maximusphotography.files.wordpress.com/2015/03/jupiter-febr-26-2015.jpg
// 2014-02-25 19:03 UT http://www.damianpe ach.com/jup1314/2014_02_25rgb0305.jpg // 2014-02-25 19:03 UT http://www.damianpe ach.com/jup1314/2014_02_25rgb0305.jpg
// 2013-05-01 10:29 UT http://astro.christ one.net/jupiter/jupiter2012/jupiter20130501.jpg // 2013-05-01 10:29 UT http://astro.christ one.net/jupiter/jupiter2012/jupiter20130501.jpg
// 2012-10-26 00:12 UT at http://www.lunar-ca ptures.com//jupiter2012_files/121026_JupiterGRS_Tar.jpg // 2012-10-26 00:12 UT at http://www.lunar-ca ptures.com//jupiter2012_files/121026_JupiterGRS_Tar.jpg
// 2011-08-28 02:29 UT at http://www.damianpe ach.com/jup1112/2011_08_28rgb.jpg // 2011-08-28 02:29 UT at http://www.damianpe ach.com/jup1112/2011_08_28rgb.jpg
// stellarium 2h too early: 2010-09-21 23:37 UT http ://www.digitalsky.org.uk/Jupiter/2010-09-21_23-37-30_R-G-B_800.jpg // stellarium 2h too early: 2010-09-21 23:37 UT http ://www.digitalsky.org.uk/Jupiter/2010-09-21_23-37-30_R-G-B_800.jpg
} }
skipping to change at line 751 skipping to change at line 843
p = p->parent; p = p->parent;
} }
} }
} }
// Compute the distance to the given position in heliocentric coordinate (i n AU) // Compute the distance to the given position in heliocentric coordinate (i n AU)
// This is called by SolarSystem::draw() // This is called by SolarSystem::draw()
double Planet::computeDistance(const Vec3d& obsHelioPos) double Planet::computeDistance(const Vec3d& obsHelioPos)
{ {
distance = (obsHelioPos-getHeliocentricEclipticPos()).length(); distance = (obsHelioPos-getHeliocentricEclipticPos()).length();
// improve fps by juggling updates for asteroids. They must be fast // improve fps by juggling updates for asteroids and other minor bod
if close to observer, but can be slow if further away. ies. They must be fast if close to observer, but can be slow if further awa
if (pType == Planet::isAsteroid) y.
deltaJD=distance*StelCore::JD_SECOND; if (pType >= Planet::isAsteroid)
deltaJDE=distance*StelCore::JD_SECOND;
return distance; return distance;
} }
// Get the phase angle (radians) for an observer at pos obsPos in heliocent ric coordinates (dist in AU) // Get the phase angle (radians) for an observer at pos obsPos in heliocent ric coordinates (dist in AU)
double Planet::getPhaseAngle(const Vec3d& obsPos) const double Planet::getPhaseAngle(const Vec3d& obsPos) const
{ {
const double observerRq = obsPos.lengthSquared(); const double observerRq = obsPos.lengthSquared();
const Vec3d& planetHelioPos = getHeliocentricEclipticPos(); const Vec3d& planetHelioPos = getHeliocentricEclipticPos();
const double planetRq = planetHelioPos.lengthSquared(); const double planetRq = planetHelioPos.lengthSquared();
const double observerPlanetRq = (obsPos - planetHelioPos).lengthSqua red(); const double observerPlanetRq = (obsPos - planetHelioPos).lengthSqua red();
skipping to change at line 793 skipping to change at line 885
const double planetRq = planetHelioPos.lengthSquared(); const double planetRq = planetHelioPos.lengthSquared();
const double observerPlanetRq = (obsPos - planetHelioPos).lengthSqua red(); const double observerPlanetRq = (obsPos - planetHelioPos).lengthSqua red();
return std::acos((observerPlanetRq + observerRq - planetRq)/(2.0*sq rt(observerPlanetRq*observerRq))); return std::acos((observerPlanetRq + observerRq - planetRq)/(2.0*sq rt(observerPlanetRq*observerRq)));
} }
// Computation of the visual magnitude (V band) of the planet. // Computation of the visual magnitude (V band) of the planet.
float Planet::getVMagnitude(const StelCore* core) const float Planet::getVMagnitude(const StelCore* core) const
{ {
if (parent == 0) if (parent == 0)
{ {
// sun, compute the apparent magnitude for the absolute mag // Sun, compute the apparent magnitude for the absolute mag
(4.83) and observer's distance (V: 4.83) and observer's distance
// Hint: Absolute Magnitude of the Sun in Several Bands: htt
p://mips.as.arizona.edu/~cnaw/sun.html
const double distParsec = std::sqrt(core->getObserverHelioce ntricEclipticPos().lengthSquared())*AU/PARSEC; const double distParsec = std::sqrt(core->getObserverHelioce ntricEclipticPos().lengthSquared())*AU/PARSEC;
return 4.83 + 5.*(std::log10(distParsec)-1.);
// check how much of it is visible
const SolarSystem* ssm = GETSTELMODULE(SolarSystem);
double shadowFactor = ssm->getEclipseFactor(core);
// See: Hughes, D. W., Brightness during a solar eclipse //
Journal of the British Astronomical Association, vol.110, no.4, p.203-205
// URL: http://adsabs.harvard.edu/abs/2000JBAA..110..203H
if(shadowFactor < 0.000128)
shadowFactor = 0.000128;
return 4.83 + 5.*(std::log10(distParsec)-1.) - 2.5*(std::log
10(shadowFactor));
} }
// Compute the angular phase // Compute the angular phase
const Vec3d& observerHelioPos = core->getObserverHeliocentricEclipti cPos(); const Vec3d& observerHelioPos = core->getObserverHeliocentricEclipti cPos();
const double observerRq = observerHelioPos.lengthSquared(); const double observerRq = observerHelioPos.lengthSquared();
const Vec3d& planetHelioPos = getHeliocentricEclipticPos(); const Vec3d& planetHelioPos = getHeliocentricEclipticPos();
const double planetRq = planetHelioPos.lengthSquared(); const double planetRq = planetHelioPos.lengthSquared();
const double observerPlanetRq = (observerHelioPos - planetHelioPos). lengthSquared(); const double observerPlanetRq = (observerHelioPos - planetHelioPos). lengthSquared();
const double cos_chi = (observerPlanetRq + planetRq - observerRq)/(2 .0*std::sqrt(observerPlanetRq*planetRq)); const double cos_chi = (observerPlanetRq + planetRq - observerRq)/(2 .0*std::sqrt(observerPlanetRq*planetRq));
const double phase = std::acos(cos_chi); const double phase = std::acos(cos_chi);
skipping to change at line 883 skipping to change at line 985
if (englishName=="Venus") if (englishName=="Venus")
return -4.29 + d + 0.09*f1 + 2.39*f1 *f1 - 0.65*f1*f1*f1; return -4.29 + d + 0.09*f1 + 2.39*f1 *f1 - 0.65*f1*f1*f1;
if (englishName=="Mars") if (englishName=="Mars")
return -1.52 + d + 0.016*phaseDeg; return -1.52 + d + 0.016*phaseDeg;
if (englishName=="Jupiter") if (englishName=="Jupiter")
return -9.25 + d + 0.005*phaseDeg; return -9.25 + d + 0.005*phaseDeg;
if (englishName=="Saturn") if (englishName=="Saturn")
{ {
// add rings computation // add rings computation
// implemented from Meeus, Astr.Alg. 1992 // implemented from Meeus, Astr.Alg. 1992
const double jd=core->getJDay(); const double jde=core->getJDE();
const double T=(jd-2451545.0)/36525. const double T=(jde-2451545.0)/36525
0; .0;
const double i=((0.000004*T-0.012998 )*T+28.075216)*M_PI/180.0; const double i=((0.000004*T-0.012998 )*T+28.075216)*M_PI/180.0;
const double Omega=((0.000412*T+1.39 4681)*T+169.508470)*M_PI/180.0; const double Omega=((0.000412*T+1.39 4681)*T+169.508470)*M_PI/180.0;
static SolarSystem *ssystem=GETSTELM ODULE(SolarSystem); static SolarSystem *ssystem=GETSTELM ODULE(SolarSystem);
const Vec3d saturnEarth=getHeliocent ricEclipticPos() - ssystem->getEarth()->getHeliocentricEclipticPos(); const Vec3d saturnEarth=getHeliocent ricEclipticPos() - ssystem->getEarth()->getHeliocentricEclipticPos();
double lambda=atan2(saturnEarth[1], saturnEarth[0]); double lambda=atan2(saturnEarth[1], saturnEarth[0]);
double beta=atan2(saturnEarth[2], st d::sqrt(saturnEarth[0]*saturnEarth[0]+saturnEarth[1]*saturnEarth[1])); double beta=atan2(saturnEarth[2], st d::sqrt(saturnEarth[0]*saturnEarth[0]+saturnEarth[1]*saturnEarth[1]));
const double sinx=sin(i)*cos(beta)*s in(lambda-Omega)-cos(i)*sin(beta); const double sinx=sin(i)*cos(beta)*s in(lambda-Omega)-cos(i)*sin(beta);
double rings = -2.6*sinx + 1.25*sinx *sinx; double rings = -2.6*sinx + 1.25*sinx *sinx;
return -8.88 + d + 0.044*phaseDeg + rings; return -8.88 + d + 0.044*phaseDeg + rings;
} }
skipping to change at line 922 skipping to change at line 1024
if (englishName=="Venus") if (englishName=="Venus")
return -4.0 + d + 0.01322*phaseDeg + 0.0000004247*phaseDeg*phaseDeg*phaseDeg; return -4.0 + d + 0.01322*phaseDeg + 0.0000004247*phaseDeg*phaseDeg*phaseDeg;
if (englishName=="Mars") if (englishName=="Mars")
return -1.3 + d + 0.01486*phaseDeg; return -1.3 + d + 0.01486*phaseDeg;
if (englishName=="Jupiter") if (englishName=="Jupiter")
return -8.93 + d; return -8.93 + d;
if (englishName=="Saturn") if (englishName=="Saturn")
{ {
// add rings computation // add rings computation
// implemented from Meeus, Astr.Alg. 1992 // implemented from Meeus, Astr.Alg. 1992
const double jd=core->getJDay(); const double jde=core->getJDE();
const double T=(jd-2451545.0)/36525. const double T=(jde-2451545.0)/36525
0; .0;
const double i=((0.000004*T-0.012998 )*T+28.075216)*M_PI/180.0; const double i=((0.000004*T-0.012998 )*T+28.075216)*M_PI/180.0;
const double Omega=((0.000412*T+1.39 4681)*T+169.508470)*M_PI/180.0; const double Omega=((0.000412*T+1.39 4681)*T+169.508470)*M_PI/180.0;
SolarSystem *ssystem=GETSTELMODULE(S olarSystem); SolarSystem *ssystem=GETSTELMODULE(S olarSystem);
const Vec3d saturnEarth=getHeliocent ricEclipticPos() - ssystem->getEarth()->getHeliocentricEclipticPos(); const Vec3d saturnEarth=getHeliocent ricEclipticPos() - ssystem->getEarth()->getHeliocentricEclipticPos();
double lambda=atan2(saturnEarth[1], saturnEarth[0]); double lambda=atan2(saturnEarth[1], saturnEarth[0]);
double beta=atan2(saturnEarth[2], st d::sqrt(saturnEarth[0]*saturnEarth[0]+saturnEarth[1]*saturnEarth[1])); double beta=atan2(saturnEarth[2], st d::sqrt(saturnEarth[0]*saturnEarth[0]+saturnEarth[1]*saturnEarth[1]));
const double sinB=sin(i)*cos(beta)*s in(lambda-Omega)-cos(i)*sin(beta); const double sinB=sin(i)*cos(beta)*s in(lambda-Omega)-cos(i)*sin(beta);
double rings = -2.6*fabs(sinB) + 1.2 5*sinB*sinB; // sinx=sinB, saturnicentric latitude of earth. longish, see M eeus. double rings = -2.6*fabs(sinB) + 1.2 5*sinB*sinB; // sinx=sinB, saturnicentric latitude of earth. longish, see M eeus.
return -8.68 + d + 0.044*phaseDeg + rings; return -8.68 + d + 0.044*phaseDeg + rings;
} }
skipping to change at line 960 skipping to change at line 1062
if (englishName=="Venus") if (englishName=="Venus")
return -4.40 + d + 0.0009*phaseDeg + 0.000239*phaseDeg*phaseDeg - 0.00000065*phaseDeg*phaseDeg*phaseDeg; return -4.40 + d + 0.0009*phaseDeg + 0.000239*phaseDeg*phaseDeg - 0.00000065*phaseDeg*phaseDeg*phaseDeg;
if (englishName=="Mars") if (englishName=="Mars")
return -1.52 + d + 0.016*phaseDeg; return -1.52 + d + 0.016*phaseDeg;
if (englishName=="Jupiter") if (englishName=="Jupiter")
return -9.40 + d + 0.005*phaseDeg; return -9.40 + d + 0.005*phaseDeg;
if (englishName=="Saturn") if (englishName=="Saturn")
{ {
// add rings computation // add rings computation
// implemented from Meeus, Astr.Alg. 1992 // implemented from Meeus, Astr.Alg. 1992
const double jd=core->getJDay(); const double jde=core->getJDE();
const double T=(jd-2451545.0)/36525. const double T=(jde-2451545.0)/36525
0; .0;
const double i=((0.000004*T-0.012998 )*T+28.075216)*M_PI/180.0; const double i=((0.000004*T-0.012998 )*T+28.075216)*M_PI/180.0;
const double Omega=((0.000412*T+1.39 4681)*T+169.508470)*M_PI/180.0; const double Omega=((0.000412*T+1.39 4681)*T+169.508470)*M_PI/180.0;
static SolarSystem *ssystem=GETSTELM ODULE(SolarSystem); static SolarSystem *ssystem=GETSTELM ODULE(SolarSystem);
const Vec3d saturnEarth=getHeliocent ricEclipticPos() - ssystem->getEarth()->getHeliocentricEclipticPos(); const Vec3d saturnEarth=getHeliocent ricEclipticPos() - ssystem->getEarth()->getHeliocentricEclipticPos();
double lambda=atan2(saturnEarth[1], saturnEarth[0]); double lambda=atan2(saturnEarth[1], saturnEarth[0]);
double beta=atan2(saturnEarth[2], st d::sqrt(saturnEarth[0]*saturnEarth[0]+saturnEarth[1]*saturnEarth[1])); double beta=atan2(saturnEarth[2], st d::sqrt(saturnEarth[0]*saturnEarth[0]+saturnEarth[1]*saturnEarth[1]));
const double sinB=sin(i)*cos(beta)*s in(lambda-Omega)-cos(i)*sin(beta); const double sinB=sin(i)*cos(beta)*s in(lambda-Omega)-cos(i)*sin(beta);
double rings = -2.6*fabs(sinB) + 1.2 5*sinB*sinB; // sinx=sinB, saturnicentric latitude of earth. longish, see M eeus. double rings = -2.6*fabs(sinB) + 1.2 5*sinB*sinB; // sinx=sinB, saturnicentric latitude of earth. longish, see M eeus.
return -8.88 + d + 0.044*phaseDeg + rings; return -8.88 + d + 0.044*phaseDeg + rings;
} }
skipping to change at line 1013 skipping to change at line 1115
double Planet::getSpheroidAngularSize(const StelCore* core) const double Planet::getSpheroidAngularSize(const StelCore* core) const
{ {
return std::atan2(radius*sphereScale,getJ2000EquatorialPos(core).len gth()) * 180./M_PI; return std::atan2(radius*sphereScale,getJ2000EquatorialPos(core).len gth()) * 180./M_PI;
} }
// Draw the Planet and all the related infos : name, circle etc.. // Draw the Planet and all the related infos : name, circle etc..
void Planet::draw(StelCore* core, float maxMagLabels, const QFont& planetNa meFont) void Planet::draw(StelCore* core, float maxMagLabels, const QFont& planetNa meFont)
{ {
if (hidden) if (hidden)
return; return;
// Exclude drawing if user set a hard limit magnitude.
if (core->getSkyDrawer()->getFlagPlanetMagnitudeLimit() && (getVMagn
itude(core) > core->getSkyDrawer()->getCustomPlanetMagnitudeLimit()))
return;
// Try to improve speed for minor planets: test if visible at all. // Try to improve speed for minor planets: test if visible at all.
// For a full catalog of NEAs (11000 objects), with this and resetti ng deltaJD according to distance, rendering time went 4.5fps->12fps. // For a full catalog of NEAs (11000 objects), with this and resetti ng deltaJD according to distance, rendering time went 4.5fps->12fps.
// TBD: Note that taking away the asteroids at this stage breaks dim -asteroid occultation of stars! // TBD: Note that taking away the asteroids at this stage breaks dim -asteroid occultation of stars!
// Maybe make another configurable flag for those interested? // Maybe make another configurable flag for those interested?
// Problematic: Early-out here of course disables the wanted hint ci rcles for dim asteroids. // Problematic: Early-out here of course disables the wanted hint ci rcles for dim asteroids.
// The line makes hints for asteroids 5 magnitudes below sky limitin g magnitude visible. // The line makes hints for asteroids 5 magnitudes below sky limitin g magnitude visible.
// If asteroid is too faint to be seen, don't bother rendering. (Mas sive speedup if people have hundreds of orbital elements!) // If asteroid is too faint to be seen, don't bother rendering. (Mas sive speedup if people have hundreds of orbital elements!)
if (((getVMagnitude(core)-5.0f) > core->getSkyDrawer()->getLimitMagn itude()) && pType>=Planet::isAsteroid) if (((getVMagnitude(core)-5.0f) > core->getSkyDrawer()->getLimitMagn itude()) && pType>=Planet::isAsteroid)
{ {
return; return;
skipping to change at line 1374 skipping to change at line 1481
delete moonShaderProgram; delete moonShaderProgram;
moonShaderProgram = NULL; moonShaderProgram = NULL;
} }
void Planet::draw3dModel(StelCore* core, StelProjector::ModelViewTranformP transfo, float screenSz, bool drawOnlyRing) void Planet::draw3dModel(StelCore* core, StelProjector::ModelViewTranformP transfo, float screenSz, bool drawOnlyRing)
{ {
// This is the main method drawing a planet 3d model // This is the main method drawing a planet 3d model
// Some work has to be done on this method to make the rendering nic er // Some work has to be done on this method to make the rendering nic er
SolarSystem* ssm = GETSTELMODULE(SolarSystem); SolarSystem* ssm = GETSTELMODULE(SolarSystem);
// Find extinction settings to change colors. The method is rather a
d-hoc.
float extinctedMag=getVMagnitudeWithExtinction(core)-getVMagnitude(c
ore); // this is net value of extinction, in mag.
float magFactorGreen=pow(0.85f, 0.6f*extinctedMag);
float magFactorBlue=pow(0.6f, 0.5f*extinctedMag);
if (screenSz>1.) if (screenSz>1.)
{ {
StelProjector::ModelViewTranformP transfo2 = transfo->clone( ); StelProjector::ModelViewTranformP transfo2 = transfo->clone( );
transfo2->combine(Mat4d::zrotation(M_PI/180*(axisRotation + 90.))); transfo2->combine(Mat4d::zrotation(M_PI/180*(axisRotation + 90.)));
StelPainter* sPainter = new StelPainter(core->getProjection( transfo2)); StelPainter* sPainter = new StelPainter(core->getProjection( transfo2));
if (flagLighting) if (flagLighting)
{ {
// Set the main source of light to be the sun // Set the main source of light to be the sun
Vec3d sunPos(0); Vec3d sunPos(0);
core->getHeliocentricEclipticModelViewTransform()->f orward(sunPos); core->getHeliocentricEclipticModelViewTransform()->f orward(sunPos);
light.position=Vec4f(sunPos[0],sunPos[1],sunPos[2],1 .f); light.position=Vec4f(sunPos[0],sunPos[1],sunPos[2],1 .f);
// Set the light parameters taking sun as the light source // Set the light parameters taking sun as the light source
light.diffuse = Vec4f(1.f,1.f,1.f); light.diffuse = Vec4f(1.f,magFactorGreen*1.f,magFact
light.ambient = Vec4f(0.02f,0.02f,0.02f); orBlue*1.f);
light.ambient = Vec4f(0.02f,magFactorGreen*0.02f,mag
FactorBlue*0.02f);
if (this==ssm->getMoon()) if (this==ssm->getMoon())
{ {
// Special case for the Moon (maybe better u // Special case for the Moon (maybe better u
se 1.5,1.5,1.5,1.0 ?) se 1.5,1.5,1.5,1.0 ?) Was 1.6, but this often is too bright.
light.diffuse = Vec4f(1.6f,1.6f,1.6f,1.f); light.diffuse = Vec4f(1.5f,magFactorGreen*1.
5f,magFactorBlue*1.5f,1.f);
} }
} }
else else
{ {
sPainter->setColor(albedo,albedo,albedo); sPainter->setColor(albedo,magFactorGreen*albedo,magF actorBlue*albedo);
} }
// possibly tint sun's color from extinction. This should de
liberately cause stronger reddening than for the other objects.
if (this==ssm->getSun())
sPainter->setColor(2.f, pow(0.75f, extinctedMag)*2.f
, pow(0.42f, 0.9f*extinctedMag)*2.f);
if (rings) if (rings)
{ {
// The planet has rings, we need to use depth buffer and adjust the clipping planes to avoid // The planet has rings, we need to use depth buffer and adjust the clipping planes to avoid
// reaching the maximum resolution of the depth buff er // reaching the maximum resolution of the depth buff er
const double dist = getEquinoxEquatorialPos(core).le ngth(); const double dist = getEquinoxEquatorialPos(core).le ngth();
double z_near = 0.9*(dist - rings->getSize()); double z_near = 0.9*(dist - rings->getSize());
double z_far = 1.1*(dist + rings->getSize()); double z_far = 1.1*(dist + rings->getSize());
if (z_near < 0.0) z_near = 0.0; if (z_near < 0.0) z_near = 0.0;
double n,f; double n,f;
core->getClippingPlanes(&n,&f); // Save clipping pla nes core->getClippingPlanes(&n,&f); // Save clipping pla nes
skipping to change at line 1436 skipping to change at line 1552
// Draw the halo if it enabled in the ssystem.ini file (+ special ca se for backward compatible for the Sun) // Draw the halo if it enabled in the ssystem.ini file (+ special ca se for backward compatible for the Sun)
if (hasHalo() || this==ssm->getSun()) if (hasHalo() || this==ssm->getSun())
{ {
// Prepare openGL lighting parameters according to luminance // Prepare openGL lighting parameters according to luminance
float surfArcMin2 = getSpheroidAngularSize(core)*60; float surfArcMin2 = getSpheroidAngularSize(core)*60;
surfArcMin2 = surfArcMin2*surfArcMin2*M_PI; // the total ill uminated area in arcmin^2 surfArcMin2 = surfArcMin2*surfArcMin2*M_PI; // the total ill uminated area in arcmin^2
StelPainter sPainter(core->getProjection(StelCore::FrameJ200 0)); StelPainter sPainter(core->getProjection(StelCore::FrameJ200 0));
Vec3d tmp = getJ2000EquatorialPos(core); Vec3d tmp = getJ2000EquatorialPos(core);
core->getSkyDrawer()->postDrawSky3dModel(&sPainter, Vec3f(tm
p[0], tmp[1], tmp[2]), surfArcMin2, getVMagnitudeWithExtinction(core), colo // Find new extincted color for halo. The method is again ra
r); ther ad-hoc, but does not look too bad.
// For the sun, we have again to use the stronger extinction
to avoid color mismatch.
Vec3f haloColorToDraw;
if (this==ssm->getSun())
haloColorToDraw.set(haloColor[0], pow(0.75f, extinct
edMag) * haloColor[1], pow(0.42f, 0.9f*extinctedMag) * haloColor[2]);
else
haloColorToDraw.set(haloColor[0], magFactorGreen * h
aloColor[1], magFactorBlue * haloColor[2]);
core->getSkyDrawer()->postDrawSky3dModel(&sPainter, Vec3f(tm
p[0], tmp[1], tmp[2]), surfArcMin2, getVMagnitudeWithExtinction(core), halo
ColorToDraw);
if ((englishName=="Sun") && (core->getCurrentLocation().plan
etName == "Earth"))
{
float eclipseFactor = ssm->getEclipseFactor(core);
// This alpha ensures 0 for complete sun, 1 for ecli
pse better 1e-10, with a strong increase towards full eclipse. We still nee
d to square it.
float alpha=-0.1f*qMax(-10.0f, (float) std::log10(ec
lipseFactor));
core->getSkyDrawer()->drawSunCorona(&sPainter, Vec3f
(tmp[0], tmp[1], tmp[2]), 512.f/192.f*screenSz, haloColorToDraw, alpha*alph
a);
}
} }
} }
struct Planet3DModel struct Planet3DModel
{ {
QVector<float> vertexArr; QVector<float> vertexArr;
QVector<float> texCoordArr; QVector<float> texCoordArr;
QVector<unsigned short> indiceArr; QVector<unsigned short> indiceArr;
}; };
skipping to change at line 1586 skipping to change at line 1719
QVector<float> projectedVertexArr; QVector<float> projectedVertexArr;
projectedVertexArr.resize(model.vertexArr.size()); projectedVertexArr.resize(model.vertexArr.size());
for (int i=0;i<model.vertexArr.size()/3;++i) for (int i=0;i<model.vertexArr.size()/3;++i)
painter->getProjector()->project(*((Vec3f*)(model.vertexArr. constData()+i*3)), *((Vec3f*)(projectedVertexArr.data()+i*3))); painter->getProjector()->project(*((Vec3f*)(model.vertexArr. constData()+i*3)), *((Vec3f*)(projectedVertexArr.data()+i*3)));
const SolarSystem* ssm = GETSTELMODULE(SolarSystem); const SolarSystem* ssm = GETSTELMODULE(SolarSystem);
if (this==ssm->getSun()) if (this==ssm->getSun())
{ {
texMap->bind(); texMap->bind();
painter->setColor(2, 2, 2); //painter->setColor(2, 2, 0.2); // This is now in draw3dMode l() to apply extinction
painter->setArrays((Vec3f*)projectedVertexArr.constData(), ( Vec2f*)model.texCoordArr.constData()); painter->setArrays((Vec3f*)projectedVertexArr.constData(), ( Vec2f*)model.texCoordArr.constData());
painter->drawFromArray(StelPainter::Triangles, model.indiceA rr.size(), 0, false, model.indiceArr.constData()); painter->drawFromArray(StelPainter::Triangles, model.indiceA rr.size(), 0, false, model.indiceArr.constData());
return; return;
} }
if (planetShaderProgram==NULL) if (planetShaderProgram==NULL)
Planet::initShader(); Planet::initShader();
Q_ASSERT(planetShaderProgram!=NULL); Q_ASSERT(planetShaderProgram!=NULL);
Q_ASSERT(ringPlanetShaderProgram!=NULL); Q_ASSERT(ringPlanetShaderProgram!=NULL);
Q_ASSERT(moonShaderProgram!=NULL); Q_ASSERT(moonShaderProgram!=NULL);
skipping to change at line 1795 skipping to change at line 1928
} }
// draw orbital path of Planet // draw orbital path of Planet
void Planet::drawOrbit(const StelCore* core) void Planet::drawOrbit(const StelCore* core)
{ {
if (!orbitFader.getInterstate()) if (!orbitFader.getInterstate())
return; return;
if (!re.siderealPeriod) if (!re.siderealPeriod)
return; return;
const StelProjectorP prj = core->getProjection(StelCore::FrameHelioc entricEcliptic); const StelProjectorP prj = core->getProjection(StelCore::FrameHelioc entricEclipticJ2000);
StelPainter sPainter(prj); StelPainter sPainter(prj);
// Normal transparency mode // Normal transparency mode
glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
glEnable(GL_BLEND); glEnable(GL_BLEND);
sPainter.setColor(orbitColor[0], orbitColor[1], orbitColor[2], orbit Fader.getInterstate()); sPainter.setColor(orbitColor[0], orbitColor[1], orbitColor[2], orbit Fader.getInterstate());
Vec3d onscreen; Vec3d onscreen;
// special case - use current Planet position as center vertex so th at draws // special case - use current Planet position as center vertex so th at draws
// on it's orbit all the time (since segmented rather than smooth cu rve) // on its orbit all the time (since segmented rather than smooth cur ve)
Vec3d savePos = orbit[ORBIT_SEGMENTS/2]; Vec3d savePos = orbit[ORBIT_SEGMENTS/2];
orbit[ORBIT_SEGMENTS/2]=getHeliocentricEclipticPos(); orbit[ORBIT_SEGMENTS/2]=getHeliocentricEclipticPos();
orbit[ORBIT_SEGMENTS]=orbit[0]; orbit[ORBIT_SEGMENTS]=orbit[0];
int nbIter = closeOrbit ? ORBIT_SEGMENTS : ORBIT_SEGMENTS-1; int nbIter = closeOrbit ? ORBIT_SEGMENTS : ORBIT_SEGMENTS-1;
QVarLengthArray<float, 1024> vertexArray; QVarLengthArray<float, 1024> vertexArray;
sPainter.enableClientStates(true, false, false); sPainter.enableClientStates(true, false, false);
for (int n=0; n<=nbIter; ++n) for (int n=0; n<=nbIter; ++n)
{ {
 End of changes. 78 change blocks. 
111 lines changed or deleted 310 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/