Changeset 5012

Timestamp:

Sep 12, 2005, 2:54:15 PM (21 years ago)

Author:

desonia

Message:

* empty log message *

Location:

trunk/psLib/src

Files:

• astro/psEarthOrientation.c (modified) (5 diffs)
• astro/psEarthOrientation.h (modified) (2 diffs)
• sys/psTrace.h (modified) (2 diffs)

trunk/psLib/src/astro/psEarthOrientation.c

@@ -8 +8 @@
 *  @author Robert Daniel DeSonia, MHPCC
 *
-*  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-07-12 19:27:27 $
+*  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-09-13 00:49:53 $
 *
 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
 */
 
+#include <math.h>
+#include <string.h>
+
+#include "psEarthOrientation.h"
+#include "psTime.h"
+#include "psArray.h"
+#include "psFunctions.h"
+#include "psVector.h"
+#include "psMetadata.h"
+#include "psMetadataConfig.h"
+#include "psError.h"
+#include "psErrorText.h"
+#include "psMemory.h"
+#include "psCoord.h"
+#include "psConstants.h"
+
+// Sun's Mass (src: Google's Calculator Service)
+#define PS_M   1.98892e30 /* kilograms */
+// Newton's Gravitational Constant (src:NIST)
+#define PS_G   6.6742e-11 /* m^3/kg/s^2 */
+// Speed of light in vacuum (src:NIST)
+#define PS_C0  299792458.0 /* m/s */
+// Average distance from earth to sun
+#define PS_AU 149597890000.0 /* meters */
+// Modified Julian Day 1/1/2000 00:00:00
+#define MJD_2000  51544.0
+// Days in Julian century
+#define JULIAN_CENTURY 36525.0
+
+static psArray* xTable = NULL;
+static psArray* yTable = NULL;
+static psArray* sTable = NULL;
+static psPolynomial1D* xPoly = NULL;
+static psPolynomial1D* yPoly = NULL;
+static psPolynomial1D* sPoly = NULL;
+static bool eocInitialized = false;
+
+static bool eocInit()
+{
+    int nFail = 0;
+
+    // Read config file
+    psMetadata* eocMetadata = psMetadataConfigParse(NULL,
+                              &nFail,
+                              p_psGetConfigFileName(),
+                              true);
+
+    if(eocMetadata == NULL) {
+        return false;
+    } else if(nFail != 0) {
+        return false;
+    }
+
+    bool success = false;
+    // Get table format
+    char* tableFormat = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.table.format");
+    if(! success || tableFormat == NULL) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
+                "psLib.eoc.precession.table.format");
+        return false;
+    }
+
+    char* xTableName = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.table.file.x");
+    if(! success || xTableName == NULL) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
+                "psLib.eoc.precession.x.file");
+        return false;
+    }
+
+    char* yTableName = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.table.file.y");
+    if(! success || yTableName == NULL) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
+                "psLib.eoc.precession.y.file");
+        return false;
+    }
+
+    char* sTableName = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.table.file.s");
+    if(! success || sTableName == NULL) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
+                "psLib.eoc.precession.s.file");
+        return false;
+    }
+    xTable = psVectorsReadFromFile(xTableName, tableFormat);
+    yTable = psVectorsReadFromFile(yTableName, tableFormat);
+    sTable = psVectorsReadFromFile(sTableName, tableFormat);
+    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;
+    }
+
+    psVector* xCoeff = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.poly.x");
+    if(! success || xCoeff == NULL || xCoeff->type.type != PS_TYPE_F64) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
+                "psLib.eoc.precession.poly.x");
+        return false;
+    }
+    psVector* yCoeff = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.poly.y");
+    if(! success || yCoeff == NULL || yCoeff->type.type != PS_TYPE_F64) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
+                "psLib.eoc.precession.poly.y");
+        return false;
+    }
+    psVector* sCoeff = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.poly.s");
+    if(! success || sCoeff == NULL || sCoeff->type.type != PS_TYPE_F64) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
+                "psLib.eoc.precession.poly.s");
+        return false;
+    }
+    xPoly = psPolynomial1DAlloc(xCoeff->n, PS_POLYNOMIAL_ORD);
+    memcpy(xPoly->coeff, xCoeff->data.F64,PSELEMTYPE_SIZEOF(PS_TYPE_F64)*xCoeff->n);
+    yPoly = psPolynomial1DAlloc(yCoeff->n, PS_POLYNOMIAL_ORD);
+    memcpy(yPoly->coeff, yCoeff->data.F64,PSELEMTYPE_SIZEOF(PS_TYPE_F64)*yCoeff->n);
+    sPoly = psPolynomial1DAlloc(sCoeff->n, PS_POLYNOMIAL_ORD);
+    memcpy(sPoly->coeff, sCoeff->data.F64,PSELEMTYPE_SIZEOF(PS_TYPE_F64)*sCoeff->n);
+
+    return true;
+}
+
+bool p_psEOCFinalize(void)
+{
+    psFree(xTable);
+    psFree(yTable);
+    psFree(sTable);
+
+    xTable = NULL;
+    yTable = NULL;
+    sTable = NULL;
+
+    psFree(xPoly);
+    psFree(yPoly);
+    psFree(sPoly);
+
+    xPoly = NULL;
+    yPoly = NULL;
+    sPoly = NULL;
+
+    eocInitialized = false;
+
+    return true;
+}
+
 psSphere *psAberration(psSphere *apparent, const psSphere *actual, const psSphere *direction, double speed)
 {
+
+
     return NULL;
 }
