Commit 152179b5 authored by Pablo de Vicente's avatar Pablo de Vicente

Fixed a bug first pointed out by Mikhail Zotov in a private communication,

who noticed that precession between B1950 to J2000 as done by KStars does not
coincide with the results given by SIMBAD (http://simbad.u-strasbg.fr/sim-fidl.pl).  M. Zotov has not opened yet a bug in bugs.kde.org, so I have no bug
number to refer to in this commit.

  The bug comes from the fact that conversion between coordinates in B1950
to J2000 (and viceversa) involves changing from (old) catalog FK4 to FK5,
and a simple precesion, as done by KStars is not enough.

To fix this bug I have used the following reference:
Smith, C. A.; Kaplan, G. H.; Hughes, J. A.; Seidelmann, P. K.; Yallop,
B. D.; Hohenkerk, C. Y.
Astronomical Journal, vol. 97, Jan. 1989, p. 265-279

The conversion between the FK4 catalog and the FK5 catalog requires 4 steps:
 - Drop E-Terms in B1950 coordinates. In the past, the mean places of
   stars published in the FK4 catalog included the contribution to the
   aberration due to the ellipticity of the orbit of the Earth. These terms,
   known as E-terms were almost constant, and in the newer FK5 catalog they
   are not included.
 - Precess from B1950 to 1984, January 1st, at 0h, using the parameters
   given by the Astronomische Rechen-Institut.
 - Apply the zero-point correction to the right ascensions to correct for the
   equinox error of the FK4. This is done for 1984, Jan 1st, at 0h.
 - Precess from 1984 Jan 1st at 0h, to J2000 using the new precessional
   parameters.

   This bug may seem not important since KStars produces a result which is
almost correct, but since the conversion between B1950 and J2000 is done
very frequently among astronomers it needs to be corrected. This bug affects
the precessional module and the galactic/equatorial conversion module. This
fix ONLY applies to the first module, although now it is straightforward to
correct the second module. I will do this later.

  Since Stephan Kulow only accepts big bug fixes in the KDE 3.2 deep freeze
period in order to ensure a fast stability of the first relase, I commit
these fixes only to the KDE 3.1.4 branch, even knowing there may be no ç
KDE 3.1.5 release in the future. As soon as the KDE 3.2 branch is opened
for small bug fixes I will port this fix.

CCMAIL: kstars-devel@kde.org

svn path=/branches/KDE_3_1_BRANCH/kdeedu/kstars/; revision=269012
parent 4e1155f8
......@@ -30,6 +30,9 @@ void KSNumbers::updateValues( long double jd ) {
//Julian Centuries since J2000.0
T = ( jd - J2000 ) / 36525.;
//Julian Centuries since B1950.0
TB = ( jd - 2433282.4235 ) / 36524.2199;
// Julian Millenia since J2000.0
jm = T / 10.0;
......@@ -115,4 +118,39 @@ void KSNumbers::updateValues( long double jd ) {
P2[0][2] = CX*SY;
P2[1][2] = -1.0*SX*SY;
P2[2][2] = CY;
// eqnCorr.setD ( (0.0775 + 0.0850 * TB)/3600. );
//Compute Precession Matrices from B1950 to 1984 using Newcomb formulae
XB.setD( 0.217697 );
YB.setD( 0.189274 );
ZB.setD( 0.217722 );
XB.SinCos( SXB, CXB );
YB.SinCos( SYB, CYB );
ZB.SinCos( SZB, CZB );
//P1B is used to precess from 1984 to B1950:
P1B[0][0] = CXB*CYB*CZB - SXB*SZB;
P1B[1][0] = CXB*CYB*SZB + SXB*CZB;
P1B[2][0] = CXB*SYB;
P1B[0][1] = -1.0*SXB*CYB*CZB - CXB*SZB;
P1B[1][1] = -1.0*SXB*CYB*SZB + CXB*CZB;
P1B[2][1] = -1.0*SXB*SYB;
P1B[0][2] = -1.0*SYB*CZB;
P1B[1][2] = -1.0*SYB*SZB;
P1B[2][2] = CYB;
//P2 is used to precess from B1950 to 1984 (it is the transpose of P1)
P2B[0][0] = CXB*CYB*CZB - SXB*SZB;
P2B[1][0] = -1.0*SXB*CYB*CZB - CXB*SZB;
P2B[2][0] = -1.0*SYB*CZB;
P2B[0][1] = CXB*CYB*SZB + SXB*CZB;
P2B[1][1] = -1.0*SXB*CYB*SZB + CXB*CZB;
P2B[2][1] = -1.0*SYB*SZB;
P2B[0][2] = CXB*SYB;
P2B[1][2] = -1.0*SXB*SYB;
P2B[2][2] = CYB;
}
......@@ -85,19 +85,26 @@ public:
/**@returns element of P1 precession array at position [i1][i2] */
double p1( int i1, int i2 ) const { return P1[i1][i2]; }
/**@returns element of P1 precession array at position [i1][i2] */
double p1b( int i1, int i2 ) const { return P1B[i1][i2]; }
/**@returns element of P2 precession array at position [i1][i2] */
double p2( int i1, int i2 ) const { return P2[i1][i2]; }
/**@returns element of P2 precession array at position [i1][i2] */
double p2b( int i1, int i2 ) const { return P2B[i1][i2]; }
/**@short update all values for the date given as an argument. */
void updateValues( long double jd );
private:
dms Obliquity, K, L, L0, LM, M, M0, O, P;
dms XP, YP, ZP;
dms XP, YP, ZP, XB, YB, ZB;
double CX, SX, CY, SY, CZ, SZ;
double P1[3][3], P2[3][3];
double CXB, SXB, CYB, SYB, CZB, SZB;
double P1[3][3], P2[3][3], P1B[3][3], P2B[3][3];
double deltaObliquity, deltaEcLong;
double e, T;
double e, T, TB;
long double days;
double jm;
......
......@@ -14,7 +14,6 @@
* (at your option) any later version. *
* *
***************************************************************************/
#include "modcalcprec.h"
#include "modcalcprec.moc"
#include "dms.h"
......@@ -47,7 +46,7 @@ modCalcPrec::modCalcPrec(QWidget *parentSplit, const char *name) : QWidget(paren
InputBox->setTitle( i18n("Original Coordinates") );
QHBox * buttonBox = new QHBox(InputBox);
QPushButton * Compute = new QPushButton( i18n( "Compute" ), buttonBox);
QPushButton * Clear = new QPushButton( i18n( "Clear" ), buttonBox );
......@@ -112,7 +111,7 @@ modCalcPrec::modCalcPrec(QWidget *parentSplit, const char *name) : QWidget(paren
epochfBox->setMargin(6);
epochfBox->setSpacing(6);
D1Lay->setMargin(14);
D1Lay->addWidget(rafdecfBox);
......@@ -167,7 +166,7 @@ double modCalcPrec::setCurrentEpoch () {
}
double modCalcPrec::getEpoch (QString eName) {
double epoch = eName.toDouble();
return epoch;
......@@ -187,7 +186,7 @@ void modCalcPrec::slotClearCoords (void) {
void modCalcPrec::slotComputeCoords (void) {
SkyPoint sp;
sp = getEquCoords();
double epoch0 = getEpoch( epoch0Name->text() );
double epochf = getEpoch( epochfName->text() );
......@@ -200,7 +199,7 @@ void modCalcPrec::showEquCoords ( SkyPoint sp ) {
rafBox->show( sp.ra() );
decfBox->show( sp.dec() );
}
double modCalcPrec::JdtoEpoch (long double jd) {
......@@ -223,74 +222,210 @@ long double modCalcPrec::epochToJd (double epoch) {
}
SkyPoint modCalcPrec::precess (dms RA0, dms Dec0, double epoch0, double epochf) {
SkyPoint modCalcPrec::B1950ToJ2000 (dms RA0, dms Dec0) {
double cosRA0, sinRA0, cosDec0, sinDec0;
dms RA, Dec;
double v[3], s[3];
double v[3], s[3];
long double jd0 = epochToJd ( epoch0 );
long double jdf = epochToJd ( epochf );
// 1984 January 1 0h
KSNumbers *num = new KSNumbers (2445700.5);
RA0.SinCos( sinRA0, cosRA0 );
Dec0.SinCos( sinDec0, cosDec0 );
// Eterms due to aberration
SkyPoint sPoint = SkyPoint (RA0, Dec0);
sPoint.addEterms();
sPoint.ra().SinCos( sinRA0, cosRA0 );
sPoint.dec().SinCos( sinDec0, cosDec0 );
//Need to first precess to J2000.0 coords
// Precession from B1950 to J1984
s[0] = cosRA0*cosDec0;
s[1] = sinRA0*cosDec0;
s[2] = sinDec0;
for ( unsigned int i=0; i<3; ++i ) {
v[i] = num->p2b( 0, i )*s[0] + num->p2b( 1, i )*s[1] +
num->p2b( 2, i )*s[2];
}
if ( jd0 != J2000 ) {
//v is a column vector representing input coordinates.
// RA zero-point correction at 1984 day 1, 0h.
RA.setRadians( atan2( v[1],v[0] ) );
Dec.setRadians( asin( v[2] ) );
KSNumbers *num = new KSNumbers (jd0);
RA.setH( RA.Hours() + 0.06390/3600. );
RA.SinCos( sinRA0, cosRA0 );
Dec.SinCos( sinDec0, cosDec0 );
v[0] = cosRA0*cosDec0;
v[1] = sinRA0*cosDec0;
v[2] = sinDec0;
s[0] = cosRA0*cosDec0;
s[1] = sinRA0*cosDec0;
s[2] = sinDec0;
//s is the product of P1 and v; s represents the
//coordinates precessed to J2000
for ( unsigned int i=0; i<3; ++i ) {
s[i] = num->p1( 0, i )*v[0] + num->p1( 1, i )*v[1] +
num->p1( 2, i )*v[2];
}
delete num;
// Precession from 1984 to J2000.
} else {
//Input coords already in J2000, set s accordingly.
s[0] = cosRA0*cosDec0;
s[1] = sinRA0*cosDec0;
s[2] = sinDec0;
for ( unsigned int i=0; i<3; ++i ) {
v[i] = num->p1( 0, i )*s[0] +
num->p1( 1, i )*s[1] +
num->p1( 2, i )*s[2];
}
KSNumbers *num = new KSNumbers (jdf);
//Multiply P2 and s to get v, the vector representing the new coords.
RA.setRadians( atan2( v[1],v[0] ) );
Dec.setRadians( asin( v[2] ) );
delete num;
sPoint.setRA(RA);
sPoint.setDec(Dec);
return sPoint;
}
SkyPoint modCalcPrec::J2000ToB1950 (dms RA0, dms Dec0) {
double cosRA0, sinRA0, cosDec0, sinDec0;
dms RA, Dec;
double v[3], s[3];
// 1984 January 1 0h
KSNumbers *num = new KSNumbers (2445700.5);
RA0.SinCos(sinRA0,cosRA0);
Dec0.SinCos(sinDec0,cosDec0);
s[0] = cosRA0*cosDec0;
s[1] = sinRA0*cosDec0;
s[2] = sinDec0;
// Precession from 1984 to J2000.
for ( unsigned int i=0; i<3; ++i ) {
v[i] = num->p2( 0, i )*s[0] + num->p2( 1, i )*s[1] +
v[i] = num->p2( 0, i )*s[0] +
num->p2( 1, i )*s[1] +
num->p2( 2, i )*s[2];
}
delete num;
//Extract RA, Dec from the vector:
// RA.setRadians( atan( v[1]/v[0] ) );
RA.setRadians( atan2( v[1],v[0] ) );
Dec.setRadians( asin( v[2] ) );
//resolve ambiguity of atan()
// if ( v[0] < 0.0 ) {
// RA.setD( RA.Degrees() + 180.0 );
// } else if( v[1] < 0.0 ) {
// RA.setD( RA.Degrees() + 360.0 );
// }
// RA zero-point correction at 1984 day 1, 0h.
RA.setH( RA.Hours() - 0.06390/3600. );
RA.SinCos( sinRA0, cosRA0 );
Dec.SinCos( sinDec0, cosDec0 );
// Precession from B1950 to J1984
s[0] = cosRA0*cosDec0;
s[1] = sinRA0*cosDec0;
s[2] = sinDec0;
for ( unsigned int i=0; i<3; ++i ) {
v[i] = num->p1b( 0, i )*s[0] + num->p1b( 1, i )*s[1] +
num->p1b( 2, i )*s[2];
}
if (RA.Degrees() < 0.0 )
RA.setD( RA.Degrees() + 360.0 );
RA.setRadians( atan2( v[1],v[0] ) );
Dec.setRadians( asin( v[2] ) );
// Eterms due to aberration
SkyPoint sPoint = SkyPoint (RA, Dec);
sPoint.subtractEterms();
delete num;
return sPoint;
}
SkyPoint modCalcPrec::precess (dms RA0, dms Dec0, double epoch0, double epochf) {
double cosRA0, sinRA0, cosDec0, sinDec0;
dms RA, Dec;
double v[3], s[3];
long double jd0 = epochToJd ( epoch0 );
long double jdf = epochToJd ( epochf );
RA0.SinCos( sinRA0, cosRA0 );
Dec0.SinCos( sinDec0, cosDec0 );
if ( jd0 == B1950) {
SkyPoint sp = SkyPoint(RA0, Dec0);
sp.B1950ToJ2000();
RA0.setD( sp.ra().Degrees() );
Dec0.setD( sp.dec().Degrees() );
jd0 = J2000;
}
if (jd0 !=jdf) {
// The original coordinate is referred to the FK5 system and
// is NOT J2000.
if ( jd0 != J2000 ) {
//v is a column vector representing input coordinates.
KSNumbers *num = new KSNumbers (jd0);
v[0] = cosRA0*cosDec0;
v[1] = sinRA0*cosDec0;
v[2] = sinDec0;
//Need to first precess to J2000.0 coords
//s is the product of P1 and v; s represents the
//coordinates precessed to J2000
for ( unsigned int i=0; i<3; ++i ) {
s[i] = num->p1( 0, i )*v[0] +
num->p1( 1, i )*v[1] +
num->p1( 2, i )*v[2];
}
delete num;
//Input coords already in J2000, set s accordingly.
} else {
s[0] = cosRA0*cosDec0;
s[1] = sinRA0*cosDec0;
s[2] = sinDec0;
}
if ( jdf == B1950) {
RA.setRadians( atan2( s[1],s[0] ) );
Dec.setRadians( asin( s[2] ) );
SkyPoint sp = SkyPoint(RA, Dec);
sp.J2000ToB1950();
return sp;
}
KSNumbers *num = new KSNumbers (jdf);
//Multiply P2 and s to get v, the vector representing i
//the new coords.
for ( unsigned int i=0; i<3; ++i ) {
v[i] = num->p2( 0, i )*s[0] +
num->p2( 1, i )*s[1] +
num->p2( 2, i )*s[2];
}
delete num;
//Extract RA, Dec from the vector:
// RA.setRadians( atan( v[1]/v[0] ) );
RA.setRadians( atan2( v[1],v[0] ) );
Dec.setRadians( asin( v[2] ) );
//resolve ambiguity of atan()
// if ( v[0] < 0.0 ) {
// RA.setD( RA.Degrees() + 180.0 );
// } else if( v[1] < 0.0 ) {
// RA.setD( RA.Degrees() + 360.0 );
// }
if (RA.Degrees() < 0.0 )
RA.setD( RA.Degrees() + 360.0 );
SkyPoint sPoint = SkyPoint (RA, Dec);
return sPoint;
} else {
SkyPoint sPoint = SkyPoint (RA0, Dec0);
return sPoint;
}
}
......@@ -48,6 +48,8 @@ public:
double JdtoEpoch (long double jd);
long double epochToJd (double epoch);
SkyPoint precess (dms ra0, dms dec0, double e0, double ef);
SkyPoint B1950ToJ2000 (dms RA0, dms Dec0);
SkyPoint J2000ToB1950 (dms RA0, dms Dec0);
public slots:
void slotClearCoords (void);
......
......@@ -213,12 +213,143 @@ void SkyPoint::updateCoords( KSNumbers *num, bool includePlanets ) {
//Correct the catalog coordinates for the time-dependent effects
//of precession, nutation and aberration
precess(num);
nutate(num);
aberrate(num);
}
void SkyPoint::B1950ToJ2000(void) {
double cosRA, sinRA, cosDec, sinDec;
double v[3], s[3];
// 1984 January 1 0h
KSNumbers *num = new KSNumbers (2445700.5);
// Eterms due to aberration
addEterms();
RA.SinCos( sinRA, cosRA );
Dec.SinCos( sinDec, cosDec );
// Precession from B1950 to J1984
s[0] = cosRA*cosDec;
s[1] = sinRA*cosDec;
s[2] = sinDec;
for ( unsigned int i=0; i<3; ++i ) {
v[i] = num->p2b( 0, i )*s[0] + num->p2b( 1, i )*s[1] +
num->p2b( 2, i )*s[2];
}
// RA zero-point correction at 1984 day 1, 0h.
RA.setRadians( atan2( v[1],v[0] ) );
Dec.setRadians( asin( v[2] ) );
RA.setH( RA.Hours() + 0.06390/3600. );
RA.SinCos( sinRA, cosRA );
Dec.SinCos( sinDec, cosDec );
s[0] = cosRA*cosDec;
s[1] = sinRA*cosDec;
s[2] = sinDec;
// Precession from 1984 to J2000.
for ( unsigned int i=0; i<3; ++i ) {
v[i] = num->p1( 0, i )*s[0] +
num->p1( 1, i )*s[1] +
num->p1( 2, i )*s[2];
}
RA.setRadians( atan2( v[1],v[0] ) );
Dec.setRadians( asin( v[2] ) );
delete num;
}
void SkyPoint::J2000ToB1950(void) {
double cosRA, sinRA, cosDec, sinDec;
double cosRA0, sinRA0, cosDec0, sinDec0;
double v[3], s[3];
// 1984 January 1 0h
KSNumbers *num = new KSNumbers (2445700.5);
RA0.SinCos( sinRA0, cosRA0 );
Dec0.SinCos( sinDec0, cosDec0 );
s[0] = cosRA0*cosDec0;
s[1] = sinRA0*cosDec0;
s[2] = sinDec0;
// Precession from J2000 to 1984 day, 0h.
for ( unsigned int i=0; i<3; ++i ) {
v[i] = num->p2( 0, i )*s[0] +
num->p2( 1, i )*s[1] +
num->p2( 2, i )*s[2];
}
RA.setRadians( atan2( v[1],v[0] ) );
Dec.setRadians( asin( v[2] ) );
// RA zero-point correction at 1984 day 1, 0h.
RA.setH( RA.Hours() - 0.06390/3600. );
RA.SinCos( sinRA, cosRA );
Dec.SinCos( sinDec, cosDec );
// Precession from B1950 to J1984
s[0] = cosRA*cosDec;
s[1] = sinRA*cosDec;
s[2] = sinDec;
for ( unsigned int i=0; i<3; ++i ) {
v[i] = num->p1b( 0, i )*s[0] + num->p1b( 1, i )*s[1] +
num->p1b( 2, i )*s[2];
}
RA.setRadians( atan2( v[1],v[0] ) );
Dec.setRadians( asin( v[2] ) );
// Eterms due to aberration
subtractEterms();
delete num;
}
SkyPoint SkyPoint::Eterms(void) {
double sd, cd, sinEterm, cosEterm;
dms raTemp, raDelta, decDelta;
Dec.SinCos(sd,cd);
raTemp.setH( RA.Hours() + 11.25);
raTemp.SinCos(sinEterm,cosEterm);
raDelta.setH( 0.0227*sinEterm/(3600.*cd) );
decDelta.setD( 0.341*cosEterm*sd/3600. + 0.029*cd/3600. );
SkyPoint spDelta = SkyPoint (raDelta, decDelta);
return spDelta;
}
void SkyPoint::addEterms(void) {
SkyPoint spd = Eterms();
RA.setD( RA.Degrees() + spd.ra().Degrees() );
Dec.setD( Dec.Degrees() + spd.dec().Degrees() );
}
void SkyPoint::subtractEterms(void) {
SkyPoint spd = Eterms();
RA.setD( RA.Degrees() - spd.ra().Degrees() );
Dec.setD( Dec.Degrees() - spd.dec().Degrees() );
}
......@@ -26,6 +26,8 @@
#include "ksnumbers.h"
#include "dms.h"
#define B1950 2433282.4235
/**
*Encapsulates the sky coordinates of a point in the sky. The
......@@ -243,6 +245,57 @@ public:
/**Determine the effects of aberration for this SkyPoint*/
void aberrate(const KSNumbers *num);
/** Determine the E-terms of aberration
*In the past, the mean places of stars published in catalogs included
*the contribution to the aberration due to the ellipticity of the orbit
*of the Earth. These terms, known as E-terms were almost constant, and
*in the newer catalogs (FK5) are not included. Therefore to convert from
*FK4 to FK5 one has to compute these E-terms.
*/
SkyPoint Eterms(void);
/** Exact precession from Besselian epoch 1950 to epoch J2000. The coordinates
*referred to the first epoch are in the FK4 catalog, while the latter are in the
*Fk5 one.
*Reference: Smith, C. A.; Kaplan, G. H.; Hughes, J. A.; Seidelmann,
*P. K.; Yallop, B. D.; Hohenkerk, C. Y.
*Astronomical Journal, vol. 97, Jan. 1989, p. 265-279
*This transformation requires 4 steps:
* - Correct E-terms
* - Precess from B1950 to 1984, January 1st, 0h, using Newcomb expressions
* - Add zero point correction in right ascension for 1984
* - Precess from 1984, January 1st, 0h to J2000
*/
void B1950ToJ2000(void);
/** Exact precession from epoch J2000 Besselian epoch 1950. The coordinates
*referred to the first epoch are in the FK4 catalog, while the latter are in the
*Fk5 one.
*Reference: Smith, C. A.; Kaplan, G. H.; Hughes, J. A.; Seidelmann,
*P. K.; Yallop, B. D.; Hohenkerk, C. Y.
*Astronomical Journal, vol. 97, Jan. 1989, p. 265-279
*This transformation requires 4 steps:
* - Precess from J2000 to 1984, January 1st, 0h
* - Add zero point correction in right ascension for 1984
* - Precess from 1984, January 1st, 0h, to B1950 using Newcomb expressions
* - Correct E-terms
*/
void J2000ToB1950(void);
/** Coordinates in the FK4 catalog include the effect of aberration due
*to the ellipticity of the orbit of the Earth. Coordinates in the FK5
*catalog do not include these terms. In order to convert from B1950 (FK4)
*to actual mean places one has to use this function.
*/
void addEterms(void);
/** Coordinates in the FK4 catalog include the effect of aberration due
*to the ellipticity of the orbit of the Earth. Coordinates in the FK5
*catalog do not include these terms. In order to convert from FK5 coordinates
* to B1950 (FK4) one has to use this function.
*/
void subtractEterms(void);
protected:
/**Precess this SkyPoint's catalog coordinates */