Changeset 6328

Timestamp:

Feb 6, 2006, 11:57:04 AM (20 years ago)

Author:

drobbin

Message:

added comments

File:

• trunk/psLib/src/astro/psEarthOrientation.c (modified) (14 diffs)

trunk/psLib/src/astro/psEarthOrientation.c

@@ -8 +8 @@
  *  @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 +18 @@
 
 #include "psEarthOrientation.h"
-//#include "psTime.h"
 #include "psArray.h"
 #include "psPolynomial.h"
@@ -43 +42 @@
 #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 +91 @@
                               p_psGetConfigFileName(),
                               true);
-
+    //Make sure reading of config file worked correctly
     if(eocMetadata == NULL) {
         return false;
@@ -87 +99 @@
 
     bool success = false;
-    // Get table format
+    // Get table formats & error if lookups fail.
     char* tableFormat = psMetadataLookupStr(&success, eocMetadata,
                                             "psLib.eoc.precession.table.format");
@@ -105 +117 @@
     }
 
-    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 +192 @@
     }
 
+    //Extract the necessary data to setup the tables
     xTable = psVectorsReadFromFile(xTableName, tableFormat);
     yTable = psVectorsReadFromFile(yTableName, tableFormat);
@@ -185 +199 @@
     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 +211 @@
     }
 
+    //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 +233 @@
         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 +278 @@
                        double speed)
 {
+    //Check for valid inputs
     PS_ASSERT_PTR_NON_NULL(actual, NULL);
     PS_ASSERT_PTR_NON_NULL(direction, NULL);
@@ -272 +287 @@
     }
 
-    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 +293 @@
     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 +326 @@
         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;
 }