@@ -21 +166 @@
 psSphere *psGravityDeflection(psSphere *apparent, psSphere *actual, psSphere *sun)
 {
+    // calculating the apparent angle from the actual angle and the sun position
+
+    // first, calculate the angle between the sun vector and the actual vector
+
+    // Moving to cartesian first:  XXX -- is this required?
+    psCube* sunVector = psSphereToCube(sun);
+    psCube* actualVector = psSphereToCube(actual);
+
+    // use dot product to calculate the angle of separation
+    // N.B., assuming the psSphereToCube function returns a unit vector.
+    double theta = acos(sunVector->x*actualVector->x +
+                        sunVector->y*actualVector->y +
+                        sunVector->z*actualVector->z);
+
+    double r0 = PS_AU * tan(theta);
+
+    double deflection = 4.0*PS_G*PS_M/(PS_C0*PS_C0*r0);
+
+    // make sure the deflection is not greater than 1.75 arcsec
+    double limit = SEC_TO_RAD(1.75);
+    if (deflection > limit) {
+        deflection = limit;
+    }
+
+    // bend the actual vector away from the sun vector by deflection angle.
+
+
     return NULL;
 }
@@ -27 +199 @@
 double psEOC_ParallaxFactor(const psSphere *coords, const psTime *time)
 {
+
+
     return NAN;
 }
@@ -32 +206 @@
 psEarthPole *psEOC_PrecessionModel(const psTime *time)
 {
-    return NULL;
+    // Check for null parameter
+    PS_ASSERT_PTR_NON_NULL(time, NULL);
+
+    // Convert psTime to MJD
+    double MJD = psTimeToMJD(time);
+
+    // Calculate number of Julian centuries since 2000
+    double t = ( MJD - MJD_2000 ) / JULIAN_CENTURY;
+    double t2 = t*t;
+    double t3 = t*t*t;
+    double t4 = t*t*t*t;
+
+    // N.B., the following formulae are from the ADD
+    // Mean Anomaly of the Moon
+    double F[14];
+    F[0] = DEG_TO_RAD(134.96340251) +
+           SEC_TO_RAD(1717915923.2178)*t +
+           SEC_TO_RAD(31.8792)*t2 +
+           SEC_TO_RAD(0.051635)*t3 -
+           SEC_TO_RAD(0.00024470)*t4;
+
+    // Mean Anomaly of the Sun
+    F[1] = DEG_TO_RAD(357.52910918) +
+           SEC_TO_RAD(129596581.0481)*t -
+           SEC_TO_RAD(0.5532)*t2 +
+           SEC_TO_RAD(0.000136)*t3 -
+           SEC_TO_RAD(0.00001149)*t4;
+
+    // L − Omega
+    F[2] = DEG_TO_RAD(93.27209062) +
+           SEC_TO_RAD(1739527262.8478)*t -
+           SEC_TO_RAD(12.7512)*t2 -
+           SEC_TO_RAD(0.001037)*t3 +
+           SEC_TO_RAD(0.00000417)*t4;
+
+    // Mean Elongation of the Moon from the Sun
+    F[3] = DEG_TO_RAD(297.85019547) +
+           SEC_TO_RAD(1602961601.2090)*t -
+           SEC_TO_RAD(6.3706)*t2 +
+           SEC_TO_RAD(0.006593)*t3 -
+           SEC_TO_RAD(0.00003169)*t4;
+
+    // Mean Longitude of the Ascending Node of the Moon
+    F[4] = DEG_TO_RAD(125.04455501) -
+           SEC_TO_RAD(6962890.5431)*t +
+           SEC_TO_RAD(7.4722)*t2 +
+           SEC_TO_RAD(0.007702)*t3 -
+           SEC_TO_RAD(0.0000593)*t4;
+
+    F[5] = 4.402608842 + 2608.7903141574*t;
+    F[6] = 3.176146697 + 1021.3285546211*t;
+    F[7] = 1.753470314 + 628.3075849991*t;
+    F[8] = 6.203480913 + 334.0612426700*t;
+    F[9] = 0.599546497 + 52.9690962641*t;
+    F[10] = 0.874016757 + 21.3299104960*t;
+    F[11] = 5.481293872 + 7.4781598567*t;
+    F[12] = 5.311886287 + 3.8133035638*t;
+    F[13] = 0.024381750*t + 0.00000538691*t2;
+
+    // Check if EOC data loaded
+    if(! eocInitialized) {
+        eocInitialized = eocInit();
+        if(!eocInitialized) {
+            // XXX: Move error message.
+            psError(PS_ERR_UNKNOWN, false,
+                    "Could not initialize EOC tables -- check data files.");
+            return NULL;
+        }
+    }
+
+    // calculate the polynomial portion first
+    double X = psPolynomial1DEval(xPoly,t);
+    double Y = psPolynomial1DEval(yPoly,t);
+    double S = psPolynomial1DEval(sPoly,t);
+
+    // now calculate the non-poly portion from the tables
+
+    psF64* cols[17];
+    for (int lcv = 0; lcv < 17; lcv++) {
+        cols[lcv] = ((psVector*)(xTable->data[lcv]))->data.F64;
+    }
+    int numRows = ((psVector*)(xTable->data[0]))->n;
+    for (int lcv = 0; lcv < numRows; lcv++) {
+        // arg = sum w_(i,j,k)*F_k
+        double arg = 0.0;
+        for (int k = 0; k < 14; k++) {
+            arg += cols[k+3][lcv]*F[k];
+        }
+        double tj = pow(t,cols[0][lcv]);
+        double as = cols[1][lcv];
+        double ac = cols[2][lcv];
+
+        X += as*tj*sin(arg) + ac*tj*cos(arg);
+    }
+
+    for (int lcv = 0; lcv < 17; lcv++) {
+        cols[lcv] = ((psVector*)(yTable->data[lcv]))->data.F64;
+    }
+    numRows = ((psVector*)(yTable->data[0]))->n;
+    for (int lcv = 0; lcv < numRows; lcv++) {
+        // arg = sum w_(i,j,k)*F_k
+        double arg = 0.0;
+        for (int k = 0; k < 14; k++) {
+            arg += cols[k+3][lcv]*F[k];
+        }
+        double tj = pow(t,cols[0][lcv]);
+        double as = cols[1][lcv];
+        double ac = cols[2][lcv];
+
+        Y += as*tj*sin(arg) + ac*tj*cos(arg);
+    }
+
+    for (int lcv = 0; lcv < 17; lcv++) {
+        cols[lcv] = ((psVector*)(sTable->data[lcv]))->data.F64;
+    }
+    numRows = ((psVector*)(sTable->data[0]))->n;
+    for (int lcv = 0; lcv < numRows; lcv++) {
+        // arg = sum w_(i,j,k)*F_k
+        double arg = 0.0;
+        for (int k = 0; k < 14; k++) {
+            arg += cols[k+3][lcv]*F[k];
+        }
+        double tj = pow(t,cols[0][lcv]);
+        double as = cols[1][lcv];
+        double ac = cols[2][lcv];
+
+        S += as*tj*sin(arg) + ac*tj*cos(arg);
+    }
+
+    // now, the tables for S actually gives S + XY/2, so let's get the real S now
+    S -= X*Y/2.0;
+
+    psEarthPole* pole = psAlloc(sizeof(psEarthPole));
+    pole->x = X;
+    pole->y = Y;
+    pole->s = S;
+
+    return pole;
 }
 
