summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAkarsh Simha <akarsh@kde.org>2016-09-28 07:33:50 (GMT)
committerAkarsh Simha <akarsh@kde.org>2016-09-29 02:20:27 (GMT)
commit66a7ea2f3829bff896128836479bdf5fb365df48 (patch)
treec03d64c3b253611d69b80b9dfe9814dc29f47a84
parent743cba78b14e7ec002484faa4e3e124deccf782b (diff)
Experiment: Explicitly write out the computations of precession
-rw-r--r--kstars/ksnumbers.cpp79
-rw-r--r--kstars/ksnumbers.h30
-rw-r--r--kstars/skyobjects/skypoint.cpp27
3 files changed, 66 insertions, 70 deletions
diff --git a/kstars/ksnumbers.cpp b/kstars/ksnumbers.cpp
index 5193472..b1e0b44 100644
--- a/kstars/ksnumbers.cpp
+++ b/kstars/ksnumbers.cpp
@@ -183,27 +183,27 @@ void KSNumbers::computeConstantValues() {
//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;
+ 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)
// FIXME: This can be optimized by taking the transpose of P1 instead of recomputing it from scratch
- 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;
+ 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;
}
KSNumbers::~KSNumbers(){
@@ -313,30 +313,29 @@ void KSNumbers::updateValues( long double jd ) {
ZP.SinCos( SZ, CZ );
//P1 is used to precess from any epoch to J2000
- // Note: P1 is a rotation matrix, and P2 is its transpose = inverse (also a rotation matrix)
- P1(0, 0) = CX*CY*CZ - SX*SZ;
- P1(1, 0) = CX*CY*SZ + SX*CZ;
- P1(2, 0) = CX*SY;
- P1(0, 1) = -1.0*SX*CY*CZ - CX*SZ;
- P1(1, 1) = -1.0*SX*CY*SZ + CX*CZ;
- P1(2, 1) = -1.0*SX*SY;
- P1(0, 2) = -1.0*SY*CZ;
- P1(1, 2) = -1.0*SY*SZ;
- P1(2, 2) = CY;
+ // FIXME: Is this a rotation matrix? If so, quaternions might be more efficient
+ // A: Yes, it has to be, because the inverse is the transpose, so the matrix is orthogonal 3x3
+ P1[0][0] = CX*CY*CZ - SX*SZ;
+ P1[1][0] = CX*CY*SZ + SX*CZ;
+ P1[2][0] = CX*SY;
+ P1[0][1] = -1.0*SX*CY*CZ - CX*SZ;
+ P1[1][1] = -1.0*SX*CY*SZ + CX*CZ;
+ P1[2][1] = -1.0*SX*SY;
+ P1[0][2] = -1.0*SY*CZ;
+ P1[1][2] = -1.0*SY*SZ;
+ P1[2][2] = CY;
//P2 is used to precess from J2000 to any other epoch (it is the transpose of P1)
- // FIXME: More optimization -- just use P1(j, i) instead of P2(i, j) in code
- P2(0, 0) = P1(0, 0);
- P2(1, 0) = P1(0, 1);
- P2(2, 0) = P1(0, 2);
- P2(0, 1) = P1(1, 0);
- P2(1, 1) = P1(1, 1);
- P2(2, 1) = P1(1, 2);
- P2(0, 2) = P1(2, 0);
- P2(1, 2) = P1(2, 1);
- P2(2, 2) = P1(2, 2);
-
-
+ // FIXME: More optimization -- just use P1[j][i] instead of P2[i][j] in code
+ P2[0][0] = P1[0][0];
+ P2[1][0] = P1[0][1];
+ P2[2][0] = P1[0][2];
+ P2[0][1] = P1[1][0];
+ P2[1][1] = P1[1][1];
+ P2[2][1] = P1[1][2];
+ P2[0][2] = P1[2][0];
+ P2[1][2] = P1[2][1];
+ P2[2][2] = P1[2][2];
// Mean longitudes for the planets. radians
//
diff --git a/kstars/ksnumbers.h b/kstars/ksnumbers.h
index aba30f4..2ce8a62 100644
--- a/kstars/ksnumbers.h
+++ b/kstars/ksnumbers.h
@@ -24,10 +24,6 @@
#include "cachingdms.h"
-#include <Eigen/Dense>
-
-using Eigen::Matrix3d;
-
/** @class KSNumbers
*
*There are several time-dependent values used in position calculations,
@@ -99,23 +95,21 @@ public:
/** @return Julian Millenia since J2000*/
inline double julianMillenia() const { return jm; }
- /** @return element of P1 precession array at position (i1, i2) */
- inline double p1( int i1, int i2 ) const { return P1(i1, i2); }
+ /** @return element of P1 precession array at position [i1][i2] */
+ inline double p1( int i1, int i2 ) const { return P1[i1][i2]; }
+
+ /** @return element of P2 precession array at position [i1][i2] */
+ inline double p2( int i1, int i2 ) const { return P2[i1][i2]; }
- /** @return element of P2 precession array at position (i1, i2) */
- inline double p2( int i1, int i2 ) const { return P2(i1, i2); }
+ /** @return element of P1B precession array at position [i1][i2] */
+ inline double p1b( int i1, int i2 ) const { return P1B[i1][i2]; }
- /** @return element of P1B precession array at position (i1, i2) */
- inline double p1b( int i1, int i2 ) const { return P1B(i1, i2); }
+ /** @return element of P2B precession array at position [i1][i2] */
+ inline double p2b( int i1, int i2 ) const { return P2B[i1][i2]; }
- /** @return element of P2B precession array at position (i1, i2) */
- inline double p2b( int i1, int i2 ) const { return P2B(i1, i2); }
+ inline const double *p1() const { return (*P1); }
- /** @return the precession matrix directly **/
- inline const Matrix3d &p2() const { return P1; }
- inline const Matrix3d &p1() const { return P2; }
- inline const Matrix3d &p1b() const { return P1B; }
- inline const Matrix3d &p2b() const { return P2B; }
+ inline const double *p2() const { return (*P2); }
/**
*@short compute constant values that need to be computed only once per instance of the application
@@ -140,7 +134,7 @@ private:
dms XP, YP, ZP, XB, YB, ZB;
double CX, SX, CY, SY, CZ, SZ;
double CXB, SXB, CYB, SYB, CZB, SZB;
- Matrix3d P1, P2, P1B, P2B;
+ double P1[3][3], P2[3][3], P1B[3][3], P2B[3][3];
double deltaObliquity, deltaEcLong;
double e, T;
long double days; // JD for which the last update was called
diff --git a/kstars/skyobjects/skypoint.cpp b/kstars/skyobjects/skypoint.cpp
index 6f9ae9a..c7f86af 100644
--- a/kstars/skyobjects/skypoint.cpp
+++ b/kstars/skyobjects/skypoint.cpp
@@ -190,8 +190,13 @@ void SkyPoint::setFromEcliptic( const CachingDms *Obliquity, const dms& EcLong,
void SkyPoint::precess( const KSNumbers *num ) {
double cosRA0, sinRA0, cosDec0, sinDec0;
- const Eigen::Matrix3d &precessionMatrix = num->p2();
- Eigen::Vector3d v, s;
+ double s[3];
+ const double *P = num->p1();
+
+ // NOTE: Since dms::sin() and dms::cos() are inline methods, one
+ // might be tempted to replace cosRA0 -> ra0().cos() etc and avoid
+ // dms::SinCos. callgrind tells me that this deteriorates performance.
+ // -- asimha
RA0.SinCos( sinRA0, cosRA0 );
Dec0.SinCos( sinDec0, cosDec0 );
@@ -208,18 +213,16 @@ void SkyPoint::precess( const KSNumbers *num ) {
// there isn't much overhead in accessing them.
//Multiply P2 and s to get v, the vector representing the new coords.
- // for ( unsigned int i=0; i<3; ++i ) {
- // v[i] = 0.0;
- // for (uint j=0; j< 3; ++j) {
- // v[i] += num->p2( j, i )*s[j];
- // }
- // }
- v.noalias() = precessionMatrix * s;
-
+#define v0 P[0] * s[0] + P[1] * s[1] + P[2] * s[2]
+#define v1 P[3] * s[0] + P[4] * s[1] + P[5] * s[2]
+#define v2 P[6] * s[0] + P[7] * s[1] + P[8] * s[2]
//Extract RA, Dec from the vector:
- RA.setUsing_atan2( v[1], v[0] );
+ RA.setUsing_atan2( v1, v0 );
RA.reduceToRange( dms::ZERO_TO_2PI );
- Dec.setUsing_asin( v[2] );
+ Dec.setUsing_asin( v2 );
+#undef v0
+#undef v1
+#undef v2
}
SkyPoint SkyPoint::deprecess( const KSNumbers *num, long double epoch ) {