StelUtils.cpp   StelUtils.cpp 
skipping to change at line 57 skipping to change at line 57
{ {
#ifdef BZR_REVISION #ifdef BZR_REVISION
return QString(PACKAGE_VERSION)+" (BZR r"+BZR_REVISION+")"; return QString(PACKAGE_VERSION)+" (BZR r"+BZR_REVISION+")";
#elif SVN_REVISION #elif SVN_REVISION
return QString(PACKAGE_VERSION)+" (SVN r"+SVN_REVISION+")"; return QString(PACKAGE_VERSION)+" (SVN r"+SVN_REVISION+")";
#else #else
return QString(PACKAGE_VERSION); return QString(PACKAGE_VERSION);
#endif #endif
} }
double hmsToRad(unsigned int h, unsigned int m, double s ) double hmsToRad(const unsigned int h, const unsigned int m, const double s )
{ {
//return (double)M_PI/24.*h*2.+(double)M_PI/12.*m/60.+s*M_PI/43200.; // Wrong formula! --AW //return (double)M_PI/24.*h*2.+(double)M_PI/12.*m/60.+s*M_PI/43200.; // Wrong formula! --AW
return (double)h*M_PI/12.+(double)m*M_PI/10800.+(double)s*M_PI/64800 0.; return (double)h*M_PI/12.+(double)m*M_PI/10800.+(double)s*M_PI/64800 0.;
} }
double dmsToRad(int d, unsigned int m, double s) double dmsToRad(const int d, const unsigned int m, const double s)
{ {
if (d>=0) if (d>=0)
return (double)M_PI/180.*d+(double)M_PI/10800.*m+s*M_PI/6480 00.; return (double)M_PI/180.*d+(double)M_PI/10800.*m+s*M_PI/6480 00.;
return (double)M_PI/180.*d-(double)M_PI/10800.*m-s*M_PI/648000.; return (double)M_PI/180.*d-(double)M_PI/10800.*m-s*M_PI/648000.;
} }
/************************************************************************* /*************************************************************************
Convert an angle in radian to hms Convert an angle in radian to hms
*************************************************************************/ *************************************************************************/
void radToHms(double angle, unsigned int& h, unsigned int& m, double& s) void radToHms(double angle, unsigned int& h, unsigned int& m, double& s)
skipping to change at line 125 skipping to change at line 125
d += 1; d += 1;
else else
d -= 1; d -= 1;
} }
} }
/************************************************************************* /*************************************************************************
Convert an angle in radian to a hms formatted string Convert an angle in radian to a hms formatted string
If the minute and second part are null are too small, don't print them If the minute and second part are null are too small, don't print them
*************************************************************************/ *************************************************************************/
QString radToHmsStrAdapt(double angle) QString radToHmsStrAdapt(const double angle)
{ {
unsigned int h,m; unsigned int h,m;
double s; double s;
QString buf; QString buf;
QTextStream ts(&buf); QTextStream ts(&buf);
StelUtils::radToHms(angle+0.005*M_PI/12/(60*60), h, m, s); StelUtils::radToHms(angle+0.005*M_PI/12/(60*60), h, m, s);
ts << h << 'h'; ts << h << 'h';
if (std::fabs(s*100-(int)s*100)>=1) if (std::fabs(s*100-(int)s*100)>=1)
{ {
ts << m << 'm'; ts << m << 'm';
skipping to change at line 161 skipping to change at line 161
} }
return buf; return buf;
} }
/************************************************************************* /*************************************************************************
Convert an angle in radian to a hms formatted string Convert an angle in radian to a hms formatted string
If decimal is true, output should be like this: " 16h29m55.3s" If decimal is true, output should be like this: " 16h29m55.3s"
If decimal is true, output should be like this: " 16h20m0.4s" If decimal is true, output should be like this: " 16h20m0.4s"
If decimal is false, output should be like this: "0h26m5s" If decimal is false, output should be like this: "0h26m5s"
*************************************************************************/ *************************************************************************/
QString radToHmsStr(double angle, bool decimal) QString radToHmsStr(const double angle, const bool decimal)
{ {
unsigned int h,m; unsigned int h,m;
double s; double s;
StelUtils::radToHms(angle+0.005*M_PI/12/(60*60), h, m, s); StelUtils::radToHms(angle+0.005*M_PI/12/(60*60), h, m, s);
int width, precision; int width, precision;
QString carry; QString carry;
if (decimal) if (decimal)
{ {
width=4; width=4;
precision=1; precision=1;
skipping to change at line 202 skipping to change at line 202
if (h==24 && m==0 && s==0) if (h==24 && m==0 && s==0)
h=0; h=0;
return QString("%1h%2m%3s").arg(h, width).arg(m,2,10,QLatin1Char('0' )).arg(s, 0, 'f', precision); return QString("%1h%2m%3s").arg(h, width).arg(m,2,10,QLatin1Char('0' )).arg(s, 0, 'f', precision);
} }
/************************************************************************* /*************************************************************************
Convert an angle in radian to a dms formatted string Convert an angle in radian to a dms formatted string
If the minute and second part are null are too small, don't print them If the minute and second part are null are too small, don't print them
*************************************************************************/ *************************************************************************/
QString radToDmsStrAdapt(double angle, bool useD) QString radToDmsStrAdapt(const double angle, const bool useD)
{ {
QChar degsign('d'); QChar degsign('d');
if (!useD) if (!useD)
{ {
degsign = 0x00B0; degsign = 0x00B0;
} }
bool sign; bool sign;
unsigned int d,m; unsigned int d,m;
double s; double s;
StelUtils::radToDms(angle+0.005*M_PI/180/(60*60)*(angle<0?-1.:1.), s ign, d, m, s); StelUtils::radToDms(angle+0.005*M_PI/180/(60*60)*(angle<0?-1.:1.), s ign, d, m, s);
skipping to change at line 236 skipping to change at line 236
{ {
os << m << '\''; os << m << '\'';
} }
//qDebug() << "radToDmsStrAdapt(" << angle << ", " << useD << ") = " << str; //qDebug() << "radToDmsStrAdapt(" << angle << ", " << useD << ") = " << str;
return str; return str;
} }
/************************************************************************* /*************************************************************************
Convert an angle in radian to a dms formatted string Convert an angle in radian to a dms formatted string
*************************************************************************/ *************************************************************************/
QString radToDmsStr(double angle, bool decimal, bool useD) QString radToDmsStr(const double angle, const bool decimal, const bool useD )
{ {
QChar degsign('d'); QChar degsign('d');
if (!useD) if (!useD)
{ {
degsign = 0x00B0; degsign = 0x00B0;
} }
bool sign; bool sign;
unsigned int d,m; unsigned int d,m;
double s; double s;
StelUtils::radToDms(angle+0.005*M_PI/180/(60*60)*(angle<0?-1.:1.), s ign, d, m, s); StelUtils::radToDms(angle+0.005*M_PI/180/(60*60)*(angle<0?-1.:1.), s ign, d, m, s);
skipping to change at line 329 skipping to change at line 329
} }
else else
{ {
v[0] = 0.; v[0] = 0.;
v[1] = 0.; v[1] = 0.;
v[2] = 0.; v[2] = 0.;
} }
return v; return v;
} }
void spheToRect(double lng, double lat, Vec3d& v) void spheToRect(const double lng, const double lat, Vec3d& v)
{ {
const double cosLat = cos(lat); const double cosLat = cos(lat);
v.set(cos(lng) * cosLat, sin(lng) * cosLat, sin(lat)); v.set(cos(lng) * cosLat, sin(lng) * cosLat, sin(lat));
} }
void spheToRect(float lng, float lat, Vec3f& v) void spheToRect(const float lng, const float lat, Vec3f& v)
{ {
const double dlng = lng, dlat = lat, cosLat = cos(dlat); const double dlng = lng, dlat = lat, cosLat = cos(dlat);
v.set(cos(dlng) * cosLat, sin(dlng) * cosLat, sin(dlat)); v.set(cos(dlng) * cosLat, sin(dlng) * cosLat, sin(dlat));
} }
void rectToSphe(double *lng, double *lat, const Vec3d& v) void rectToSphe(double *lng, double *lat, const Vec3d& v)
{ {
double r = v.length(); double r = v.length();
*lat = asin(v[2]/r); *lat = asin(v[2]/r);
*lng = atan2(v[1],v[0]); *lng = atan2(v[1],v[0]);
skipping to change at line 372 skipping to change at line 372
void rectToSphe(double *lng, double *lat, const Vec3f& v) void rectToSphe(double *lng, double *lat, const Vec3f& v)
{ {
double r = v.length(); double r = v.length();
*lat = asin(v[2]/r); *lat = asin(v[2]/r);
*lng = atan2(v[1],v[0]); *lng = atan2(v[1],v[0]);
} }
// GZ: some additions. I need those just for quick conversions for text dis play. // GZ: some additions. I need those just for quick conversions for text dis play.
void ctRadec2Ecl(const double raRad, const double decRad, const double eclR ad, double *lambdaRad, double *betaRad) void ctRadec2Ecl(const double raRad, const double decRad, const double eclR ad, double *lambdaRad, double *betaRad)
{ {
*lambdaRad=std::atan2(std::sin(raRad)*std::cos(eclRad)+std::tan(decRad)*s *lambdaRad=std::atan2(std::sin(raRad)*std::cos(eclRad)+std::tan(decR
td::sin(eclRad), std::cos(raRad)); ad)*std::sin(eclRad), std::cos(raRad));
*betaRad=std::asin(std::sin(decRad)*std::cos(eclRad)-std::cos(decRad)*std *betaRad=std::asin(std::sin(decRad)*std::cos(eclRad)-std::cos(decRad
::sin(eclRad)*std::sin(raRad)); )*std::sin(eclRad)*std::sin(raRad));
} }
// GZ: done // GZ: done
double getDecAngle(const QString& str) double getDecAngle(const QString& str)
{ {
QRegExp re1("^\\s*([\\+\\-])?\\s*(\\d+)\\s*([hHDd\xBA])\\s*(\\d+)\\s *['Mm]\\s*(\\d+(\\.\\d+)?)\\s*[\"Ss]\\s*([NSEWnsew])?\\s*$"); // DMS/HMS QRegExp re1("^\\s*([\\+\\-])?\\s*(\\d+)\\s*([hHDd\xBA])\\s*(\\d+)\\s *['Mm]\\s*(\\d+(\\.\\d+)?)\\s*[\"Ss]\\s*([NSEWnsew])?\\s*$"); // DMS/HMS
QRegExp re2("^\\s*([\\+\\-])?\\s*(\\d+(\\.\\d+)?).?([NSEWnsew])?\\s* $"); // Decimal QRegExp re2("^\\s*([\\+\\-])?\\s*(\\d+(\\.\\d+)?).?([NSEWnsew])?\\s* $"); // Decimal
if (re1.exactMatch(str)) if (re1.exactMatch(str))
{ {
skipping to change at line 415 skipping to change at line 415
if (cardinal.toLower() == "s" || cardinal.toLower() == "w" | | neg) if (cardinal.toLower() == "s" || cardinal.toLower() == "w" | | neg)
deg *= -1.; deg *= -1.;
return (deg * 2 * M_PI / 360.); return (deg * 2 * M_PI / 360.);
} }
qDebug() << "getDecAngle failed to parse angle string:" << str; qDebug() << "getDecAngle failed to parse angle string:" << str;
return -0.0; return -0.0;
} }
// Check if a number is a power of 2 // Check if a number is a power of 2
bool isPowerOfTwo(int value) bool isPowerOfTwo(const int value)
{ {
return (value & -value) == value; return (value & -value) == value;
} }
// Return the first power of two greater or equal to the given value // Return the first power of two greater or equal to the given value
int smallestPowerOfTwoGreaterOrEqualTo(int value) int smallestPowerOfTwoGreaterOrEqualTo(const int value)
{ {
#ifndef NDEBUG #ifndef NDEBUG
const int twoTo30 = 1073741824; const int twoTo30 = 1073741824;
Q_ASSERT_X(value <= twoTo30, Q_FUNC_INFO, Q_ASSERT_X(value <= twoTo30, Q_FUNC_INFO,
"Value too large - smallest greater/equal power-of-2 is o ut of range"); "Value too large - smallest greater/equal power-of-2 is o ut of range");
#endif #endif
if(value == 0){return 0;} if(value == 0){return 0;}
int pot=1; int pot=1;
skipping to change at line 443 skipping to change at line 443
return pot; return pot;
} }
QSize smallestPowerOfTwoSizeGreaterOrEqualTo(const QSize base) QSize smallestPowerOfTwoSizeGreaterOrEqualTo(const QSize base)
{ {
return QSize(smallestPowerOfTwoGreaterOrEqualTo(base.width()), return QSize(smallestPowerOfTwoGreaterOrEqualTo(base.width()),
smallestPowerOfTwoGreaterOrEqualTo(base.height())); smallestPowerOfTwoGreaterOrEqualTo(base.height()));
} }
// Return the inverse sinus hyperbolic of z // Return the inverse sinus hyperbolic of z
double asinh(double z) double asinh(const double z)
{ {
return std::log(z+std::sqrt(z*z+1)); return std::log(z+std::sqrt(z*z+1));
} }
/************************************************************************* /*************************************************************************
Convert a QT QDateTime class to julian day Convert a QT QDateTime class to julian day
*************************************************************************/ *************************************************************************/
double qDateTimeToJd(const QDateTime& dateTime) double qDateTimeToJd(const QDateTime& dateTime)
{ {
return (double)(dateTime.date().toJulianDay())+(double)1./(24*60*60* 1000)*QTime().msecsTo(dateTime.time())-0.5; return (double)(dateTime.date().toJulianDay())+(double)1./(24*60*60* 1000)*QTime().msecsTo(dateTime.time())-0.5;
skipping to change at line 465 skipping to change at line 465
QDateTime jdToQDateTime(const double& jd) QDateTime jdToQDateTime(const double& jd)
{ {
int year, month, day; int year, month, day;
getDateFromJulianDay(jd, &year, &month, &day); getDateFromJulianDay(jd, &year, &month, &day);
QDateTime result = QDateTime::fromString(QString("%1.%2.%3").arg(yea r, 4, 10, QLatin1Char('0')).arg(month).arg(day), "yyyy.M.d"); QDateTime result = QDateTime::fromString(QString("%1.%2.%3").arg(yea r, 4, 10, QLatin1Char('0')).arg(month).arg(day), "yyyy.M.d");
result.setTime(jdFractionToQTime(jd)); result.setTime(jdFractionToQTime(jd));
return result; return result;
} }
void getDateFromJulianDay(double jd, int *yy, int *mm, int *dd) void getDateFromJulianDay(const double jd, int *yy, int *mm, int *dd)
{ {
/* /*
* This algorithm is taken from * This algorithm is taken from
* "Numerical Recipes in c, 2nd Ed." (1992), pp. 14-15 * "Numerical Recipes in c, 2nd Ed." (1992), pp. 14-15
* and converted to integer math. * and converted to integer math.
* The electronic version of the book is freely available * The electronic version of the book is freely available
* at http://www.nr.com/ , under "Obsolete Versions - Older * at http://www.nr.com/ , under "Obsolete Versions - Older
* book and code versions. * book and code versions.
*/ */
skipping to change at line 528 skipping to change at line 528
if (*mm > 2) if (*mm > 2)
{ {
--(*yy); --(*yy);
} }
if (julian < 0) if (julian < 0)
{ {
*yy -= 100 * (1 - julian / 36525); *yy -= 100 * (1 - julian / 36525);
} }
} }
void getTimeFromJulianDay(double julianDay, int *hour, int *minute, int *se cond) void getTimeFromJulianDay(const double julianDay, int *hour, int *minute, i nt *second)
{ {
double frac = julianDay - (floor(julianDay)); double frac = julianDay - (floor(julianDay));
int s = (int)floor((frac * 24.0 * 60.0 * 60.0) + 0.0001); // add co nstant to fix floating-point truncation error int s = (int)floor((frac * 24.0 * 60.0 * 60.0) + 0.0001); // add co nstant to fix floating-point truncation error
*hour = ((s / (60 * 60))+12)%24; *hour = ((s / (60 * 60))+12)%24;
*minute = (s/(60))%60; *minute = (s/(60))%60;
*second = s % 60; *second = s % 60;
} }
QString julianDayToISO8601String(double jd) QString julianDayToISO8601String(const double jd)
{ {
int year, month, day, hour, minute, second; int year, month, day, hour, minute, second;
getDateFromJulianDay(jd, &year, &month, &day); getDateFromJulianDay(jd, &year, &month, &day);
getTimeFromJulianDay(jd, &hour, &minute, &second); getTimeFromJulianDay(jd, &hour, &minute, &second);
QString res = QString("%1-%2-%3T%4:%5:%6") QString res = QString("%1-%2-%3T%4:%5:%6")
.arg((year >= 0 ? year : -1* year),4,10,QLa tin1Char('0')) .arg((year >= 0 ? year : -1* year),4,10,QLa tin1Char('0'))
.arg(month,2,10,QLatin1Char('0')) .arg(month,2,10,QLatin1Char('0'))
.arg(day,2,10,QLatin1Char('0')) .arg(day,2,10,QLatin1Char('0'))
.arg(hour,2,10,QLatin1Char('0')) .arg(hour,2,10,QLatin1Char('0'))
.arg(minute,2,10,QLatin1Char('0')) .arg(minute,2,10,QLatin1Char('0'))
.arg(second,2,10,QLatin1Char('0')); .arg(second,2,10,QLatin1Char('0'));
if (year < 0) if (year < 0)
{ {
res.prepend("-"); res.prepend("-");
} }
return res; return res;
} }
// Format the date per the fmt. // Format the date per the fmt.
QString localeDateString(int year, int month, int day, int dayOfWeek, QStri ng fmt) QString localeDateString(const int year, const int month, const int day, co nst int dayOfWeek, const QString fmt)
{ {
/* we have to handle the year zero, and the years before qdatetime c an represent. */ /* we have to handle the year zero, and the years before qdatetime c an represent. */
const QLatin1Char quote('\''); const QLatin1Char quote('\'');
QString out; QString out;
int quotestartedat = -1; int quotestartedat = -1;
for (int i = 0; i < (int)fmt.length(); i++) for (int i = 0; i < (int)fmt.length(); i++)
{ {
if (fmt.at(i) == quote) if (fmt.at(i) == quote)
{ {
skipping to change at line 675 skipping to change at line 675
out += fmt.at(i); out += fmt.at(i);
} }
} }
return out; return out;
} }
//! try to get a reasonable locale date string from the system, trying to w ork around //! try to get a reasonable locale date string from the system, trying to w ork around
//! limitations of qdatetime for large dates in the past. see QDateTime::t oString(). //! limitations of qdatetime for large dates in the past. see QDateTime::t oString().
QString localeDateString(int year, int month, int day, int dayOfWeek) QString localeDateString(const int year, const int month, const int day, co nst int dayOfWeek)
{ {
// try the QDateTime first // try the QDateTime first
QDate test(year, month, day); QDate test(year, month, day);
// try to avoid QDate's non-astronomical time here, don't do BCE or year 0. // try to avoid QDate's non-astronomical time here, don't do BCE or year 0.
if (year > 0 && test.isValid() && !test.toString(Qt::LocaleDate).isE mpty()) if (year > 0 && test.isValid() && !test.toString(Qt::LocaleDate).isE mpty())
{ {
return test.toString(Qt::LocaleDate); return test.toString(Qt::LocaleDate);
} }
skipping to change at line 705 skipping to change at line 705
double getJDFromSystem() double getJDFromSystem()
{ {
return qDateTimeToJd(QDateTime::currentDateTime().toUTC()); return qDateTimeToJd(QDateTime::currentDateTime().toUTC());
} }
double qTimeToJDFraction(const QTime& time) double qTimeToJDFraction(const QTime& time)
{ {
return (double)1./(24*60*60*1000)*QTime().msecsTo(time)-0.5; return (double)1./(24*60*60*1000)*QTime().msecsTo(time)-0.5;
} }
QTime jdFractionToQTime(double jd) QTime jdFractionToQTime(const double jd)
{ {
double decHours = std::fmod(jd+0.5, 1.0); double decHours = std::fmod(jd+0.5, 1.0);
int hours = (int)(decHours/0.041666666666666666666); int hours = (int)(decHours/0.041666666666666666666);
int mins = (int)((decHours-(hours*0.041666666666666666666))/0.000694 44444444444444444); int mins = (int)((decHours-(hours*0.041666666666666666666))/0.000694 44444444444444444);
return QTime::fromString(QString("%1.%2").arg(hours).arg(mins), "h.m "); return QTime::fromString(QString("%1.%2").arg(hours).arg(mins), "h.m ");
} }
// Use Qt's own sense of time and offset instead of platform specific code. // Use Qt's own sense of time and offset instead of platform specific code.
float getGMTShiftFromQT(double JD) float getGMTShiftFromQT(const double JD)
{ {
int year, month, day, hour, minute, second; int year, month, day, hour, minute, second;
getDateFromJulianDay(JD, &year, &month, &day); getDateFromJulianDay(JD, &year, &month, &day);
getTimeFromJulianDay(JD, &hour, &minute, &second); getTimeFromJulianDay(JD, &hour, &minute, &second);
// as analogous to second statement in getJDFromDate, nkerr // as analogous to second statement in getJDFromDate, nkerr
if ( year <= 0 ) if ( year <= 0 )
{ {
year = year - 1; year = year - 1;
} }
//getTime/DateFromJulianDay returns UTC time, not local time //getTime/DateFromJulianDay returns UTC time, not local time
skipping to change at line 743 skipping to change at line 743
//Both timezones should be interpreted as UTC because secsTo() conve rts both //Both timezones should be interpreted as UTC because secsTo() conve rts both
//times to UTC if their zones have different daylight saving time ru les. //times to UTC if their zones have different daylight saving time ru les.
local.setTimeSpec(Qt::UTC); local.setTimeSpec(Qt::UTC);
int shiftInSeconds = universal.secsTo(local); int shiftInSeconds = universal.secsTo(local);
float shiftInHours = shiftInSeconds / 3600.0f; float shiftInHours = shiftInSeconds / 3600.0f;
return shiftInHours; return shiftInHours;
} }
// UTC ! // UTC !
bool getJDFromDate(double* newjd, int y, int m, int d, int h, int min, int s) bool getJDFromDate(double* newjd, const int y, const int m, const int d, co nst int h, const int min, const int s)
{ {
static const long IGREG2 = 15+31L*(10+12L*1582); static const long IGREG2 = 15+31L*(10+12L*1582);
double deltaTime = (h / 24.0) + (min / (24.0*60.0)) + (s / (24.0 * 6 0.0 * 60.0)) - 0.5; double deltaTime = (h / 24.0) + (min / (24.0*60.0)) + (s / (24.0 * 6 0.0 * 60.0)) - 0.5;
QDate test((y <= 0 ? y-1 : y), m, d); QDate test((y <= 0 ? y-1 : y), m, d);
// if QDate will oblige, do so. // if QDate will oblige, do so.
if ( test.isValid() ) if ( test.isValid() )
{ {
double qdjd = (double)test.toJulianDay(); double qdjd = (double)test.toJulianDay();
qdjd += deltaTime; qdjd += deltaTime;
*newjd = qdjd; *newjd = qdjd;
skipping to change at line 806 skipping to change at line 806
ljul += 2 - lcc + lee; ljul += 2 - lcc + lee;
} }
double jd = (double)ljul; double jd = (double)ljul;
jd += deltaTime; jd += deltaTime;
*newjd = jd; *newjd = jd;
return true; return true;
} }
return false; return false;
} }
double getJDFromDate_alg2(int y, int m, int d, int h, int min, int s) double getJDFromDate_alg2(const int y, const int m, const int d, const int h, const int min, const int s)
{ {
double extra = (100.0* y) + m - 190002.5; double extra = (100.0* y) + m - 190002.5;
double rjd = 367.0 * y; double rjd = 367.0 * y;
rjd -= floor(7.0*(y+floor((m+9.0)/12.0))/4.0); rjd -= floor(7.0*(y+floor((m+9.0)/12.0))/4.0);
rjd += floor(275.0*m/9.0) ; rjd += floor(275.0*m/9.0) ;
rjd += d; rjd += d;
rjd += (h + (min + s/60.0)/60.)/24.0; rjd += (h + (min + s/60.0)/60.)/24.0;
rjd += 1721013.5; rjd += 1721013.5;
rjd -= 0.5*extra/std::fabs(extra); rjd -= 0.5*extra/std::fabs(extra);
rjd += 0.5; rjd += 0.5;
return rjd; return rjd;
} }
int numberOfDaysInMonthInYear(int month, int year) int numberOfDaysInMonthInYear(const int month, const int year)
{ {
switch(month) switch(month)
{ {
case 1: case 1:
case 3: case 3:
case 5: case 5:
case 7: case 7:
case 8: case 8:
case 10: case 10:
case 12: case 12:
skipping to change at line 1050 skipping to change at line 1050
error = error || !ok; error = error || !ok;
*s = finalRe.capturedTexts().at(6).toFloat(&ok); *s = finalRe.capturedTexts().at(6).toFloat(&ok);
error = error || !ok; error = error || !ok;
if (!error) if (!error)
return true; return true;
} }
return false; return false;
} }
// Calculate and getting sidereal period in days from semi-major axis // Calculate and getting sidereal period in days from semi-major axis
double calculateSiderealPeriod(double SemiMajorAxis) double calculateSiderealPeriod(const double SemiMajorAxis)
{ {
// Calculate semi-major axis in meters // Calculate semi-major axis in meters
double a = AU*1000*SemiMajorAxis; double a = AU*1000*SemiMajorAxis;
// Calculate orbital period in seconds // Calculate orbital period in seconds
// Here 1.32712440018e20 is heliocentric gravitational constant // Here 1.32712440018e20 is heliocentric gravitational constant
double period = 2*M_PI*std::sqrt(a*a*a/1.32712440018e20); double period = 2*M_PI*std::sqrt(a*a*a/1.32712440018e20);
return period/86400; // return period in days return period/86400; // return period in days
} }
// Calculate duration of mean solar day QString hoursToHmsStr(const double hours)
double calculateSolarDay(double siderealPeriod, double siderealDay, bool fo
rwardDirection)
{
double coeff;
if (forwardDirection)
coeff = (siderealPeriod + 1)/siderealPeriod;
else
coeff = -1 * (siderealPeriod - 1)/siderealPeriod;
return siderealDay/coeff;
}
QString hoursToHmsStr(double hours)
{ {
int h = (int)hours; int h = (int)hours;
int m = (int)((hours-h)*60); int m = (int)((std::abs(hours)-std::abs(h))*60);
float s = (((hours-h)*60)-m)*60; float s = (((std::abs(hours)-std::abs(h))*60)-m)*60;
return QString("%1h%2m%3s").arg(h).arg(m).arg(QString::number(s, 'f' , 1)); return QString("%1h%2m%3s").arg(h).arg(m).arg(QString::number(s, 'f' , 1));
} }
#if defined(Q_OS_MAC) || defined(Q_OS_LINUX) || defined(Q_OS_FREEBSD) #if defined(Q_OS_MAC) || defined(Q_OS_LINUX) || defined(Q_OS_FREEBSD)
#include <sys/time.h> #include <sys/time.h>
#else #else
#include <time.h> #include <time.h>
#endif #endif
skipping to change at line 1116 skipping to change at line 1104
} }
//! Time when the program execution started. //! Time when the program execution started.
long double startTime = getTime(); long double startTime = getTime();
long double secondsSinceStart() long double secondsSinceStart()
{ {
return getTime() - startTime; return getTime() - startTime;
} }
double decYear2DeltaT(double y) /* /////////////////// DELTA T VARIANTS
// For the standard epochs for many formulae, we use
// J2000.0=2000-jan-1.5=2451545.0,
// 1900.0=1900-jan-0.5=2415020.0
// 1820.0=1820-jan-0.5=2385800.0
// 1810.0=1810-jan-0.5=2382148.0
// 1800.0=1800-jan-0.5=2378496.0
// 1735.0=1735-jan-0.5=2354755.0
// 1625.0=1625-jan-0.5=2314579.0
*/
double decYear2DeltaT(const double y)
{ {
// Note: the method here is adapted from Adapted from // Note: the method here is adapted from
// "Five Millennium Canon of Solar Eclipses" [Espenak and Meeus] // "Five Millennium Canon of Solar Eclipses" [Espenak and Meeus, 200
6]
// A summary is described here: // A summary is described here:
// http://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html // http://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html
// GZ: I replaced the std::pow() calls by Horner's scheme with rever
sed factors, it's more accurate and efficient.
// Old code left for readability, but can also be deleted.
// set the default value for Delta T // set the default value for Delta T
double u = (y-1820)/100.; double u = (y-1820)/100.;
double r = (-20 + 32 * std::pow(u,2.0)); double r = (-20 + 32 * u * u);
if (y < -500) if (y < -500)
{ {
u = (y-1820)/100.; // values are equal to defaults!
r = (-20 + 32 * std::pow(u,2.0)); // u = (y-1820)/100.;
// r = (-20 + 32 * u * u);
} }
else if (y < 500) else if (y < 500)
{ {
u = y/100; u = y/100;
r = (10583.6 - 1014.41 * u + 33.78311 * std::pow(u,2) - 5.95 //r = (10583.6 - 1014.41 * u + 33.78311 * std::pow(u,2) - 5.
2053 * std::pow(u,3) 952053 * std::pow(u,3)
- 0.1798452 * std::pow(u,4) + 0.022174192 * std::pow( // - 0.1798452 * std::pow(u,4) + 0.022174192 * std::po
u,5) + 0.0090316521 * std::pow(u,6)); w(u,5) + 0.0090316521 * std::pow(u,6));
r = (((((0.0090316521*u + 0.022174192)*u - 0.1798452)*u - 5.
952053)*u + 33.78311)*u -1014.41)*u +10583.6;
} }
else if (y < 1600) else if (y < 1600)
{ {
u = (y-1000)/100; u = (y-1000)/100;
r = (1574.2 - 556.01 * u + 71.23472 * std::pow(u,2) + 0.3197 //r = (1574.2 - 556.01 * u + 71.23472 * std::pow(u,2) + 0.31
81 * std::pow(u,3) 9781 * std::pow(u,3)
- 0.8503463 * std::pow(u,4) - 0.005050998 * std::pow( // - 0.8503463 * std::pow(u,4) - 0.005050998 * std::po
u,5) + 0.0083572073 * std::pow(u,6)); w(u,5) + 0.0083572073 * std::pow(u,6));
r = (((((0.0083572073*u - 0.005050998)*u -0.8503463)*u +0.31
9781)*u + 71.23472)*u -556.01)*u + 1574.2;
} }
else if (y < 1700) else if (y < 1700)
{ {
double t = y - 1600; double t = y - 1600;
r = (120 - 0.9808 * t - 0.01532 * std::pow(t,2) + std::pow(t //r = (120 - 0.9808 * t - 0.01532 * std::pow(t,2) + std::pow
,3) / 7129); (t,3) / 7129);
r = ((t/7129.0 - 0.01532)*t - 0.9808)*t +120.0;
} }
else if (y < 1800) else if (y < 1800)
{ {
double t = y - 1700; double t = y - 1700;
r = (8.83 + 0.1603 * t - 0.0059285 * std::pow(t,2) + 0.00013 //r = (8.83 + 0.1603 * t - 0.0059285 * std::pow(t,2) + 0.000
336 * std::pow(t,3) - std::pow(t,4) / 1174000); 13336 * std::pow(t,3) - std::pow(t,4) / 1174000);
r = (((-t/1174000.0 + 0.00013336)*t - 0.0059285)*t + 0.1603)
*t +8.83;
} }
else if (y < 1860) else if (y < 1860)
{ {
double t = y - 1800; double t = y - 1800;
r = (13.72 - 0.332447 * t + 0.0068612 * std::pow(t,2) + 0.00 //r = (13.72 - 0.332447 * t + 0.0068612 * std::pow(t,2) + 0.
41116 * std::pow(t,3) - 0.00037436 * std::pow(t,4) 0041116 * std::pow(t,3) - 0.00037436 * std::pow(t,4)
+ 0.0000121272 * std::pow(t,5) - 0.0000001699 * std:: // + 0.0000121272 * std::pow(t,5) - 0.0000001699 * std
pow(t,6) + 0.000000000875 * std::pow(t,7)); ::pow(t,6) + 0.000000000875 * std::pow(t,7));
r = ((((((.000000000875*t -.0000001699)*t + 0.0000121272)*t
- 0.00037436)*t + 0.0041116)*t + 0.0068612)*t - 0.332447)*t +13.72;
} }
else if (y < 1900) else if (y < 1900)
{ {
double t = y - 1860; double t = y - 1860;
r = (7.62 + 0.5737 * t - 0.251754 * std::pow(t,2) + 0.016806 //r = (7.62 + 0.5737 * t - 0.251754 * std::pow(t,2) + 0.0168
68 * std::pow(t,3) 0668 * std::pow(t,3)
-0.0004473624 * std::pow(t,4) + std::pow(t,5) / 2331 // -0.0004473624 * std::pow(t,4) + std::pow(t,5) / 2331
74); 74);
r = ((((t/233174.0 -0.0004473624)*t + 0.01680668)*t - 0.2517
54)*t + 0.5737)*t + 7.62;
} }
else if (y < 1920) else if (y < 1920)
{ {
double t = y - 1900; double t = y - 1900;
r = (-2.79 + 1.494119 * t - 0.0598939 * std::pow(t,2) + 0.00 //r = (-2.79 + 1.494119 * t - 0.0598939 * std::pow(t,2) + 0.
61966 * std::pow(t,3) - 0.000197 * std::pow(t,4)); 0061966 * std::pow(t,3) - 0.000197 * std::pow(t,4));
r = (((-0.000197*t + 0.0061966)*t - 0.0598939)*t + 1.494119)
*t -2.79;
} }
else if (y < 1941) else if (y < 1941)
{ {
double t = y - 1920; double t = y - 1920;
r = (21.20 + 0.84493*t - 0.076100 * std::pow(t,2) + 0.002093 //r = (21.20 + 0.84493*t - 0.076100 * std::pow(t,2) + 0.0020
6 * std::pow(t,3)); 936 * std::pow(t,3));
r = ((0.0020936*t - 0.076100)*t+ 0.84493)*t +21.20;
} }
else if (y < 1961) else if (y < 1961)
{ {
double t = y - 1950; double t = y - 1950;
r = (29.07 + 0.407*t - std::pow(t,2)/233 + std::pow(t,3) / 2 //r = (29.07 + 0.407*t - std::pow(t,2)/233 + std::pow(t,3) /
547); 2547);
r = ((t/2547.0 -1.0/233.0)*t + 0.407)*t +29.07;
} }
else if (y < 1986) else if (y < 1986)
{ {
double t = y - 1975; double t = y - 1975;
r = (45.45 + 1.067*t - std::pow(t,2)/260 - std::pow(t,3) / 7 //r = (45.45 + 1.067*t - std::pow(t,2)/260 - std::pow(t,3) /
18); 718);
r = ((-t/718.0 +1/260.0)*t + 1.067)*t +45.45;
} }
else if (y < 2005) else if (y < 2005)
{ {
double t = y - 2000; double t = y - 2000;
r = (63.86 + 0.3345 * t - 0.060374 * std::pow(t,2) + 0.00172 //r = (63.86 + 0.3345 * t - 0.060374 * std::pow(t,2) + 0.001
75 * std::pow(t,3) + 0.000651814 * std::pow(t,4) 7275 * std::pow(t,3) + 0.000651814 * std::pow(t,4) + 0.00002373599 * std::p
+ 0.00002373599 * std::pow(t,5)); ow(t,5));
r = ((((0.00002373599*t + 0.000651814)*t + 0.0017275)*t - 0.
060374)*t + 0.3345)*t +63.86;
} }
else if (y < 2050) else if (y < 2050)
{ {
double t = y - 2000; double t = y - 2000;
r = (62.92 + 0.32217 * t + 0.005589 * std::pow(t,2)); //r = (62.92 + 0.32217 * t + 0.005589 * std::pow(t,2));
r = (0.005589*t +0.32217)*t + 62.92;
} }
else if (y < 2150) else if (y < 2150)
{ {
r = (-20 + 32 * std::pow((y-1820)/100,2) - 0.5628 * (2150 - //r = (-20 + 32 * std::pow((y-1820)/100,2) - 0.5628 * (2150
y)); - y));
// r has been precomputed before, just add the term patching
the discontinuity
r -= 0.5628*(2150.0-y);
} }
return r; return r;
} }
double getDeltaT(double jDay) // Implementation of algorithm by Espenak & Meeus (2006) for DeltaT computa
tion
double getDeltaTByEspenakMeeus(const double jDay)
{ {
int year, month, day; int year, month, day;
double moon = 0.;
getDateFromJulianDay(jDay, &year, &month, &day); getDateFromJulianDay(jDay, &year, &month, &day);
if (year<1955 or year>2005)
moon = getMoonSecularAcceleration(jDay);
// approximate "decimal year" = year + (month - 0.5)/12 // approximate "decimal year" = year + (month - 0.5)/12
return decYear2DeltaT(year + (month - 0.5)/12)+moon; //return decYear2DeltaT(year + (month - 0.5)/12);
return decYear2DeltaT(year+((month-1)*30.5+day/31*30.5)/366);
}
// Implementation of algorithm by Schoch (1931) for DeltaT computation
double getDeltaTBySchoch(const double jDay)
{
double u=(jDay-2378496.0)/36525.0; // (1800-jan-0.5)
return -36.28 + 36.28*std::pow(u,2);
}
// Implementation of algorithm by Clemence (1948) for DeltaT computation
double getDeltaTByClemence(const double jDay)
{
double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5)
return +8.72 + 26.75*u + 11.22*std::pow(u,2);
}
// Implementation of algorithm by IAU (1952) for DeltaT computation
double getDeltaTByIAU(const double jDay)
{
double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5)
// TODO: Calculate Moon's longitude fluctuation
return (29.949*u +72.3165)*u +24.349 /* + 1.821*b*/ ;
}
// Implementation of algorithm by Astronomical Ephemeris (1960) for DeltaT
computation, also used by Mucke&Meeus, Canon of Solar Eclipses, Vienna 1983
double getDeltaTByAstronomicalEphemeris(const double jDay)
{
double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5)
// TODO: Calculate Moon's longitude fluctuation
// Note: also Mucke&Meeus 1983 ignore b
return (29.950*u +72.318)*u +24.349 /* + 1.82144*b */ ;
}
// Implementation of algorithm by Tuckerman (1962, 1964) & Goldstine (1973)
for DeltaT computation
double getDeltaTByTuckermanGoldstine(const double jDay)
{
double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5)
return (36.79*u +35.06)*u + 4.87;
}
// Implementation of algorithm by Muller & Stephenson (1975) for DeltaT com
putation
double getDeltaTByMullerStephenson(const double jDay)
{
double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5)
return (45.78*u +120.38)*u +66.0;
}
// Implementation of algorithm by Stephenson (1978) for DeltaT computation
double getDeltaTByStephenson1978(const double jDay)
{
double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5)
return (38.30*u +114.0)*u +20.0;
}
// Implementation of algorithm by Stephenson (1997) for DeltaT computation
double getDeltaTByStephenson1997(const double jDay)
{
double u=(jDay-2354755.0)/36525.0; // (1735-jan-0.5)
return -20.0 + 35.0*u*u;
}
// Implementation of algorithm by Schmadel & Zech (1979) for DeltaT computa
tion
double getDeltaTBySchmadelZech1979(const double jDay)
{
double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5)
double deltaT=(((((((((((-0.089491*u -0.117389)*u + 0.185489)*u + 0.
247433)*u - 0.159732)*u - 0.200097)*u + 0.075456)*u
+ 0.076929)*u - 0.020446)*u - 0.013867)*u + 0.003081
)*u + 0.001233)*u -0.000029;
return deltaT * 86400.0;
}
// Implementation of algorithm by Morrison & Stephenson (1982) for DeltaT c
omputation
double getDeltaTByMorrisonStephenson1982(const double jDay)
{
double u=(jDay-2382148.0)/36525.0; // (1810-jan-0.5)
return -15.0+32.50*std::pow(u,2);
}
// Implementation of algorithm by Stephenson & Morrison (1984) for DeltaT c
omputation
double getDeltaTByStephensonMorrison1984(const double jDay)
{
int year, month, day;
double deltaT = 0.;
getDateFromJulianDay(jDay, &year, &month, &day);
double yeardec=year+((month-1)*30.5+day/31*30.5)/366;
double u = (yeardec-1800)/100;
if (-391 < year and year <= 948)
deltaT = (44.3*u +320.0)*u +1360.0;
if (948 < year and year <= 1600)
deltaT = 25.5*u*u;
return deltaT;
}
// Implementation of algorithm by Stephenson & Morrison (1995) for DeltaT c
omputation
double getDeltaTByStephensonMorrison1995(const double jDay)
{
double u=(jDay-2385800.0)/36525.0; // (1820-jan-0.5)
return -20.0 + 31.0*u*u;
}
// Implementation of algorithm by Stephenson & Houlden (1986) for DeltaT co
mputation
double getDeltaTByStephensonHoulden(const double jDay)
{
int year, month, day;
double u;
double deltaT = 0.;
getDateFromJulianDay(jDay, &year, &month, &day);
double yeardec=year+((month-1)*30.5+day/31*30.5)/366;
if (year <= 948)
{
u = (yeardec-948)/100;
deltaT = (46.5*u -405.0)*u + 1830.0; // GZ: TODO: merge conf
lict, did I mess orig -405->+405? Or AW?
}
if (948 < year and year <= 1600)
{
u = (yeardec-1850)/100;
deltaT = 25.5*u*u;
}
return deltaT;
}
// Implementation of algorithm by Espenak (1987, 1989) for DeltaT computati
on
double getDeltaTByEspenak(const double jDay)
{
double u=(jDay-2451545.0)/36525.0; // (2000-jan-1.5)
return (64.3*u +61.0)*u +67.0;
}
// Implementation of algorithm by Borkowski (1988) for DeltaT computation
// This is explicitly compatible with ELP2000-85.
double getDeltaTByBorkowski(const double jDay)
{
double u=(jDay-2451545.0)/36525.0 + 3.75; // (2000-jan-1.5), deviati
on from 1625 as given in the paper.
return 40.0 + 35.0*u*u;
}
// Implementation of algorithm by Schmadel & Zech (1988) for DeltaT computa
tion
double getDeltaTBySchmadelZech1988(const double jDay)
{
double u=(jDay-2415020.0)/36525.0; // (1900-jan-0.5)
double deltaT = (((((((((((-0.058091*u -0.067471)*u +.145932)*u +.16
1416)*u -.149279)*u -.146960)*u +.079441)*u +.062971)*u -.022542)*u -.01246
2)*u +.003357)*u +.001148)*u-.000014;
return deltaT * 86400.0;
}
// Implementation of algorithm by Chapront-Touzé & Chapront (1991) for Delt
aT computation
double getDeltaTByChaprontTouze(const double jDay)
{
int year, month, day;
double deltaT = 0.;
getDateFromJulianDay(jDay, &year, &month, &day);
double u=(jDay-2451545.0)/36525.0; // (2000-jan-1.5)
if (-391 < year and year <= 948)
deltaT = (42.4*u -495.0)*u + 2177.0;
if (948 < year and year <= 1600)
deltaT = (23.6*u +100.0)*u + 102.0;
return deltaT;
}
// Implementation of algorithm by JPL Horizons for DeltaT computation
double getDeltaTByJPLHorizons(const double jDay)
{ // GZ: TODO: FIXME! It does not make sense to have zeros after 1620 in a
JPL Horizons compatible implementation!
int year, month, day;
double u;
double deltaT = 0.;
getDateFromJulianDay(jDay, &year, &month, &day);
if (-2999 < year and year < 948)
{
u=(jDay-2385800.0)/36525.0; // (1820-jan-1.5)
deltaT = 31.0*u*u;
}
if (948 < year and year <= 1620)
{
u=(jDay-2451545.0)/36525.0; // (2000-jan-1.5)
deltaT = (22.5*u +67.5)*u + 50.6;
}
return deltaT;
}
// Implementation of algorithm by Morrison & Stephenson (2004, 2005) for De
ltaT computation
double getDeltaTByMorrisonStephenson2004(const double jDay)
{
double u=(jDay-2385800.0)/36525.0; // (1820-jan-0.5)
return -20.0 + 32.0 * u*u;
}
// Implementation of algorithm by Reijs (2006) for DeltaT computation
double getDeltaTByReijs(const double jDay)
{
double OffSetYear = (2385800.0 - jDay)/365.25;
return ((1.8 * std::pow(OffSetYear,2)/200 + 1443*3.76/(2*M_PI)*(std:
:cos(2*M_PI*OffSetYear/1443)-1))*365.25)/1000;
}
// Implementation of algorithm by Chapront, Chapront-Touze & Francou (1997)
& Meeus (1998) for DeltaT computation
// 191 values in tenths of second for interpolation table 1620..2000, every
2nd year.
static const int MeeusDeltaTTable[] = { 1210, 1120, 1030, 950, 880, 820, 77
0, 720, 680, 630, 600, 560, 530, 510, 480,
460, 440, 420, 400, 380, 350, 330, 3
10, 290, 260, 240, 220, 200, 180, 160, 140, 120, 110, 100, 90, 80, 70,
70, 70, 70, // before 1700
70, 70, 80, 80, 90, 90, 90, 9
0, 90, 100, 100, 100, 100, 100, 100, 100, 100, 110, 110, 110, 110, 110, 12
0, 120, 120, // before 1750
120, 130, 130, 130, 140, 140, 140, 1
40, 150, 150, 150, 150, 150, 160, 160, 160, 160, 160, 160, 160, 160, 150, 1
50, 140, 130, // before 1800
131, 125, 122, 120, 120, 120, 120, 1
20, 120, 119, 116, 110, 102, 92, 82, 71, 62, 56, 54, 53, 54, 56,
59, 62, 65, // before 1850
68, 71, 73, 75, 76, 77, 73, 6
2, 52, 27, 14, -12, -28, -38, -48, -55, -53, -56, -57, -59, -60, -63, -6
5, -62, -47, // before 1900
-28, -1, 26, 53, 77, 104, 133, 1
60, 182, 202, 211, 224, 235, 238, 243, 240, 239, 239, 237, 240, 243, 253, 2
62, 273, 282, // before 1950
291, 300, 307, 314, 322, 331, 340, 3
50, 365, 383, 402, 422, 445, 465, 485, 505, 522, 538, 549, 558, 569, 583, 6
00, 616, 630, 650}; //closing: 2000
double getDeltaTByChaprontMeeus(const double jDay)
{
int year, month, day;
double deltaT = 0.;
getDateFromJulianDay(jDay, &year, &month, &day);
//const double u19=(jDay-2415020.0)/36525.0; // (1900-jan-0.5) Only
for approx form.
const double u20=(jDay-2451545.0)/36525.0; // (2000-jan-1.5)
if (year<948)
deltaT=(44.1*u20 +497.0)*u20+2177.0;
else if (year<1620)
deltaT=(25.3*u20 + 102.0)*u20+102.0;
// The next terms are the short approximative formulae. But Meeus gi
ves exact values, table, etc. for 1620..2000.
//else if (year<1900)
// deltaT= ((((((((( 123563.95*u19 + 727058.63)*u19 + 1818961
.41)*u19 + 2513807.78)*u19 + 2087298.89)*u19 + 1061660.75)*u19 +
// 324011.78)*u19 + 56282.84)*u19 + 5218.61)*u19
+ 228.95)*u19 - 2.50;
//else if (year<1997)
// deltaT= (((((((( 58353.42*u19 -232424.66)*u19 +372919.88)*
u19 - 303191.19)*u19 + 124906.15)*u19 - 18756.33)*u19 - 2637.80)*u19 + 815.
20)*u19 + 87.24)*u19 - 2.44;
else if (year <2000)
{
double yeardec=year+((month-1)*30.5+day/31*30.5)/366;
int pos=(year-1620)/2; // this is a deliberate integer divis
ion! 2->1, 3->1, 4->2, 5->2 etc.
deltaT= MeeusDeltaTTable[pos]+ (yeardec-(2*pos+1620))*0.5 *
(MeeusDeltaTTable[pos+1]-MeeusDeltaTTable[pos]);
deltaT /= 10.0;
}
else if (year<2100)
deltaT=(25.3*u20 + 102.0)*u20+102.0 + 0.37*(year-2100);
else // year > 2100
deltaT=(25.3*u20 + 102.0)*u20+102.0;
return deltaT;
}
// Implementation of algorithm by Montenbruck & Pfleger (2000) for DeltaT c
omputation
// Their implementation explicitly returns 0 for out-of-range data, so must
ours!
// Note: the polynomes do not form a well-behaved continuous line.
// The source of their data is likely the data table from Meeus, Astr. Alg.
1991.
double getDeltaTByMontenbruckPfleger(const double jDay)
{
double deltaT = 0.;
const double T=(jDay-2451545)/36525;
double t;
if (jDay<2387627.5) // 1825-01-01 0:00 ...
deltaT=0.0;
else if (jDay < 2396758.5) { // 1850-01-01 0:00
t=T+1.75;
deltaT=(( -572.3*t+413.9)*t -80.8)*t +10.4;
} else if (jDay < 2405889.5) { // 1875-01-01 0:00
t=T+1.50;
deltaT=(( 18.8*t-358.4)*t +46.3)*t + 6.6;
} else if (jDay < 2415020.5) { // 1900-01-01 0:00
t=T+1.25;
deltaT=(( 867.4*t-166.2)*t -10.8)*t - 3.9;
} else if (jDay < 2424151.5) { // 1925-01-01 0:00
t=T+1.00;
deltaT=((-1467.4*t+327.5)*t +114.1)*t - 2.6;
} else if (jDay < 2433282.5) { // 1950-01-01 0:00
t=T+0.75;
deltaT=(( 483.4*t - 8.2)*t - 6.3)*t +24.2;
} else if (jDay < 2442413.5) { // 1975-01-01 0:00
t=T+0.50;
deltaT=(( 550.7*t - 3.8)*t +32.5)*t +29.3;
} else if (jDay < 2453736.5) { // 2006-01-01 0:00
t=T+0.25;
deltaT=(( 1516.7*t-570.5)*t +130.5)*t +45.3;
}
return deltaT;
}
// Implementation of algorithm by Meeus & Simons (2000) for DeltaT computat
ion
// Zero outside defined range 1620..2000!
double getDeltaTByMeeusSimons(const double jDay)
{
int year, month, day;
double u;
double deltaT = 0.;
getDateFromJulianDay(jDay, &year, &month, &day);
//double yeardec=year+((month-1)*30.5+day/31*30.5)/366;
//double ub = (yeardec-2000)/100;
const double ub=(jDay-2451545.0)/36525.0; // (2000-jan-1.5)
if (year <1620)
deltaT=0.0;
else if (year < 1690)
{
u = 3.45 + ub;
//deltaT = +40.3 - 107.0*u + 50.0*std::pow(u,2) - 454.0*std:
:pow(u,3) + 1244.0*std::pow(u,4);
deltaT = (((1244.0*u -454.0)*u + 50.0)*u -107.0)*u +40.3;
}
else if (year < 1770)
{
u = 2.70 + ub;
//deltaT = +10.2 + 11.3*u - std::pow(u,2) - 16.0*std::pow(u,
3) + 70.0*std::pow(u,4);
deltaT = (((70.0*u -16.0)*u -1.0)*u +11.3)*u +10.2;
}
else if (year < 1820)
{
u = 2.05 + ub;
//deltaT = +14.7 - 18.8*u - 22.0*std::pow(u,2) + 173.0*std::
pow(u,3) + 6.0*std::pow(u,4);
deltaT = (((6.0*u -173.0)*u -22.0)*u -18.8)*u +14.7;
}
else if (year < 1870)
{
u = 1.55 + ub;
//deltaT = +5.7 + 12.7*u + 111.0*std::pow(u,2) - 534.0*std::
pow(u,3) + 1654.0*std::pow(u,4);
deltaT = (((1654.0*u -534.0)*u +111)*u +12.7)*u +5.7;
}
else if (year < 1900)
{
u = 1.15 + ub;
//deltaT = -5.8 - 14.6*u + 27.0*std::pow(u,2) + 101.0*std::p
ow(u,3) + 8234.0*std::pow(u,4);
deltaT = (((8234.0*u +101.0)*u +27.0)*u - 14.6)*u -5.8;
}
else if (year < 1940)
{
u = 0.80 + ub;
//deltaT = +21.4 + 67.0*u + 443.0*std::pow(u,2) + 19.0*std::
pow(u,3) + 4441.0*std::pow(u,4);
deltaT = (((4441.0*u + 19.0)*u +443.0)*u +67.0)*u +21.4;
}
else if (year < 1990)
{
u = 0.35 + ub;
//deltaT = +36.2 + 74.0*u + 189.0*std::pow(u,2) - 140.0*std:
:pow(u,3) - 1883.0*std::pow(u,4);
deltaT = (((-1883.0*u -140.0)*u +189.0)*u +74.0)*u +36.2;
}
else if (year <= 2000)
{
u = 0.05 + ub;
//deltaT = +60.8 + 82.0*u + 188.0*std::pow(u,2) - 5034.0*std
::pow(u,3);
deltaT = ((-5034.0*u +188.0)*u +82.0)*u +60.8;
}
return deltaT;
}
// Implementation of algorithm by Reingold & Dershowitz (Cal. Calc. 1997, 2
001, 2007, Cal. Tab. 2002) for DeltaT computation.
// GZ: Created as yet another multi-segment polynomial fit through the tabl
e in Meeus: Astronomical Algorithms (1991).
// GZ: Note that only the Third edition (2007) adds the 1700-1799 term.
// GZ: More efficient reimplementation with stricter adherence to the sourc
e.
double getDeltaTByReingoldDershowitz(const double jDay)
{
int year, month, day;
getDateFromJulianDay(jDay, &year, &month, &day);
// GZ: R&D don't use a float-fraction year, but explicitly only the
integer year! And R&D use a proleptic Gregorian year before 1582.
// GZ: We cannot do that, but the difference is negligible.
// GZ: FIXME: why are displayed values so far off the computed value
s? It seems currently broken!
double deltaT=0.0; // If it returns 0, there is a bug!
if ((year >= 2019) || (year < 1620))
{
double jdYear_0; getJDFromDate(&jdYear_0, year, 1, 1, 0, 0,
0);
double jd1810_0; getJDFromDate(&jd1810_0, 1810, 1, 1, 0, 0,
0);
double x = (jdYear_0-jd1810_0+0.5);
deltaT = x*x/41048480.0 - 15.0;
}
else if (year >= 1988)
{
deltaT = year-1933.0;
}
else if (year >= 1800)
{
double jd1900_0; getJDFromDate(&jd1900_0, 1900, 1, 1, 0, 0,
0);
double jdYear_5; getJDFromDate(&jdYear_5, year, 7, 1, 0, 0,
0);
double c = (jdYear_5-jd1900_0) * (1/36525.0);
if (year >= 1900)
{
deltaT = (((((((-0.212591*c +0.677066)*c -0.861938)*
c +0.553040)*c -0.181133)*c +0.025184)*c +0.000297)*c -0.00002) * 86400.0;
}
else //if (year >= 1800)
{
deltaT = ((((((((((2.043794*c +11.636204)*c +28.3162
89)*c +38.291999)*c +31.332267)*c +15.845535)*c +4.867575)*c +0.865736)*c +
0.083563)*c +0.003844)*c -0.000009) * 86400.0;
}
}
else if (year >= 1700)
{ // This term was added in the third edition (2007), its omission w
as a fault of the authors!
double yDiff1700 = year-1700.0;
deltaT = ((-0.0000266484*yDiff1700 +0.003336121)*yDiff1700 +
0.005092142)*yDiff1700 + 8.118780842;
}
else if (year >= 1620)
{
double yDiff1600 = year-1600.0;
deltaT = (0.0219167*yDiff1600 -4.0675)*yDiff1600 +196.58333;
}
return deltaT;
} }
double getMoonSecularAcceleration(double jDay) double getMoonSecularAcceleration(const double jDay, const double nd)
{ {
// Method described is here: http://eclipse.gsfc.nasa.gov/SEcat5/sec ular.html
int year, month, day; int year, month, day;
getDateFromJulianDay(jDay, &year, &month, &day); getDateFromJulianDay(jDay, &year, &month, &day);
double t = (year-1995)/100;
return -0.12932224 * t * t; double yeardec=year+((month-1)*30.5+day/31*30.5)/366.0;
double t = (yeardec-1955.5)/100.0;
// n.dot for secular acceleration of the Moon in ELP2000-82B
// have value -23.8946 "/cy/cy
return -0.91072 * (-23.8946 + std::abs(nd))*t*t;
}
double getDeltaTStandardError(const double jDay)
{
int year, month, day;
getDateFromJulianDay(jDay, &year, &month, &day);
//double yeardec=year+((month-1)*30.5+day/31*30.5)/366;
double sigma = -1.;
if (-1000 <= year and year <= 1600)
{
double cDiff1820= (jDay-2385800.0)/36525.0; // 1820.0=182
0-jan-0.5=2385800.0
// sigma = std::pow((yeardec-1820.0)/100,2); // sigma(DeltaT
) = 0.8*u^2
sigma = 0.8 * cDiff1820 * cDiff1820;
}
return sigma;
} }
} // end of the StelUtils namespace } // end of the StelUtils namespace
 End of changes. 50 change blocks. 
91 lines changed or deleted 610 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/