@@ -42 +353 @@
 
 
-psSphereRot *psSphereRot_CEOtoGCRS(const psEarthPole *pole)
-{
-    return NULL;
-}
-
-
-psSphereRot *psSphereRot_TEOtoCEO(const psTime *time)
-{
-    return NULL;
-}
-
-
-psEarthPole *psEOC_GetPolarMotion(const psTime *time, psTimeBulletin bulletin)
-{
-    return NULL;
-}
-
-
-psEarthPole *psEOC_PolarTideCorr(const psTime *time)
-{
-    return NULL;
-}
-
-
-psEarthPole *psEOC_NutationCorr(psTime *time)
-{
-    return NULL;
-}
-
-
-psSphereRot *psSphereRot_ITRStoTEO(const psEarthPole *motion)
-{
-    return NULL;
-}
-
-
+psSphereRot* psSphereRot_CEOtoGCRS(const psEarthPole *pole)
+{
+    return NULL;
+}
+
+
+psSphereRot* psSphereRot_TEOtoCEO(const psTime *time)
+{
+    return NULL;
+}
+
+
+psEarthPole* psEOC_GetPolarMotion(const psTime *time, psTimeBulletin bulletin)
+{
+    return NULL;
+}
+
+
+psEarthPole* psEOC_PolarTideCorr(const psTime *time)
+{
+    return NULL;
+}
+
+
+psEarthPole* psEOC_NutationCorr(psTime *time)
+{
+    return NULL;
+}
+
+
+psSphereRot* psSphereRot_ITRStoTEO(const psEarthPole* motion)
+{
+    return NULL;
+}

trunk/psLib/src/astro/psEarthOrientation.h

@@ -8 +8 @@
 *  @author Robert Daniel DeSonia, MHPCC
 *
-*  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-07-12 19:27:27 $
+*  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-09-13 00:50:04 $
 *
 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -19 +19 @@
 #include "psCoord.h"
 #include "psTime.h"
+#include "psSphereOps.h"
 
 typedef struct

trunk/psLib/src/sys/psTrace.h

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-08 00:27:50 $
+ *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-13 00:54:15 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -16 +16 @@
 #if !defined(PS_TRACE_H)
 #define PS_TRACE_H 1
+#include <stdarg.h>
 
 #define PS_UNKNOWN_TRACE_LEVEL -9999   // we don't know this name's level