Index: /trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- /trunk/psLib/src/astro/psEarthOrientation.c (revision 5482) +++ /trunk/psLib/src/astro/psEarthOrientation.c (revision 5483) @@ -9,6 +9,6 @@ * @author Robert Daniel DeSonia, MHPCC * -* @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-03 04:31:22 $ +* @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-11-07 20:52:43 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -52,4 +52,18 @@ static bool eocInitialized = false; +static psSphereRot *rotMatrix_To_Quat(double A[3][3]); + +static void earthPoleFree(psEarthPole *pole) +{ + // There are non dynamic allocated items +} + +psEarthPole *psEarthPoleAlloc(void) +{ + psEarthPole* pole = psAlloc(sizeof(psEarthPole)); + psMemSetDeallocator(pole, (psFreeFunc) earthPoleFree); + return pole; +} + static bool eocInit() { @@ -139,5 +153,6 @@ psFree(xTable); psFree(yTable); - psFree(sTable); + psFree(sTable); // There are non dynamic allocated items + xTable = NULL; @@ -179,6 +194,8 @@ apparent->dErr = 0.0; } - psSphere *rp = psSphereAlloc(); - psSphere *r_p = psSphereAlloc(); + // psSphere *rp = psSphereAlloc(); + psCube *rp = psCubeAlloc(); + // psSphere *r_p = psSphereAlloc(); + psCube *r_p = psCubeAlloc(); double mu = 0.0; double mu_p = 0.0; @@ -188,5 +205,9 @@ // mu = apparent->r * direction->r + apparent->d * direction->d; psCube* directionVector = psSphereToCube(direction); + if (directionVector == NULL) + printf("directionVector is null\n"); psCube* actualVector = psSphereToCube(actual); + if (directionVector == NULL) + printf("actualVector is null\n"); mu = acos(directionVector->x*actualVector->x + directionVector->y*actualVector->y + @@ -194,21 +215,41 @@ //rp = apparent - mu * direction; - rp->r = actual->r - mu * direction->r; - rp->d = actual->d - mu * direction->d; + // rp->r = actual->r - mu * direction->r; + // rp->d = actual->d - mu * direction->d; + rp->x = actualVector->x - mu * directionVector->x; + rp->y = actualVector->y - mu * directionVector->y; + rp->z = actualVector->z - mu * directionVector->z; mu_p = mu + speed * ((mu * mu - 1.0) / (1.0 - speed * mu)); //Not sure if this is right. ADD gets a scalar from division of rp (a vector) - psCube* rpVector = psSphereToCube(rp); - a = sqrt( (1.0 - mu_p * mu_p) / - (sqrt(rpVector->x*rpVector->x - + rpVector->y*rpVector->y - + rpVector->z*rpVector->z)) ); - + //Paul claims rp is the modulus or square of the vector components + /* psCube* rpVector = psSphereToCube(rp); + a = sqrt( (1.0 - mu_p * mu_p) / ( rpVector->x*rpVector->x + + rpVector->y*rpVector->y + rpVector->z*rpVector->z) ); + */ + a = (1.0 - mu_p * mu_p) / (rp->x * rp->x + rp->y * rp->y + rp->z * rp->z); + if (a < 0.0) { + printf("a is negative\n"); + psFree(rp); + psFree(r_p); + psFree(directionVector); + psFree(actualVector); + psFree(apparent) + return NULL; + } else + a = sqrt(a); //r_p = mu_p * direction + a * rp; - r_p->r = mu_p * direction->r + a * rp->r; - r_p->d = mu_p * direction->d + a * rp->d; - - *apparent = *r_p; + // r_p->r = mu_p * direction->r + a * rp->r; + // r_p->d = mu_p * direction->d + a * rp->d; + r_p->x = mu_p * directionVector->x + a * rp->x; + r_p->y = mu_p * directionVector->y + a * rp->y; + r_p->z = mu_p * directionVector->z + a * rp->z; + + psSphere *r_pSphere = psCubeToSphere(r_p); + if (r_pSphere == NULL) + printf("r_pSphere is null\n"); + + *apparent = *r_pSphere; /* psSphereRot *rot = NULL; @@ -223,7 +264,8 @@ psFree(rp); psFree(r_p); + psFree(r_pSphere) psFree(directionVector); psFree(actualVector); - psFree(rpVector); + // psFree(rpVector); return apparent; } @@ -280,5 +322,6 @@ deflection = SEC_TO_RAD(deflection); theta = atan(r0/PS_AU) * tan(deflection); - phi = sqrt( deflection*deflection - theta*theta ); + // phi = sqrt( deflection*deflection - theta*theta ); + phi = deflection * cos(asin(theta/deflection)); apparent->r = theta; apparent->d = phi; @@ -447,8 +490,74 @@ } +static psSphereRot *rotMatrix_To_Quat(double A[3][3]) +{ + int i; + psSphereRot *new = (psSphereRot*)psAlloc(sizeof(psSphereRot)); + // psMemSetDeallocator(new, (psFreeFunc)sphereRotFree); + //Convert rotation matrix to quaternions + double diag_sum[3]; + int maxi; + double recip; + diag_sum[0] = 1.0 + A[0][0] - A[1][1] - A[2][2]; + diag_sum[1] = 1.0 - A[0][0] + A[1][1] - A[2][2]; + diag_sum[2] = 1.0 - A[0][0] - A[1][1] + A[2][2]; + diag_sum[3] = 1.0 + A[0][0] + A[1][1] + A[2][2]; + + maxi = 0; + for (i = 1; i < 4; ++i) { + if (diag_sum[i] > diag_sum[maxi]) { + maxi = i; + } + } + + double p = 0.5 * sqrt(diag_sum[maxi]); + recip = 1.0 / (4.0 * p); + + if (maxi == 0) { + new->q0 = p; + new->q1 = recip * (A[0][1] + A[1][0]); + new->q2 = recip * (A[2][0] + A[0][2]); + new->q3 = recip * (A[1][2] - A[2][1]); + } else if (maxi == 1) { + new->q0 = recip * (A[0][1] + A[1][0]); + new->q1 = p; + new->q2 = recip * (A[1][2] + A[2][1]); + new->q3 = recip * (A[2][0] - A[0][2]); + } else if (maxi == 2) { + new->q0 = recip * (A[2][0] + A[0][2]); + new->q1 = recip * (A[1][2] + A[2][1]); + new->q2 = p; + new->q3 = recip * (A[0][1] - A[1][0]); + } else if (maxi == 3) { + new->q0 = recip * (A[1][2] - A[2][1]); + new->q1 = recip * (A[2][0] - A[0][2]); + new->q2 = recip * (A[0][1] - A[1][0]); + new->q3 = p; + } + return new; +} psSphereRot* psSphereRot_CEOtoGCRS(const psEarthPole *pole) { - return NULL; + PS_ASSERT_PTR_NON_NULL(pole,NULL); + double A[3][3]; + psSphereRot *out = NULL; + + //Setup the rotation matrix and scalar value, a, as outlined by the ADD + //XXX: Used formula for rotation matrix D from mathworld for z-axis rotation + double a = 1.0 / (1.0 + sqrt(1.0 - (pole->x*pole->x + pole->y*pole->y) ) ); + A[0][0] = (1.0 - a*pole->x*pole->x)*cos(pole->s) - a*pole->x*pole->y*sin(pole->s); + A[0][1] = -a*pole->x*pole->y*cos(pole->s) + (1.0 - a*pole->y*pole->y)*sin(pole->s); + A[0][2] = pole->x*cos(pole->s) + pole->y*sin(pole->s); + A[1][0] = -(1.0 - a*pole->x*pole->x)*sin(pole->s) - a*pole->x*pole->y*cos(pole->s); + A[1][1] = a*pole->x*pole->y*sin(pole->s) + (1.0 - a*pole->y*pole->y)*cos(pole->s); + A[1][2] = -pole->x*sin(pole->s) + pole->y*cos(pole->s); + A[2][0] = -pole->x; + A[2][1] = -pole->y; + A[2][2] = 1.0 - a*(pole->x*pole->x + pole->y*pole->y); + + out = rotMatrix_To_Quat(A); + + return out; } @@ -495,4 +604,61 @@ psSphereRot* psSphereRot_ITRStoTEO(const psEarthPole* motion) { - return NULL; -} + PS_ASSERT_PTR_NON_NULL(motion,NULL); + double A[3][3]; + psSphereRot *out = NULL; + + //Setup Rotation Matrix for transformation (x,y,z rotation) + //XXX: May need to be (z,y,x as in Mathworld?) + A[0][0] = cos(motion->x)*cos(-motion->s); + A[0][1] = sin(motion->y)*sin(motion->x)*cos(-motion->s) + cos(motion->y)*sin(-motion->s); + A[0][2] = -cos(motion->y)*sin(motion->x)*cos(-motion->s) + sin(motion->y)*sin(-motion->s); + A[1][0] = -cos(motion->x)*sin(-motion->s); + A[1][1] = -sin(motion->y)*sin(motion->x)*sin(-motion->s) + cos(motion->y)*cos(-motion->s); + A[1][2] = cos(motion->y)*sin(motion->x)*sin(-motion->s) + sin(motion->y)*cos(-motion->s); + A[2][0] = sin(motion->x); + A[2][1] = -sin(motion->y)*cos(motion->x); + A[2][2] = cos(motion->y)*cos(motion->x); + + //Convert rotation matrix to quaternions + double diag_sum[3]; + int maxi; + double recip; + diag_sum[0] = 1.0 + A[0][0] - A[1][1] - A[2][2]; + diag_sum[1] = 1.0 - A[0][0] + A[1][1] - A[2][2]; + diag_sum[2] = 1.0 - A[0][0] - A[1][1] + A[2][2]; + diag_sum[3] = 1.0 + A[0][0] + A[1][1] + A[2][2]; + + maxi = 0; + for (int i = 1; i < 4; ++i) { + if (diag_sum[i] > diag_sum[maxi]) { + maxi = i; + } + } + + double p = 0.5 * sqrt(diag_sum[maxi]); + recip = 1.0 / (4.0 * p); + + if (maxi == 0) { + out->q0 = p; + out->q1 = recip * (A[0][1] + A[1][0]); + out->q2 = recip * (A[2][0] + A[0][2]); + out->q3 = recip * (A[1][2] - A[2][1]); + } else if (maxi == 1) { + out->q0 = recip * (A[0][1] + A[1][0]); + out->q1 = p; + out->q2 = recip * (A[1][2] + A[2][1]); + out->q3 = recip * (A[2][0] - A[0][2]); + } else if (maxi == 2) { + out->q0 = recip * (A[2][0] + A[0][2]); + out->q1 = recip * (A[1][2] + A[2][1]); + out->q2 = p; + out->q3 = recip * (A[0][1] - A[1][0]); + } else if (maxi == 3) { + out->q0 = recip * (A[1][2] - A[2][1]); + out->q1 = recip * (A[2][0] - A[0][2]); + out->q2 = recip * (A[0][1] - A[1][0]); + out->q3 = p; + } + + return out; +} Index: /trunk/psLib/src/astro/psEarthOrientation.h =================================================================== --- /trunk/psLib/src/astro/psEarthOrientation.h (revision 5482) +++ /trunk/psLib/src/astro/psEarthOrientation.h (revision 5483) @@ -9,6 +9,6 @@ * @author Robert Daniel DeSonia, MHPCC * -* @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-03 04:31:22 $ +* @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-11-07 20:52:43 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -44,4 +44,7 @@ } psPrecessMethod; + +/** Allocates a new psEarthPole structure. */ +psEarthPole *psEarthPoleAlloc(void); /** Calculates the actual position of a star, given its apparent position and the Index: /trunk/psLib/src/astro/psSphereOps.h =================================================================== --- /trunk/psLib/src/astro/psSphereOps.h (revision 5482) +++ /trunk/psLib/src/astro/psSphereOps.h (revision 5483) @@ -7,6 +7,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-10-25 00:38:00 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-07 20:52:43 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -186,5 +186,5 @@ * coordinate systems. * - * @return psSphereTransform* transform for ICRS->Ecliptic coordinate systems + * @return psSphereRot* transform for ICRS->Ecliptic coordinate systems */ psSphereRot* psSphereRotICRSToEcliptic( @@ -195,5 +195,5 @@ * coordinate systems. * - * @return psSphereTransform* transform for Ecliptic->ICRS coordinate systems + * @return psSphereRot* transform for Ecliptic->ICRS coordinate systems */ psSphereRot* psSphereRotEclipticToICRS( @@ -204,4 +204,5 @@ * coordinate systems. * + * @return psSphereRot* new sphere rotation for ICRS to Galactic transformations. */ psSphereRot* psSphereRotICRSToGalactic(void); @@ -210,4 +211,5 @@ * coordinate systems. * + * @return psSphereRot* new sphere rotation for Galactic to ICRS transformations. */ psSphereRot* psSphereRotGalacticToICRS(void); Index: /trunk/psLib/test/astro/tst_psEarthOrientation.c =================================================================== --- /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5482) +++ /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5483) @@ -5,6 +5,6 @@ * @author d-Rob, MHPCC * -* @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-03 04:31:22 $ +* @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-11-07 20:52:43 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -63,8 +63,12 @@ psSphere *empty = NULL; - actual->r = 0.2; - actual->d = 0.2; - direction->r = 0.2035; - direction->d = 0.2035; + // actual->r = 0.2; + // actual->d = 0.2; + actual->r = DEG_TO_RAD(45.0); + actual->d = DEG_TO_RAD(30.0); + // direction->r = 0.2035; + // direction->d = 0.2035; + direction->r = DEG_TO_RAD(48.0); + direction->d = DEG_TO_RAD(20.7072); empty = psAberration(empty, apparent, direction, 0.1); @@ -82,5 +86,8 @@ apparent = psAberration(apparent, actual, direction, 0.4); - printf("\napparent = r,d = %.8g, %.8g\n", apparent->r, apparent->d); + if (apparent == NULL) + printf("\nApparent is NULL\n"); + else + printf("\napparent = r,d = %.8g, %.8g\n", apparent->r, apparent->d); psFree(apparent); Index: /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stderr =================================================================== --- /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stderr (revision 5482) +++ /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stderr (revision 5483) @@ -24,13 +24,4 @@ ---> TESTPOINT PASSED (psEarthOrientation{psGravityDeflect()} | tst_psEarthOrientation.c) - -/***************************** TESTPOINT ******************************************\ -* TestFile: tst_psEarthOrientation.c * -* TestPoint: psEarthOrientation{psEOCParallax()} * -* TestType: Positive * -\**********************************************************************************/ - - ----> TESTPOINT PASSED (psEarthOrientation{psEOCParallax()} | tst_psEarthOrientation.c) /***************************** TESTPOINT ******************************************\ @@ -63,9 +54,27 @@ /***************************** TESTPOINT ******************************************\ * TestFile: tst_psEarthOrientation.c * -* TestPoint: psEarthOrientation{psSphereRotTransforms()} * +* TestPoint: psEarthOrientation{psSphereRot_TEOtoCEO()} * * TestType: Positive * \**********************************************************************************/ ----> TESTPOINT PASSED (psEarthOrientation{psSphereRotTransforms()} | tst_psEarthOrientation.c) +---> TESTPOINT PASSED (psEarthOrientation{psSphereRot_TEOtoCEO()} | tst_psEarthOrientation.c) +/***************************** TESTPOINT ******************************************\ +* TestFile: tst_psEarthOrientation.c * +* TestPoint: psEarthOrientation{psSphereRot_CEOtoGCRS()} * +* TestType: Positive * +\**********************************************************************************/ + + +---> TESTPOINT PASSED (psEarthOrientation{psSphereRot_CEOtoGCRS()} | tst_psEarthOrientation.c) + +/***************************** TESTPOINT ******************************************\ +* TestFile: tst_psEarthOrientation.c * +* TestPoint: psEarthOrientation{psSphereRRot_ITRStoTEO()} * +* TestType: Positive * +\**********************************************************************************/ + + +---> TESTPOINT PASSED (psEarthOrientation{psSphereRRot_ITRStoTEO()} | tst_psEarthOrientation.c) + Index: /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stdout =================================================================== --- /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stdout (revision 5482) +++ /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stdout (revision 5483) @@ -1,4 +1,9 @@ -apparent = r,d = 0.10218436, 0.10218436 +apparent = r,d = 1.2163241, 1.0535409 Actual r,d = 0.2035,0.2035 Apparent r,d = 0.2035000000002, 0.2035000000391 + +The Value of T is = 1128530000.931 + +The Value of T is = 1126078455.931 +Output sphere rotation = 0.000000,0.000000,0.770554,0.637375