Observability.cpp   Observability.cpp 
skipping to change at line 27 skipping to change at line 27
*/ */
#include <QDebug> #include <QDebug>
#include <QKeyEvent> #include <QKeyEvent>
#include <QMouseEvent> #include <QMouseEvent>
#include <QPixmap> #include <QPixmap>
#include <QSettings> #include <QSettings>
#include <QString> #include <QString>
#include <QTimer> #include <QTimer>
//#include <QtNetwork> // Why do we need a full part of the framwork again?
#include "Observability.hpp" #include "Observability.hpp"
#include "ObservabilityDialog.hpp" #include "ObservabilityDialog.hpp"
#include "Planet.hpp" #include "Planet.hpp"
#include "SolarSystem.hpp" #include "SolarSystem.hpp"
#include "StarMgr.hpp" #include "StarMgr.hpp"
#include "StelActionMgr.hpp" #include "StelActionMgr.hpp"
#include "StelApp.hpp" #include "StelApp.hpp"
#include "StelCore.hpp" #include "StelCore.hpp"
#include "StelFader.hpp" #include "StelFader.hpp"
skipping to change at line 74 skipping to change at line 72
StelPluginInfo info; StelPluginInfo info;
info.id = "Observability"; info.id = "Observability";
info.displayedName = N_("Observability Analysis"); info.displayedName = N_("Observability Analysis");
info.authors = "Ivan Marti-Vidal (Onsala Space Observatory)"; // non -translatable field info.authors = "Ivan Marti-Vidal (Onsala Space Observatory)"; // non -translatable field
info.contact = "i.martividal@gmail.com"; info.contact = "i.martividal@gmail.com";
info.description = N_("Displays an analysis of a selected object's o bservability (rise, set, and transit times) for the current date, as well a s when it is observable through the year. An object is assumed to be observ able if it is above the horizon during a fraction of the night. Also includ ed are the dates of the largest separation from the Sun and acronychal and cosmical rising and setting. (Explanations are provided in the 'About' tab of the plugin's configuration window.)"); info.description = N_("Displays an analysis of a selected object's o bservability (rise, set, and transit times) for the current date, as well a s when it is observable through the year. An object is assumed to be observ able if it is above the horizon during a fraction of the night. Also includ ed are the dates of the largest separation from the Sun and acronychal and cosmical rising and setting. (Explanations are provided in the 'About' tab of the plugin's configuration window.)");
info.version = OBSERVABILITY_PLUGIN_VERSION; info.version = OBSERVABILITY_PLUGIN_VERSION;
return info; return info;
} }
// TODO: Migrate to static const? --BM ==> GZ during JDfix for 0.14: SURE!
// Some useful constants:
const double Observability::Rad2Deg = 180./M_PI; // Convert degrees
into radians
const double Observability::Rad2Hr = 12./M_PI; // Convert hours i
nto radians
const double Observability::UA = AU; // 1.4958e+8;
// Astronomical Unit in Km. ==> HAS BEEN DEFINED IN StelUtils.hpp!
const double Observability::TFrac = 0.9972677595628414; // Convert siderea
l time into Solar time
const double Observability::JDsec = 1./86400.; // A second in day
s. ==> TODO USE StelCore::JD_SECOND instead
const double Observability::halfpi = M_PI * 0.5; // 1.57079632675;
// pi/2
const double Observability::MoonT = 29.530588; // Moon synodic pe
riod (used as first estimate of Full Moon). ==> FIND MORE DEC. PLACES!
const double Observability::RefFullMoon = 2451564.696; // Reference Julian
date of a Full Moon.
const double Observability::MoonPerilune = 0.0024236308; // Smallest Earth-
Moon distance (in AU).
Observability::Observability() Observability::Observability()
: Jan1stJD(0.) : configDialog(new ObservabilityDialog())
, nextFullMoon(0.)
, prevFullMoon(0.)
, GMTShift(0.) , GMTShift(0.)
, Jan1stJD(0.)
, twilightAltRad(0.) , twilightAltRad(0.)
, twilightAltDeg(0.) , twilightAltDeg(0.)
, refractedHorizonAlt(0.)
, horizonAltitude(0.) , horizonAltitude(0.)
, horizonAltDeg(0.) , horizonAltDeg(0.)
, selRA(0.) , selRA(0.)
, selDec(0.) , selDec(0.)
, alti(0.) , alti(0.)
, horizH(0.) , horizH(0.)
, culmAlt(0.) , culmAlt(0.)
, MoonRise(0.) , MoonRise(0.)
, MoonSet(0.) , MoonSet(0.)
, MoonCulm(0.) , MoonCulm(0.)
skipping to change at line 108 skipping to change at line 122
, show_AcroCos(false) , show_AcroCos(false)
, show_Good_Nights(false) , show_Good_Nights(false)
, show_Best_Night(false) , show_Best_Night(false)
, show_Today(false) , show_Today(false)
, show_FullMoon(false) , show_FullMoon(false)
, flagShowReport(false) , flagShowReport(false)
, fontSize(14) , fontSize(14)
, button(NULL) , button(NULL)
{ {
setObjectName("Observability"); setObjectName("Observability");
configDialog = new ObservabilityDialog();
// TODO: Migrate to static const? --BM
// Some useful constants:
Rad2Deg = 180./3.1415927; // Convert degrees into radians
Rad2Hr = 12./3.1415927; // Convert hours into radians
UA = 1.4958e+8; // Astronomical Unit in Km.
TFrac = 0.9972677595628414; // Convert sidereal time into Solar tim
e
JDsec = 1./86400.; // A second in days.
halfpi = 1.57079632675; // pi/2
MoonT = 29.530588; // Moon synodic period (used as first estimate of
Full Moon).
RefFullMoon = 2451564.696; // Reference Julian date of a Full Moon.
MoonPerilune = 0.0024236308; // Smallest Earth-Moon distance (in AU)
.
nextFullMoon = 0.0;
prevFullMoon = 0.0;
refractedHorizonAlt = 0.0;
selName = "";
// Dummy initial values for parameters and data vectors: // Dummy initial values for parameters and data vectors:
mylat = 1000.; mylon = 1000.; mylat = 1000.;
myJD = 0.0; mylon = 1000.;
myJD.first = 0.;
myJD.second = 0.;
curYear = 0; curYear = 0;
isStar = true; isStar = true;
isMoon = false; isMoon = false;
isSun = false; isSun = false;
isScreen = true; isScreen = true;
//Get pointer to the Earth: //Get pointer to the Earth:
PlanetP Earth = GETSTELMODULE(SolarSystem)->getEarth(); PlanetP Earth = GETSTELMODULE(SolarSystem)->getEarth();
myEarth = Earth.data(); myEarth = Earth.data();
// Get pointer to the Moon/Sun: // Get pointer to the Moon/Sun:
PlanetP Moon = GETSTELMODULE(SolarSystem)->getMoon(); PlanetP Moon = GETSTELMODULE(SolarSystem)->getMoon();
myMoon = Moon.data(); myMoon = Moon.data();
// I think this can be done in a more simple way...--BM // I think this can be done in a more simple way...--BM
for (int i=0;i<366;i++) { for (int i=0;i<366;i++) {
sunRA[i] = 0.0; sunDec[i] = 0.0; // sunRA[i] = 0.0; sunDec[i] = 0.0;
objectRA[i] = 0.0; objectDec[i]=0.0; // objectRA[i] = 0.0; objectDec[i]=0.0;
sunSidT[0][i]=0.0; sunSidT[1][i]=0.0; // sunSidT[0][i]=0.0; sunSidT[1][i]=0.0;
objectSidT[0][i]=0.0; objectSidT[1][i]=0.0; // objectSidT[0][i]=0.0; objectSidT[1][i]=0.0;
objectH0[i] = 0.0; // objectH0[i] = 0.0;
yearJD[i] = 0.0; yearJD[i]=QPair<double, double>(0.0, 0.0);
}; };
// GZ Sure:
memset(sunRA, 0, 366*sizeof(double));
memset(sunDec, 0, 366*sizeof(double));
memset(objectRA, 0, 366*sizeof(double));
memset(objectDec, 0, 366*sizeof(double));
memset(sunSidT, 0, 2*366*sizeof(double));
memset(objectSidT, 0, 2*366*sizeof(double));
memset(objectH0, 0, 366*sizeof(double));
} }
Observability::~Observability() Observability::~Observability()
{ {
// Shouldn't this be in the deinit()? --BM // Shouldn't this be in the deinit()? --BM
if (configDialog != NULL) if (configDialog != NULL)
delete configDialog; delete configDialog;
} }
skipping to change at line 281 skipping to change at line 287
// Set the painter: // Set the painter:
StelPainter painter(core->getProjection2d()); StelPainter painter(core->getProjection2d());
painter.setColor(fontColor[0],fontColor[1],fontColor[2],1); painter.setColor(fontColor[0],fontColor[1],fontColor[2],1);
font.setPixelSize(fontSize); font.setPixelSize(fontSize);
painter.setFont(font); painter.setFont(font);
// Get current date, location, and check if there is something selected. // Get current date, location, and check if there is something selected.
double currlat = (core->getCurrentLocation().latitude)/Rad2Deg; double currlat = (core->getCurrentLocation().latitude)/Rad2Deg;
double currlon = (core->getCurrentLocation().longitude)/Rad2Deg; double currlon = (core->getCurrentLocation().longitude)/Rad2Deg;
double currheight = (6371.+(core->getCurrentLocation().altitude)/100 0.)/UA; double currheight = (6371.+(core->getCurrentLocation().altitude)/100 0.)/UA;
double currJD = core->getJDay(); double currJD = core->getJD();
double currJDint; double currJDint;
// GMTShift = StelUtils::getGMTShiftFromQT(currJD)/24.0; // GMTShift = StelUtils::getGMTShiftFromQT(currJD)/24.0;
GMTShift = StelApp::getInstance().getLocaleMgr().getGMTShift(currJD) /24.0; GMTShift = StelApp::getInstance().getLocaleMgr().getGMTShift(currJD) /24.0;
// qDebug() << QString("%1%2 ").arg(GMTShift); // qDebug() << QString("%1%2 ").arg(GMTShift);
double currLocalT = 24.*modf(currJD + GMTShift,&currJDint); double currLocalT = 24.*modf(currJD + GMTShift,&currJDint);
int auxm, auxd, auxy; int auxm, auxd, auxy;
StelUtils::getDateFromJulianDay(currJD, &auxy, &auxm, &auxd); StelUtils::getDateFromJulianDay(currJD, &auxy, &auxm, &auxd);
bool isSource = StelApp::getInstance().getStelObjectMgr().getWasSele cted(); bool isSource = StelApp::getInstance().getStelObjectMgr().getWasSele cted();
bool show_Year = show_Best_Night || show_Good_Nights || show_AcroCos ; bool show_Year = show_Best_Night || show_Good_Nights || show_AcroCos ;
////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////
// NOW WE CHECK THE CHANGED PARAMETERS W.R.T. THE PREVIOUS FRAME: // NOW WE CHECK THE CHANGED PARAMETERS W.R.T. THE PREVIOUS FRAME:
// Update JD. // Update JD.
myJD = currJD; myJD.first = currJD;
myJD.second = currJD + core->computeDeltaT(currJD)/86400.;
// If the year changed, we must recompute the Sun's position for each new d ay: // If the year changed, we must recompute the Sun's position for each new d ay:
if (auxy != curYear) if (auxy != curYear)
{ {
yearChanged = true; yearChanged = true;
curYear = auxy; curYear = auxy;
updateSunData(core); updateSunData(core);
} }
else else
{ {
skipping to change at line 484 skipping to change at line 491
{ {
// Today's ephemeris (rise, set, and transit times) // Today's ephemeris (rise, set, and transit times)
if (!isStar) if (!isStar)
{ {
int type = (isSun) ? 1:0; int type = (isSun) ? 1:0;
type += (isMoon) ? 2:0; type += (isMoon) ? 2:0;
type += (!isSun && !isMoon) ? 3:0; type += (!isSun && !isMoon) ? 3:0;
// Returns false if the calculation fails... // Returns false if the calculation fails...
solvedMoon = calculateSolarSystemEvents(core, type); solvedMoon = calculateSolarSystemEvents(core, type);
currH = qAbs(24.*(MoonCulm-myJD)/TFrac); currH = qAbs(24.*(MoonCulm-myJD.first)/TFrac);
transit = MoonCulm-myJD<0.0; transit = MoonCulm-myJD.first<0.0;
if (solvedMoon) if (solvedMoon)
{ // If failed, Set and Rise will be dummy. { // If failed, Set and Rise will be dummy.
settingTime = qAbs(24.*(MoonSet-myJD)/TFrac) settingTime = qAbs(24.*(MoonSet-myJD.first)/
; TFrac);
risingTime = qAbs(24.*(MoonRise-myJD)/TFrac) risingTime = qAbs(24.*(MoonRise-myJD.first)/
; TFrac);
} }
} }
else if (horizH>0.0) else if (horizH>0.0)
{ // The source is not circumpolar and can be seen from this latitude. { // The source is not circumpolar and can be seen from this latitude.
if ( LocPos[1]>0.0 ) // The source is at the eastern side... if ( LocPos[1]>0.0 ) // The source is at the eastern side...
{ {
if ( currH>horizH ) // ... and below the hor izon. if ( currH>horizH ) // ... and below the hor izon.
{ {
settingTime = 24.-currH-horizH; settingTime = 24.-currH-horizH;
skipping to change at line 926 skipping to change at line 933
}; };
double auxRA = 24.*std::modf(RA/24.,&tempRA); double auxRA = 24.*std::modf(RA/24.,&tempRA);
auxRA += (auxRA<0.0)?24.0:((auxRA>24.0)?-24.0:0.0); auxRA += (auxRA<0.0)?24.0:((auxRA>24.0)?-24.0:0.0);
return auxRA; return auxRA;
} }
//////////////////////////////////// ////////////////////////////////////
QString Observability::formatAsDate(int dayNumber) QString Observability::formatAsDate(int dayNumber)
{ {
int day, month, year; int day, month, year;
StelUtils::getDateFromJulianDay(yearJD[dayNumber], &year, &month, &d ay); StelUtils::getDateFromJulianDay(yearJD[dayNumber].first, &year, &mon th, &day);
QString formatString = (getDateFormat()) ? "%1 %2" : "%2 %1"; QString formatString = (getDateFormat()) ? "%1 %2" : "%2 %1";
QString result = formatString.arg(day).arg(monthNames[month-1]); QString result = formatString.arg(day).arg(monthNames[month-1]);
return result; return result;
} }
/////////////////////////////////////////////// ///////////////////////////////////////////////
// Returns the day and month of year (to put it in format '25 Apr') // Returns the day and month of year (to put it in format '25 Apr')
QString Observability::formatAsDateRange(int startDay, int endDay) QString Observability::formatAsDateRange(int startDay, int endDay)
{ {
int sDay, sMonth, sYear, eDay, eMonth, eYear; int sDay, sMonth, sYear, eDay, eMonth, eYear;
QString range; QString range;
StelUtils::getDateFromJulianDay(yearJD[startDay], &sYear, &sMonth, & StelUtils::getDateFromJulianDay(yearJD[startDay].first, &sYear, &sMo
sDay); nth, &sDay);
StelUtils::getDateFromJulianDay(yearJD[endDay], &eYear, &eMonth, &eD StelUtils::getDateFromJulianDay(yearJD[endDay].first, &eYear, &eMont
ay); h, &eDay);
if (endDay == 0) if (endDay == 0)
{ {
eDay = 31; eDay = 31;
eMonth = 12; eMonth = 12;
} }
if (startDay == 0) if (startDay == 0)
{ {
sDay = 1; sDay = 1;
sMonth = 1; sMonth = 1;
} }
skipping to change at line 996 skipping to change at line 1003
getPlanetCoords(core, myJD, objectRA[0], objectDec[0], true); getPlanetCoords(core, myJD, objectRA[0], objectDec[0], true);
} }
///////////////////////////////////////////////// /////////////////////////////////////////////////
// Computes the Sun's RA and Dec (and the JD) for // Computes the Sun's RA and Dec (and the JD) for
// each day of the current year. // each day of the current year.
void Observability::updateSunData(StelCore* core) void Observability::updateSunData(StelCore* core)
{ {
int day, month, year, sameYear; int day, month, year, sameYear;
// Get current date: // Get current date:
StelUtils::getDateFromJulianDay(myJD,&year,&month,&day); StelUtils::getDateFromJulianDay(myJD.first,&year,&month,&day);
// Get JD for the Jan 1 of current year: // Get JD for the Jan 1 of current year:
StelUtils::getJDFromDate(&Jan1stJD,year,1,1,0,0,0); StelUtils::getJDFromDate(&Jan1stJD,year,1,1,0,0,0);
// Check if we are on a leap year: // Check if we are on a leap year:
StelUtils::getDateFromJulianDay(Jan1stJD+365., &sameYear, &month, &d ay); StelUtils::getDateFromJulianDay(Jan1stJD+365., &sameYear, &month, &d ay);
nDays = (year==sameYear)?366:365; nDays = (year==sameYear)?366:365;
// Compute Earth's position throughout the year: // Compute Earth's position throughout the year:
Vec3d pos, sunPos; Vec3d pos, sunPos;
for (int i=0; i<nDays; i++) for (int i=0; i<nDays; i++)
{ {
yearJD[i] = Jan1stJD + (double)i; yearJD[i].first = Jan1stJD + (double)i;
myEarth->computePosition(yearJD[i]); yearJD[i].second = yearJD[i].first+core->computeDeltaT(yearJ
myEarth->computeTransMatrix(yearJD[i]); D[i].first)/86400.0;
myEarth->computePosition(yearJD[i].second);
myEarth->computeTransMatrix(yearJD[i].first, yearJD[i].secon
d);
pos = myEarth->getHeliocentricEclipticPos(); pos = myEarth->getHeliocentricEclipticPos();
sunPos = core->j2000ToEquinoxEqu((core->matVsop87ToJ2000)*(- pos)); sunPos = core->j2000ToEquinoxEqu((core->matVsop87ToJ2000)*(- pos));
EarthPos[i] = -pos; EarthPos[i] = -pos;
toRADec(sunPos,sunRA[i],sunDec[i]); toRADec(sunPos,sunRA[i],sunDec[i]);
}; };
//Return the Earth to its current time: //Return the Earth to its current time:
myEarth->computePosition(myJD); myEarth->computePosition(myJD.second);
myEarth->computeTransMatrix(myJD); myEarth->computeTransMatrix(myJD.first, myJD.second);
} }
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
//////////////////////////////////////////// ////////////////////////////////////////////
// Computes Sun's Sidereal Times at twilight and culmination: // Computes Sun's Sidereal Times at twilight and culmination:
void Observability::updateSunH() void Observability::updateSunH()
{ {
double tempH, tempH00; double tempH, tempH00;
for (int i=0; i<nDays; i++) for (int i=0; i<nDays; i++)
skipping to change at line 1225 skipping to change at line 1233
// Just return the sign of a double // Just return the sign of a double
double Observability::sign(double d) double Observability::sign(double d)
{ {
return (d<0.0)?-1.0:1.0; return (d<0.0)?-1.0:1.0;
} }
////////////////////////// //////////////////////////
////////////////////////// //////////////////////////
// Get the coordinates of Sun or Moon for a given JD: // Get the coordinates of Sun or Moon for a given JD:
// getBack controls whether Earth and Moon must be returned to their origin al positions after computation. // getBack controls whether Earth and Moon must be returned to their origin al positions after computation.
void Observability::getSunMoonCoords(StelCore *core, double jd, void Observability::getSunMoonCoords(StelCore *core, QPair<double, double>
double &raSun, double &decSun, JD,
double &raMoon, double &decMoon, double &raSun, double &decSun,
double &eclLon, bool getBack) double &raMoon, double &decMoon,
double &eclLon, bool getBack)
//, Vec3d &AltAzVector) //, Vec3d &AltAzVector)
{ {
if (getBack) // Return the Moon and Earth to their current position: if (getBack) // Return the Moon and Earth to their current position:
{ {
myEarth->computePosition(jd); myEarth->computePosition(JD.second);
myEarth->computeTransMatrix(jd); myEarth->computeTransMatrix(JD.first, JD.second);
myMoon->computePosition(jd); myMoon->computePosition(JD.second);
myMoon->computeTransMatrix(jd); myMoon->computeTransMatrix(JD.first, JD.second);
} }
else // Compute coordinates: else // Compute coordinates:
{ {
myEarth->computePosition(jd); myEarth->computePosition(JD.second);
myEarth->computeTransMatrix(jd); myEarth->computeTransMatrix(JD.first, JD.second);
Vec3d earthPos = myEarth->getHeliocentricEclipticPos(); Vec3d earthPos = myEarth->getHeliocentricEclipticPos();
double curSidT; double curSidT;
// Sun coordinates: // Sun coordinates:
Vec3d sunPos = core->j2000ToEquinoxEqu((core->matVsop87ToJ20 00)*(-earthPos)); Vec3d sunPos = core->j2000ToEquinoxEqu((core->matVsop87ToJ20 00)*(-earthPos));
toRADec(sunPos, raSun, decSun); toRADec(sunPos, raSun, decSun);
// Moon coordinates: // Moon coordinates:
curSidT = myEarth->getSiderealTime(jd)/Rad2Deg; curSidT = myEarth->getSiderealTime(JD.first, JD.second)/Rad2 Deg;
RotObserver = (Mat4d::zrotation(curSidT))*ObserverLoc; RotObserver = (Mat4d::zrotation(curSidT))*ObserverLoc;
LocTrans = (core->matVsop87ToJ2000)*(Mat4d::translation(-ear thPos)); LocTrans = (core->matVsop87ToJ2000)*(Mat4d::translation(-ear thPos));
myMoon->computePosition(jd); myMoon->computePosition(JD.second);
myMoon->computeTransMatrix(jd); myMoon->computeTransMatrix(JD.first, JD.second);
Vec3d moonPos = myMoon->getHeliocentricEclipticPos(); Vec3d moonPos = myMoon->getHeliocentricEclipticPos();
sunPos = (core->j2000ToEquinoxEqu(LocTrans*moonPos))-RotObse rver; sunPos = (core->j2000ToEquinoxEqu(LocTrans*moonPos))-RotObse rver;
eclLon = moonPos[0]*earthPos[1] - moonPos[1]*earthPos[0]; eclLon = moonPos[0]*earthPos[1] - moonPos[1]*earthPos[0];
toRADec(sunPos,raMoon,decMoon); toRADec(sunPos,raMoon,decMoon);
}; };
} }
////////////////////////////////////////////// //////////////////////////////////////////////
////////////////////////// //////////////////////////
// Get the Observer-to-Moon distance JD: // Get the Observer-to-Moon distance JD:
// getBack controls whether Earth and Moon must be returned to their origin al positions after computation. // getBack controls whether Earth and Moon must be returned to their origin al positions after computation.
void Observability::getMoonDistance(StelCore *core, double jd, double &dist ance, bool getBack) void Observability::getMoonDistance(StelCore *core, QPair<double, double> J D, double &distance, bool getBack)
{ {
if (getBack) // Return the Moon and Earth to their current position: if (getBack) // Return the Moon and Earth to their current position:
{ {
myEarth->computePosition(jd); myEarth->computePosition(JD.second);
myEarth->computeTransMatrix(jd); myEarth->computeTransMatrix(JD.first, JD.second);
myMoon->computePosition(jd); myMoon->computePosition(JD.second);
myMoon->computeTransMatrix(jd); myMoon->computeTransMatrix(JD.first, JD.second);
} }
else else
{ // Compute coordinates: { // Compute coordinates:
myEarth->computePosition(jd); myEarth->computePosition(JD.second);
myEarth->computeTransMatrix(jd); myEarth->computeTransMatrix(JD.first, JD.second);
Vec3d earthPos = myEarth->getHeliocentricEclipticPos(); Vec3d earthPos = myEarth->getHeliocentricEclipticPos();
// double curSidT; // double curSidT;
// Sun coordinates: // Sun coordinates:
// Pos2 = core->j2000ToEquinoxEqu((core->matVsop87ToJ2000)*(-Po s0)); // Pos2 = core->j2000ToEquinoxEqu((core->matVsop87ToJ2000)*(-Po s0));
// toRADec(Pos2,RASun,DecSun); // toRADec(Pos2,RASun,DecSun);
// Moon coordinates: // Moon coordinates:
// curSidT = myEarth->getSiderealTime(JD)/Rad2Deg; // curSidT = myEarth->getSiderealTime(JD)/Rad2Deg;
// RotObserver = (Mat4d::zrotation(curSidT))*ObserverLoc; // RotObserver = (Mat4d::zrotation(curSidT))*ObserverLoc;
LocTrans = (core->matVsop87ToJ2000)*(Mat4d::translation(-ear thPos)); LocTrans = (core->matVsop87ToJ2000)*(Mat4d::translation(-ear thPos));
myMoon->computePosition(jd); myMoon->computePosition(JD.second);
myMoon->computeTransMatrix(jd); myMoon->computeTransMatrix(JD.first, JD.second);
Pos1 = myMoon->getHeliocentricEclipticPos(); Pos1 = myMoon->getHeliocentricEclipticPos();
Pos2 = (core->j2000ToEquinoxEqu(LocTrans*Pos1)); //-RotObser ver; Pos2 = (core->j2000ToEquinoxEqu(LocTrans*Pos1)); //-RotObser ver;
distance = std::sqrt(Pos2*Pos2); distance = std::sqrt(Pos2*Pos2);
// toRADec(Pos2,RAMoon,DecMoon); // toRADec(Pos2,RAMoon,DecMoon);
}; };
} }
////////////////////////////////////////////// //////////////////////////////////////////////
////////////////////////////////////////////// //////////////////////////////////////////////
// Get the Coords of a planet: // Get the Coords of a planet:
void Observability::getPlanetCoords(StelCore *core, double JD, double &RA, double &Dec, bool getBack) void Observability::getPlanetCoords(StelCore *core, QPair<double, double> J D, double &RA, double &Dec, bool getBack)
{ {
if (getBack) if (getBack)
{ {
// Return the planet to its current time: // Return the planet to its current time:
myPlanet->computePosition(JD); myPlanet->computePosition(JD.second);
myPlanet->computeTransMatrix(JD); myPlanet->computeTransMatrix(JD.first, JD.second);
myEarth->computePosition(JD); myEarth->computePosition(JD.second);
myEarth->computeTransMatrix(JD); myEarth->computeTransMatrix(JD.first, JD.second);
} else } else
{ {
// Compute planet's position: // Compute planet's position:
myPlanet->computePosition(JD); myPlanet->computePosition(JD.second);
myPlanet->computeTransMatrix(JD); myPlanet->computeTransMatrix(JD.first, JD.second);
Pos1 = myPlanet->getHeliocentricEclipticPos(); Pos1 = myPlanet->getHeliocentricEclipticPos();
myEarth->computePosition(JD); myEarth->computePosition(JD.second);
myEarth->computeTransMatrix(JD); myEarth->computeTransMatrix(JD.first, JD.second);
Pos2 = myEarth->getHeliocentricEclipticPos(); Pos2 = myEarth->getHeliocentricEclipticPos();
LocTrans = (core->matVsop87ToJ2000)*(Mat4d::translation(-Pos 2)); LocTrans = (core->matVsop87ToJ2000)*(Mat4d::translation(-Pos 2));
Pos2 = core->j2000ToEquinoxEqu(LocTrans*Pos1); Pos2 = core->j2000ToEquinoxEqu(LocTrans*Pos1);
toRADec(Pos2,RA,Dec); toRADec(Pos2,RA,Dec);
}; };
} }
////////////////////////////////////////////// //////////////////////////////////////////////
////////////////////////////////////////////// //////////////////////////////////////////////
// Solves Moon's, Sun's, or Planet's ephemeris by bissection. // Solves Moon's, Sun's, or Planet's ephemeris by bissection.
bool Observability::calculateSolarSystemEvents(StelCore* core, int bodyType ) bool Observability::calculateSolarSystemEvents(StelCore* core, int bodyType )
{ {
const int NUM_ITER = 100; const int NUM_ITER = 100;
int i; int i;
double hHoriz, ra, dec, raSun, decSun, tempH, tempJd, tempEphH, curS double hHoriz, ra, dec, raSun, decSun, tempH, /* tempJd, */ tempEphH
idT, eclLon; , curSidT, eclLon;
QPair<double, double> tempJd;
//Vec3d Observer; //Vec3d Observer;
hHoriz = calculateHourAngle(mylat, refractedHorizonAlt, selDec); hHoriz = calculateHourAngle(mylat, refractedHorizonAlt, selDec);
bool raises = hHoriz > 0.0; bool raises = hHoriz > 0.0;
// Only recompute ephemeris from second to second (at least) // Only recompute ephemeris from second to second (at least)
// or if the source has changed (i.e., Sun <-> Moon). This saves resources: // or if the source has changed (i.e., Sun <-> Moon). This saves resources:
if (qAbs(myJD-lastJDMoon)>JDsec || lastType!=bodyType || souChanged) if (qAbs(myJD.first-lastJDMoon)>JDsec || lastType!=bodyType || souCh anged)
{ {
// qDebug() << q_("%1 %2 %3 %4").arg(Kind).arg(LastObject) .arg(myJD,0,'f',5).arg(lastJDMoon,0,'f',5); // qDebug() << q_("%1 %2 %3 %4").arg(Kind).arg(LastObject) .arg(myJD,0,'f',5).arg(lastJDMoon,0,'f',5);
lastType = bodyType; lastType = bodyType;
myEarth->computePosition(myJD); myEarth->computePosition(myJD.second);
myEarth->computeTransMatrix(myJD); myEarth->computeTransMatrix(myJD.first, myJD.second);
Vec3d earthPos = myEarth->getHeliocentricEclipticPos(); Vec3d earthPos = myEarth->getHeliocentricEclipticPos();
if (bodyType == 1) // Sun position if (bodyType == 1) // Sun position
{ {
Pos2 = core->j2000ToEquinoxEqu((core->matVsop87ToJ20 00)*(-earthPos)); Pos2 = core->j2000ToEquinoxEqu((core->matVsop87ToJ20 00)*(-earthPos));
} }
else if (bodyType==2) // Moon position else if (bodyType==2) // Moon position
{ {
curSidT = myEarth->getSiderealTime(myJD)/Rad2Deg; curSidT = myEarth->getSiderealTime(myJD.first, myJD. second)/Rad2Deg;
RotObserver = (Mat4d::zrotation(curSidT))*ObserverLo c; RotObserver = (Mat4d::zrotation(curSidT))*ObserverLo c;
LocTrans = (core->matVsop87ToJ2000)*(Mat4d::translat ion(-earthPos)); LocTrans = (core->matVsop87ToJ2000)*(Mat4d::translat ion(-earthPos));
myMoon->computePosition(myJD); myMoon->computePosition(myJD.second);
myMoon->computeTransMatrix(myJD); myMoon->computeTransMatrix(myJD.first, myJD.second);
Pos1 = myMoon->getHeliocentricEclipticPos(); Pos1 = myMoon->getHeliocentricEclipticPos();
Pos2 = (core->j2000ToEquinoxEqu(LocTrans*Pos1))-RotO bserver; Pos2 = (core->j2000ToEquinoxEqu(LocTrans*Pos1))-RotO bserver;
} }
else // Planet position else // Planet position
{ {
myPlanet->computePosition(myJD); myPlanet->computePosition(myJD.second);
myPlanet->computeTransMatrix(myJD); myPlanet->computeTransMatrix(myJD.first, myJD.second
);
Pos1 = myPlanet->getHeliocentricEclipticPos(); Pos1 = myPlanet->getHeliocentricEclipticPos();
LocTrans = (core->matVsop87ToJ2000)*(Mat4d::translat ion(-earthPos)); LocTrans = (core->matVsop87ToJ2000)*(Mat4d::translat ion(-earthPos));
Pos2 = core->j2000ToEquinoxEqu(LocTrans*Pos1); Pos2 = core->j2000ToEquinoxEqu(LocTrans*Pos1);
}; };
toRADec(Pos2,ra,dec); toRADec(Pos2,ra,dec);
Vec3d moonAltAz = core->equinoxEquToAltAz(Pos2, StelCore::Re fractionOff); Vec3d moonAltAz = core->equinoxEquToAltAz(Pos2, StelCore::Re fractionOff);
hasRisen = moonAltAz[2] > refractedHorizonAlt; hasRisen = moonAltAz[2] > refractedHorizonAlt;
// Initial guesses of rise/set/transit times. // Initial guesses of rise/set/transit times.
skipping to change at line 1408 skipping to change at line 1417
if (raises) if (raises)
{ {
if (!hasRisen) if (!hasRisen)
{ {
MoonRise += (MoonRise<0.0)?24.0:0.0; MoonRise += (MoonRise<0.0)?24.0:0.0;
MoonSet -= (MoonSet>0.0)?24.0:0.0; MoonSet -= (MoonSet>0.0)?24.0:0.0;
} }
// Rise time: // Rise time:
tempEphH = MoonRise*TFrac; tempEphH = MoonRise*TFrac;
MoonRise = myJD + (MoonRise/24.); MoonRise = myJD.first + (MoonRise/24.);
for (i=0; i<NUM_ITER; i++) for (i=0; i<NUM_ITER; i++)
{ {
// Get modified coordinates: // Get modified coordinates:
tempJd = MoonRise; tempJd.first = MoonRise;
tempJd.second=tempJd.first+core->computeDelt
aT(tempJd.first)/86400.0;
if (bodyType<3) if (bodyType<3)
{ {
getSunMoonCoords(core, tempJd, getSunMoonCoords(core, tempJd,
raSun, decSun, raSun, decSun,
ra, dec, ra, dec,
eclLon, false); eclLon, false);
} else } else
{ {
getPlanetCoords(core, tempJd, ra, de c, false); getPlanetCoords(core, tempJd, ra, de c, false);
skipping to change at line 1441 skipping to change at line 1451
// H at horizon for mod. coordinates: // H at horizon for mod. coordinates:
hHoriz = calculateHourAngle(mylat,refractedH orizonAlt,dec); hHoriz = calculateHourAngle(mylat,refractedH orizonAlt,dec);
// Compute eph. times for mod. coordinates: // Compute eph. times for mod. coordinates:
tempH = (-hHoriz-Hcurr)*TFrac; tempH = (-hHoriz-Hcurr)*TFrac;
if (hasRisen==false) tempH += (tempH<0.0)?24 .0:0.0; if (hasRisen==false) tempH += (tempH<0.0)?24 .0:0.0;
// Check convergence: // Check convergence:
if (qAbs(tempH-tempEphH)<JDsec) break; if (qAbs(tempH-tempEphH)<JDsec) break;
// Update rise-time estimate: // Update rise-time estimate:
tempEphH = tempH; tempEphH = tempH;
MoonRise = myJD + (tempEphH/24.); MoonRise = myJD.first + (tempEphH/24.);
}; };
// Set time: // Set time:
tempEphH = MoonSet; tempEphH = MoonSet;
MoonSet = myJD + (MoonSet/24.); MoonSet = myJD.first + (MoonSet/24.);
for (i=0; i<NUM_ITER; i++) for (i=0; i<NUM_ITER; i++)
{ {
// Get modified coordinates: // Get modified coordinates:
tempJd = MoonSet; tempJd.first = MoonSet;
tempJd.second=tempJd.first+core->computeDelt
aT(tempJd.first)/86400.0;
if (bodyType < 3) if (bodyType < 3)
getSunMoonCoords(core, tempJd, getSunMoonCoords(core, tempJd,
raSun, decSun, raSun, decSun,
ra, dec, ra, dec,
eclLon, false); eclLon, false);
else else
getPlanetCoords(core, tempJd, ra, de c, false); getPlanetCoords(core, tempJd, ra, de c, false);
if (bodyType==1) {ra = raSun; dec = decSun;} ; if (bodyType==1) {ra = raSun; dec = decSun;} ;
skipping to change at line 1477 skipping to change at line 1488
hHoriz = calculateHourAngle(mylat, refracted HorizonAlt, dec); hHoriz = calculateHourAngle(mylat, refracted HorizonAlt, dec);
// Compute eph. times for mod. coordinates: // Compute eph. times for mod. coordinates:
tempH = (hHoriz-Hcurr)*TFrac; tempH = (hHoriz-Hcurr)*TFrac;
if (!hasRisen) if (!hasRisen)
tempH -= (tempH>0.0)?24.0:0.0; tempH -= (tempH>0.0)?24.0:0.0;
// Check convergence: // Check convergence:
if (qAbs(tempH-tempEphH)<JDsec) if (qAbs(tempH-tempEphH)<JDsec)
break; break;
// Update set-time estimate: // Update set-time estimate:
tempEphH = tempH; tempEphH = tempH;
MoonSet = myJD + (tempEphH/24.); MoonSet = myJD.first + (tempEphH/24.);
}; };
} }
else // Comes from if(raises) else // Comes from if(raises)
{ {
MoonSet = -1.0; MoonSet = -1.0;
MoonRise = -1.0; MoonRise = -1.0;
}; };
// Culmination time: // Culmination time:
tempEphH = MoonCulm; tempEphH = MoonCulm;
MoonCulm = myJD + (MoonCulm/24.); MoonCulm = myJD.first + (MoonCulm/24.);
for (i=0; i<NUM_ITER; i++) for (i=0; i<NUM_ITER; i++)
{ {
// Get modified coordinates: // Get modified coordinates:
tempJd = MoonCulm; tempJd.first = MoonCulm;
tempJd.second=tempJd.first+core->computeDeltaT(tempJ
d.first)/86400.0;
if (bodyType<3) if (bodyType<3)
{ {
getSunMoonCoords(core,tempJd,raSun,decSun,ra ,dec,eclLon,false); getSunMoonCoords(core,tempJd,raSun,decSun,ra ,dec,eclLon,false);
} else } else
{ {
getPlanetCoords(core,tempJd,ra,dec,false); getPlanetCoords(core,tempJd,ra,dec,false);
}; };
if (bodyType==1) {ra = raSun; dec = decSun;}; if (bodyType==1) {ra = raSun; dec = decSun;};
skipping to change at line 1515 skipping to change at line 1527
// Current hour angle at mod. coordinates: // Current hour angle at mod. coordinates:
Hcurr = toUnsignedRA(SidT-ra); Hcurr = toUnsignedRA(SidT-ra);
Hcurr += (LocPos[1]<0.0)?24.0:-24.0; Hcurr += (LocPos[1]<0.0)?24.0:-24.0;
Hcurr -= (Hcurr>12.)?24.0:0.0; Hcurr -= (Hcurr>12.)?24.0:0.0;
// Compute eph. times for mod. coordinates: // Compute eph. times for mod. coordinates:
tempH = -Hcurr*TFrac; tempH = -Hcurr*TFrac;
// Check convergence: // Check convergence:
if (qAbs(tempH-tempEphH)<JDsec) break; if (qAbs(tempH-tempEphH)<JDsec) break;
tempEphH = tempH; tempEphH = tempH;
MoonCulm = myJD + (tempEphH/24.); MoonCulm = myJD.first + (tempEphH/24.);
culmAlt = qAbs(mylat-dec); // 90 - altitude at trans it. culmAlt = qAbs(mylat-dec); // 90 - altitude at trans it.
}; };
// qDebug() << q_("%1").arg(MoonCulm,0,'f',5); // qDebug() << q_("%1").arg(MoonCulm,0,'f',5);
lastJDMoon = myJD; lastJDMoon = myJD.first;
}; // Comes from if (qAbs(myJD-lastJDMoon)>JDsec || LastObject!=Kind ) }; // Comes from if (qAbs(myJD.first-lastJDMoon)>JDsec || LastObject !=Kind)
// Find out the days of Full Moon: // Find out the days of Full Moon:
if (bodyType==2 && show_FullMoon) // || show_SuperMoon)) if (bodyType==2 && show_FullMoon) // || show_SuperMoon))
{ {
// Only estimate date of Full Moon if we have changed Lunar month: // Only estimate date of Full Moon if we have changed Lunar month:
if (myJD > nextFullMoon || myJD < prevFullMoon) if (myJD.first > nextFullMoon || myJD.first < prevFullMoon)
{ {
// Estimate the nearest (in time) Full Moon: // Estimate the nearest (in time) Full Moon:
double nT; double nT;
double dT = std::modf((myJD-RefFullMoon)/MoonT,&nT); double dT = std::modf((myJD.first-RefFullMoon)/MoonT ,&nT);
if (dT>0.5) {nT += 1.0;}; if (dT>0.5) {nT += 1.0;};
if (dT<-0.5) {nT -= 1.0;}; if (dT<-0.5) {nT -= 1.0;};
double TempFullMoon = RefFullMoon + nT*MoonT; double TempFullMoon = RefFullMoon + nT*MoonT;
// Improve the estimate iteratively (Secant method over Lunar-phase vs. time): // Improve the estimate iteratively (Secant method over Lunar-phase vs. time):
dT = 0.1/1440.; // 6 seconds. Our time span for the finite-difference derivative estimate. dT = 0.1/1440.; // 6 seconds. Our time span for the finite-difference derivative estimate.
// double Deriv1, Deriv2; // Variables for temporal use . // double Deriv1, Deriv2; // Variables for temporal use .
double Sec1, Sec2, Temp1, Temp2; // Variables for te QPair<double, double>Sec1, Sec2; // Variables for te
mporal use. mporal use. Contrary to the other cases, we have Sec[12].second=DeltaT[days
], not JDE. This avoids needless re-computing.
double Temp1, Temp2; // Variables for temporal use.
double iniEst1, iniEst2; // JD values that MUST inc lude the solution within them. double iniEst1, iniEst2; // JD values that MUST inc lude the solution within them.
double Phase1; double Phase1;
for (int j=0; j<2; j++) for (int j=0; j<2; j++)
{ // Two steps: one for the previos Full Moon and th e other for the next one. { // Two steps: one for the previos Full Moon and th e other for the next one.
iniEst1 = TempFullMoon - 0.25*MoonT; iniEst1 = TempFullMoon - 0.25*MoonT;
iniEst2 = TempFullMoon + 0.25*MoonT; iniEst2 = TempFullMoon + 0.25*MoonT;
Sec1 = iniEst1; // TempFullMoon - 0.05*MoonT Sec1.first = iniEst1; // TempFullMoon - 0.05
; // Initial estimates of Full-Moon dates *MoonT; // Initial estimates of Full-Moon dates
Sec2 = iniEst2; // TempFullMoon + 0.05*MoonT Sec1.second= core->computeDeltaT(Sec1.first)
; /86400.0; // enough to compute this once.
Sec2.first = iniEst2; // TempFullMoon + 0.05
*MoonT;
Sec2.second= core->computeDeltaT(Sec2.first)
/86400.0; // enough to compute this once.
getSunMoonCoords(core,Sec1,raSun,decSun,ra,d // for the computation calls, we need tempor
ec,eclLon,false); ary QPairs here!
getSunMoonCoords(core,QPair<double, double>(
Sec1.first, Sec1.first+Sec1.second),raSun,decSun,ra,dec,eclLon,false);
Temp1 = eclLon; //Lambda(RA,Dec,RAS,DecS); Temp1 = eclLon; //Lambda(RA,Dec,RAS,DecS);
getSunMoonCoords(core,Sec2,raSun,decSun,ra,d ec,eclLon,false); getSunMoonCoords(core,QPair<double, double>( Sec2.first, Sec2.first+Sec2.second),raSun,decSun,ra,dec,eclLon,false);
Temp2 = eclLon; //Lambda(RA,Dec,RAS,DecS); Temp2 = eclLon; //Lambda(RA,Dec,RAS,DecS);
for (int i=0; i<100; i++) // A limit of 100 iterations. for (int i=0; i<100; i++) // A limit of 100 iterations.
{ {
Phase1 = (Sec2-Sec1)/(Temp1-Temp2)*T Phase1 = (Sec2.first-Sec1.first)/(Te
emp1+Sec1; mp1-Temp2)*Temp1+Sec1.first;
getSunMoonCoords(core,Phase1,raSun,d // The ad-hoc pair needs a DeltaT, u
ecSun,ra,dec,eclLon,false); se the one of Sec1
getSunMoonCoords(core,QPair<double,
double>(Phase1, Phase1+Sec1.second),raSun,decSun,ra,dec,eclLon,false);
if (Temp1*eclLon < 0.0) if (Temp1*eclLon < 0.0)
{ {
Sec2 = Phase1; Sec2.first = Phase1;
Temp2 = eclLon; Temp2 = eclLon;
} else { } else {
Sec1 = Phase1; Sec1.first = Phase1;
Temp1 = eclLon; Temp1 = eclLon;
}; };
// qDebug() << QString("%1 %2 %3 %4 "). arg(Sec1).arg(Sec2).arg(Temp1).arg(Temp2); // qDebug() << QString("%1 %2 %3 %4 "). arg(Sec1).arg(Sec2).arg(Temp1).arg(Temp2);
if (qAbs(Sec2-Sec1) < 10.*dT) // 1 minute accuracy; convergence. if (qAbs(Sec2.first-Sec1.first) < 10 .*dT) // 1 minute accuracy; convergence.
{ {
TempFullMoon = (Sec1+Sec2)/2 .; TempFullMoon = (Sec1.first+S ec2.first)/2.;
// qDebug() << QString("%1%2 ") .arg(TempFullMoon); // qDebug() << QString("%1%2 ") .arg(TempFullMoon);
break; break;
}; };
}; };
if (TempFullMoon > myJD) if (TempFullMoon > myJD.first)
{ {
nextFullMoon = TempFullMoon; nextFullMoon = TempFullMoon;
TempFullMoon -= MoonT; TempFullMoon -= MoonT;
} else } else
{ {
prevFullMoon = TempFullMoon; prevFullMoon = TempFullMoon;
TempFullMoon += MoonT; TempFullMoon += MoonT;
}; };
}; };
 End of changes. 60 change blocks. 
123 lines changed or deleted 158 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/