summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAkarsh Simha <akarsh@kde.org>2016-09-28 07:30:25 (GMT)
committerAkarsh Simha <akarsh@kde.org>2016-09-29 02:20:27 (GMT)
commitd9158d4d986b29da2950a402612f01a919480fc6 (patch)
tree5b33145536f1ef993327df44430fd4aed45239e5
parent32de27498abdaf3d001fee75c5a0b11a2e056178 (diff)
Experiment: Is using Eigen::Matrix3d faster than 2D array and for loop?
Results: It is not, but it produces cleaner code.
-rw-r--r--kstars/ksnumbers.cpp79
-rw-r--r--kstars/ksnumbers.h28
-rw-r--r--kstars/skyobjects/skypoint.cpp24
3 files changed, 69 insertions, 62 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 47cb5cd..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,17 +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 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]; }
+ /** @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
@@ -130,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 fea1ba3..4f810f2 100644
--- a/kstars/skyobjects/skypoint.cpp
+++ b/kstars/skyobjects/skypoint.cpp
@@ -190,7 +190,8 @@ void SkyPoint::setFromEcliptic( const CachingDms *Obliquity, const dms& EcLong,
void SkyPoint::precess( const KSNumbers *num ) {
double cosRA0, sinRA0, cosDec0, sinDec0;
- double v[3], s[3];
+ const Eigen::Matrix3d &precessionMatrix = num->p2();
+ Eigen::Vector3d v, s;
RA0.SinCos( sinRA0, cosRA0 );
Dec0.SinCos( sinDec0, cosDec0 );
@@ -201,24 +202,19 @@ void SkyPoint::precess( const KSNumbers *num ) {
// FIXME: 1. We should be using eigen / better algorithms for
// matrix multiplication
- // 2. We should NOT be using matrix multiplication, but use
- // quaternion multiplication instead!
- // 3. But then, since we are using spherical coordinates
- // (ra,dec) as our primary representation inside of
- // KStars, we should instead be using spherical trig for
- // best performance!
- //
+
// According to callgrind, the call KSNumbers::p2( int, int ),
// which is repeated 9 times per precess, has a similar cycle cost
// as atan2(), so this could be important to fix!
//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];
- }
- }
+ // 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( v[1], v[0] );