Index: trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- trunk/psLib/src/astro/psEarthOrientation.c (revision 6326) +++ trunk/psLib/src/astro/psEarthOrientation.c (revision 6328) @@ -8,6 +8,6 @@ * @author Robert Daniel DeSonia, MHPCC * - * @version $Revision: 1.33 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-03 00:11:54 $ + * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-02-06 21:57:04 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -18,5 +18,4 @@ #include "psEarthOrientation.h" -//#include "psTime.h" #include "psArray.h" #include "psPolynomial.h" @@ -43,14 +42,27 @@ #define JULIAN_CENTURY 36525.0 +//The arrays used for storage of x,y,&s-values extracted from the precession data tables. +// Created in p_psEOCInit from finalsTable and used in psPrecessionModel. static psArray* xTable = NULL; static psArray* yTable = NULL; static psArray* sTable = NULL; + +//used for storage of the IERS precession table. Should contain MJD date followed by the +// corresponding dX & dY values from the finals2000A table. (Used in psPrecessionCorr.) static psArray* iersTable = NULL; +//used for precession data. Should contain MJD date followed by the corresponding X, Y, & +// UT1-UTC values from the finals2000A table. (Used in psPrecessionModel, psGetPolarMotion). static psArray* finalsTable = NULL; + +//1D Polynomials created from the coefficients given by the tab5.2?.txt files. Used +// by psPrecessionModel to determine psEarthPole components. static psPolynomial1D* xPoly = NULL; static psPolynomial1D* yPoly = NULL; static psPolynomial1D* sPoly = NULL; + +//Boolean variable to tell whether the above EOC data has been initialized. static bool eocInitialized = false; +//Internal function used for conversion from a 3x3 rotation matrix to a quaternion (psSphereRot). static psSphereRot *rotMatrix_To_Quat(double A[3][3]); @@ -79,5 +91,5 @@ p_psGetConfigFileName(), true); - + //Make sure reading of config file worked correctly if(eocMetadata == NULL) { return false; @@ -87,5 +99,5 @@ bool success = false; - // Get table format + // Get table formats & error if lookups fail. char* tableFormat = psMetadataLookupStr(&success, eocMetadata, "psLib.eoc.precession.table.format"); @@ -105,5 +117,6 @@ } - char* yTableName = psMetadataLookupStr(&success, eocMetadata, "psLib.eoc.precession.table.file.y"); + char* yTableName = psMetadataLookupStr(&success, eocMetadata, + "psLib.eoc.precession.table.file.y"); if(! success || yTableName == NULL) { psFree(tableFormat); @@ -179,4 +192,5 @@ } + //Extract the necessary data to setup the tables xTable = psVectorsReadFromFile(xTableName, tableFormat); yTable = psVectorsReadFromFile(yTableName, tableFormat); @@ -185,13 +199,11 @@ finalsTable = psVectorsReadFromFile(finalsTableName, finalsTableFormat); + //Check that the data extraction was performed successfully. if(xTable == NULL || yTable == NULL || sTable == NULL) { - // XXX: need to move the error message to the error messages file. psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Failed to read the precession-nutation model tables."); return false; } - if(iersTable == NULL || finalsTable == NULL) { - // XXX: need to move the error message to the error messages file. psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Failed to read the IERS/finals tables."); @@ -199,4 +211,5 @@ } + //Load the coefficient data to setup the polynomials. psVector* xCoeff = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.poly.x"); if(! success || xCoeff == NULL || xCoeff->type.type != PS_TYPE_F64) { @@ -220,4 +233,5 @@ return false; } + //Allocate the X,Y,&S polynomials to be used for earthpole calculations xPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, xCoeff->n); memcpy(xPoly->coeff, xCoeff->data.F64,PSELEMTYPE_SIZEOF(PS_TYPE_F64)*xCoeff->n); @@ -264,4 +278,5 @@ double speed) { + //Check for valid inputs PS_ASSERT_PTR_NON_NULL(actual, NULL); PS_ASSERT_PTR_NON_NULL(direction, NULL); @@ -272,15 +287,5 @@ } - if (apparent == NULL) { - apparent = psSphereAlloc(); - } else { - apparent->r = 0.0; - apparent->d = 0.0; - apparent->rErr = 0.0; - apparent->dErr = 0.0; - } - // psSphere *rp = psSphereAlloc(); psCube *rp = psCubeAlloc(); - // psSphere *r_p = psSphereAlloc(); psCube *r_p = psCubeAlloc(); double mu = 0.0; @@ -288,29 +293,31 @@ double a = 0.0; - psCube* directionVector = psSphereToCube(direction); - if (directionVector == NULL) - printf("directionVector is null\n"); - psCube* actualVector = psSphereToCube(actual); - if (actualVector == NULL) - printf("actualVector is null\n"); - // mu = acos(directionVector->x*actualVector->x + + //Convert the spherical input coords direction and actual to cubic coords for computations. + psCube* directionVector = NULL; + directionVector = psSphereToCube(direction); + if (directionVector == NULL) { + PS_ASSERT_PTR_NON_NULL(directionVector, NULL); + } + psCube* actualVector = NULL; + actualVector = psSphereToCube(actual); + if (actualVector == NULL) { + PS_ASSERT_PTR_NON_NULL(actualVector, NULL); + } + + //Calculate mu as the dot-product of direction and actual (from ADD) mu = (directionVector->x*actualVector->x + directionVector->y*actualVector->y + directionVector->z*actualVector->z); + //Calculate r-perpendicular (as in sec 3.5.3.2 of ADD). actual = r-hat, direction = beta-hat. rp->x = actualVector->x - mu * directionVector->x; rp->y = actualVector->y - mu * directionVector->y; rp->z = actualVector->z - mu * directionVector->z; - + //Calculate mu-prime. (Eqn. 129 of ADD). 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) - //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) ); - */ + //Calculate a-value. (a = sqrt[ (1 - mu-prime^2) / modulus(r-perpendicular) ] ) 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"); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Invalid parameter. a-value cannot be negative.\n"); psFree(rp); psFree(r_p); @@ -319,28 +326,27 @@ psFree(apparent) return NULL; - } else + } else { a = sqrt(a); - + } + //Calculate r-prime values. r-prime = mu-prime * beta-hat + a * r-perpendicular 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; - - //XXX: Must be a sign error somewhere above? Magnitude of change is correct but wrong way... - //This will fix the problem but is somewhat of a hack. - // r_p->x = actualVector->x - (r_p->x - actualVector->x); - // r_p->y = actualVector->y - (r_p->y - actualVector->y); - // r_p->z = actualVector->z - (r_p->z - actualVector->z); - - psSphere *r_pSphere = psCubeToSphere(r_p); - if (r_pSphere == NULL) - printf("r_pSphere is null\n"); - - *apparent = *r_pSphere; + //If apparent is non-NULL, deallocate it for use with psCubeToSphere + if (apparent != NULL) { + psFree(apparent); + } + //Convert r-prime to spherical coordinates for output. + apparent = psCubeToSphere(r_p); + if (apparent == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + "psCubeToSphere returned a NULL sphere in psAberration.\n"); + return NULL; + } psFree(rp); psFree(r_p); - psFree(r_pSphere) psFree(directionVector); psFree(actualVector); - // psFree(rpVector); + return apparent; }