Commit 66a7ea2f authored by Akarsh Simha's avatar Akarsh Simha
Browse files

Experiment: Explicitly write out the computations of precession

parent 743cba78
......@@ -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
//
......
......@@ -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
......
......@@ -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 ) {
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment