summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAkarsh Simha <akarsh@kde.org>2016-09-28 11:55:22 (GMT)
committerAkarsh Simha <akarsh@kde.org>2016-09-29 02:20:27 (GMT)
commit4fd4a0067604c52b7a52da90fbbf2646daf53304 (patch)
treeb1614736bcd3009c461131a8f3134d2199a08fc1
parent66a7ea2f3829bff896128836479bdf5fb365df48 (diff)
Revert "Experiment: Explicitly write out the computations of precession"
This reverts commit 1babe81f4b721f82663f84fbd5c15703b2ec661d. The improvement afforded by explicitly writing out the computation, given the degree of compiler optimizations, is negligible. In my profiling, it turned out that the commit improved the performance by 8 CPU cycles on average per function call at the expense of substantially reduced readability! Not worth it.
-rw-r--r--kstars/ksnumbers.cpp79
-rw-r--r--kstars/ksnumbers.h30
-rw-r--r--kstars/skyobjects/skypoint.cpp27
3 files changed, 70 insertions, 66 deletions
diff --git a/kstars/ksnumbers.cpp b/kstars/ksnumbers.cpp
index b1e0b44..5193472 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,29 +313,30 @@ void KSNumbers::updateValues( long double jd ) {
ZP.SinCos( SZ, CZ );
//P1 is used to precess from any epoch to J2000
- // 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;
+ // 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;
//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 2ce8a62..aba30f4 100644
--- a/kstars/ksnumbers.h
+++ b/kstars/ksnumbers.h
@@ -24,6 +24,10 @@
#include "cachingdms.h"
+#include <Eigen/Dense>
+
+using Eigen::Matrix3d;
+
/** @class KSNumbers
*
*There are several time-dependent values used in position calculations,
@@ -95,21 +99,23 @@ 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 P2 precession array at position [i1][i2] */
- inline double p2( int i1, int i2 ) const { return P2[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 P1B precession array at position [i1][i2] */
- inline double p1b( int i1, int i2 ) const { return P1B[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 P2B precession array at position [i1][i2] */
- inline double p2b( int i1, int i2 ) const { return P2B[i1][i2]; }
+ /** @return element of P1B precession array at position (i1, i2) */
+ inline double p1b( int i1, int i2 ) const { return P1B(i1, i2); }
- inline const double *p1() const { return (*P1); }
+ /** @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 *p2() const { return (*P2); }
+ /** @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; }
/**
*@short compute constant values that need to be computed only once per instance of the application
@@ -134,7 +140,7 @@ private:
dms XP, YP, ZP, XB, YB, ZB;
double CX, SX, CY, SY, CZ, SZ;
double CXB, SXB, CYB, SYB, CZB, SZB;
- double P1[3][3], P2[3][3], P1B[3][3], P2B[3][3];
+ Matrix3d P1, P2, P1B, P2B;
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 c7f86af..6f9ae9a 100644
--- a/kstars/skyobjects/skypoint.cpp
+++ b/kstars/skyobjects/skypoint.cpp
@@ -190,13 +190,8 @@ void SkyPoint::setFromEcliptic( const CachingDms *Obliquity, const dms& EcLong,
void SkyPoint::precess( const KSNumbers *num ) {
double cosRA0, sinRA0, cosDec0, sinDec0;
- 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
+ const Eigen::Matrix3d &precessionMatrix = num->p2();
+ Eigen::Vector3d v, s;
RA0.SinCos( sinRA0, cosRA0 );
Dec0.SinCos( sinDec0, cosDec0 );
@@ -213,16 +208,18 @@ 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.
-#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]
+ // 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;
+
//Extract RA, Dec from the vector:
- RA.setUsing_atan2( v1, v0 );
+ RA.setUsing_atan2( v[1], v[0] );
RA.reduceToRange( dms::ZERO_TO_2PI );
- Dec.setUsing_asin( v2 );
-#undef v0
-#undef v1
-#undef v2
+ Dec.setUsing_asin( v[2] );
}
SkyPoint SkyPoint::deprecess( const KSNumbers *num, long double epoch ) {