IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5533


Ignore:
Timestamp:
Nov 16, 2005, 5:59:05 PM (21 years ago)
Author:
drobbin
Message:

Implemented function NutationCorr and test. Has yet to be verified.

Location:
trunk/psLib
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/astro/psEarthOrientation.c

    r5531 r5533  
    88 *  @author Robert Daniel DeSonia, MHPCC
    99 *
    10  *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2005-11-17 00:07:25 $
     10 *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2005-11-17 03:59:05 $
    1212 *
    1313 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    493493    S -= X*Y/2.0;
    494494
    495     psEarthPole* pole = psAlloc(sizeof(psEarthPole));
     495    //    psEarthPole* pole = psAlloc(sizeof(psEarthPole));
     496    psEarthPole* pole = psEarthPoleAlloc();
    496497    pole->x = X;
    497498    pole->y = Y;
    … …  
    850851psEarthPole* psEOC_NutationCorr(psTime *time)
    851852{
    852     return NULL;
     853    // Check for null parameter
     854    PS_ASSERT_PTR_NON_NULL(time, NULL);
     855    if (time->type == PS_TIME_UT1) {
     856        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     857                "Invalid time input.  Time cannot be of type UT1.\n");
     858        return NULL;
     859    }
     860
     861    // Convert psTime to MJD
     862    double MJD = psTimeToMJD(time);
     863
     864    // Calculate number of Julian centuries since 2000
     865    double t = ( MJD - MJD_2000 ) / JULIAN_CENTURY;
     866    double t2 = t*t;
     867    double t3 = t*t*t;
     868    double t4 = t*t*t*t;
     869
     870    //XXX: I think the t's should be inside of the SEC_TO_RAD conversion.
     871    //Check this and for Precession Model as well!
     872    double F[5];
     873    F[0] = DEG_TO_RAD(134.96340251) +
     874           SEC_TO_RAD(1717915923.2178)*t +
     875           SEC_TO_RAD(31.8792)*t2 +
     876           SEC_TO_RAD(0.051635)*t3 -
     877           SEC_TO_RAD(0.00024470)*t4;
     878
     879    // Mean Anomaly of the Sun
     880    F[1] = DEG_TO_RAD(357.52910918) +
     881           SEC_TO_RAD(129596581.0481)*t -
     882           SEC_TO_RAD(0.5532)*t2 +
     883           SEC_TO_RAD(0.000136)*t3 -
     884           SEC_TO_RAD(0.00001149)*t4;
     885
     886    // L − Omega
     887    F[2] = DEG_TO_RAD(93.27209062) +
     888           SEC_TO_RAD(1739527262.8478)*t -
     889           SEC_TO_RAD(12.7512)*t2 -
     890           SEC_TO_RAD(0.001037)*t3 +
     891           SEC_TO_RAD(0.00000417)*t4;
     892
     893    // Mean Elongation of the Moon from the Sun
     894    F[3] = DEG_TO_RAD(297.85019547) +
     895           SEC_TO_RAD(1602961601.2090)*t -
     896           SEC_TO_RAD(6.3706)*t2 +
     897           SEC_TO_RAD(0.006593)*t3 -
     898           SEC_TO_RAD(0.00003169)*t4;
     899
     900    // Mean Longitude of the Ascending Node of the Moon
     901    F[4] = DEG_TO_RAD(125.04455501) -
     902           SEC_TO_RAD(6962890.5431)*t +
     903           SEC_TO_RAD(7.4722)*t2 +
     904           SEC_TO_RAD(0.007702)*t3 -
     905           SEC_TO_RAD(0.0000593)*t4;
     906
     907    //argument values taken from table 5.1 in IERS techical note No.32
     908    //http://maia.usno.navy.mil/conv2000/chapter5/tn32_c5.pdf, p38
     909    //Units are in micro-arcseconds here and must be converted to radians before using
     910    double w_l[10] = {SEC_TO_RAD(-1.0), SEC_TO_RAD(-1.0), SEC_TO_RAD(1.0), 0.0, 0.0,
     911                      SEC_TO_RAD(-1.0), 0.0, 0.0, 0.0, SEC_TO_RAD(1.0)};
     912    double w_l_p[10] = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
     913    double w_F[10] = {SEC_TO_RAD(-2.0), SEC_TO_RAD(-2.0), SEC_TO_RAD(-2.0),
     914                      SEC_TO_RAD(-2.0), SEC_TO_RAD(-2.0), 0.0, SEC_TO_RAD(-2.0), 0.0,0.0,0.0};
     915    double w_D[10] = {0.0,0.0, SEC_TO_RAD(-2.0), 0.0,0.0,0.0, SEC_TO_RAD(2.0), 0.0,0.0,0.0};
     916    double w_Omega[10] = {SEC_TO_RAD(-1.0), SEC_TO_RAD(-2.0), SEC_TO_RAD(-2.0),
     917                          SEC_TO_RAD(-1.0), SEC_TO_RAD(-2.0), 0.0, SEC_TO_RAD(-2.0), 0.0, SEC_TO_RAD(-1.0), 0.0};
     918    double xp_sin[10] = {SEC_TO_RAD(-0.44), SEC_TO_RAD(-2.31), SEC_TO_RAD(-0.44),
     919                         SEC_TO_RAD(-2.14), SEC_TO_RAD(-11.36), SEC_TO_RAD(0.84), SEC_TO_RAD(-4.76),
     920                         SEC_TO_RAD(14.27), SEC_TO_RAD(1.93), SEC_TO_RAD(0.76)};
     921    double xp_cos[10] = {SEC_TO_RAD(0.25), SEC_TO_RAD(1.32), SEC_TO_RAD(0.25), SEC_TO_RAD(1.23),
     922                         SEC_TO_RAD(6.52), SEC_TO_RAD(-0.48), SEC_TO_RAD(2.73), SEC_TO_RAD(-8.19),
     923                         SEC_TO_RAD(-1.11), SEC_TO_RAD(-0.43)};
     924    double yp_sin[10] = {SEC_TO_RAD(-0.25), SEC_TO_RAD(-1.32), SEC_TO_RAD(-0.25),
     925                         SEC_TO_RAD(-1.23), SEC_TO_RAD(-6.52), SEC_TO_RAD(0.48), SEC_TO_RAD(-2.73),
     926                         SEC_TO_RAD(8.19), SEC_TO_RAD(1.11), SEC_TO_RAD(0.43)};
     927    double yp_cos[10] = {SEC_TO_RAD(-0.44), SEC_TO_RAD(-2.31), SEC_TO_RAD(-0.44),
     928                         SEC_TO_RAD(-2.14), SEC_TO_RAD(-11.36), SEC_TO_RAD(0.84), SEC_TO_RAD(-4.76),
     929                         SEC_TO_RAD(14.27), SEC_TO_RAD(1.93), SEC_TO_RAD(0.76)};
     930
     931    double X = 0.0;
     932    double Y = 0.0;
     933    double arg = 0.0;
     934    //This is from eqn 131 in the ADD - Note: pj_tj isn't included the first time.
     935    //XXX: The xp_sin, yp_cos, etc. may need to be multiplied by pow(t,i) here? adding now...
     936    double tj = 0.0;
     937
     938    // calculate the polynomial portion first - the pj * t^j (poly coeff's)
     939    // Check if EOC data loaded
     940    if(! eocInitialized) {
     941        eocInitialized = p_psEOCInit();
     942        if(!eocInitialized) {
     943            // XXX: Move error message.
     944            psError(PS_ERR_UNKNOWN, false,
     945                    "Could not initialize EOC tables -- check data files.");
     946            return NULL;
     947        }
     948    }
     949    X = psPolynomial1DEval(xPoly,t);
     950    Y = psPolynomial1DEval(yPoly,t);
     951    for (int i = 0; i < 10; i++) {
     952        //        tj = SEC_TO_RAD(pow(t, i));
     953        tj = pow(t, i);
     954        arg = w_l[i]*F[0] + w_l_p[i]*F[1] + w_F[i]*F[2] + w_D[i]*F[3] + w_Omega[i]*F[4];
     955        X += xp_sin[i] * tj * sin(arg) + xp_cos[i] * tj * cos(arg);
     956        Y += yp_sin[i] * tj * sin(arg) + yp_cos[i] * tj * cos(arg);
     957    }
     958
     959    psEarthPole *pole = psEarthPoleAlloc();
     960    pole->x = X;
     961    pole->y = Y;
     962    pole->s = SEC_TO_RAD(4.7e-5) * t;   //XXX:  This conv. should include t?
     963
     964    return pole;
    853965}
    854966

  • trunk/psLib/src/astro/psEarthOrientation.h

    r5522 r5533  
    99*  @author Robert Daniel DeSonia, MHPCC
    1010*
    11 *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
    12 *  @date $Date: 2005-11-16 02:40:13 $
     11*  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
     12*  @date $Date: 2005-11-17 03:59:05 $
    1313*
    1414*  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    4545psPrecessMethod;
    4646
     47/** Initialize EOC.
     48 *
     49 *  Parses IERS table data and creates the polynomials used by certain EOC functions.
     50 *
     51 *  @return bool:       True if successful, otherwise false.
     52 */
    4753bool p_psEOCInit(void);
    4854

  • trunk/psLib/test/astro/tst_psEarthOrientation.c

    r5524 r5533  
    55*  @author d-Rob, MHPCC
    66*
    7 *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
    8 *  @date $Date: 2005-11-16 20:52:23 $
     7*  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
     8*  @date $Date: 2005-11-17 03:59:05 $
    99*
    1010*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    1414#include "pslib_strict.h"
    1515
    16 //static psS32 testEOCParallax(void);
    17 static psS32 testEOCInit_Final(void);
    1816static psS32 testAberration(void);
    1917static psS32 testGravityDeflect(void);
    … …  
    2624
    2725testDescription tests[] = {
    28                               {testEOCInit_Final, 665, "psEOCInit_Final", 0, false},
    2926                              {testAberration, 666, "psAberration()", 0, false},
    3027                              {testGravityDeflect, 667, "psGravityDeflect()", 0, false},
    31                               //                              {testEOCParallax, 668, "psEOCParallax()", 0, false},
    3228                              {testEOCPrecession, 669, "psEOCPrecession()", 0, false},
    3329                              {testEOCPolar, 670, "psEOCPolar()", 0, false},
    … …  
    5652
    5753    //    return ( ! runTestSuite( stderr, "psEarthOrientation", tests, argc, argv ) );
    58 }
    59 
    60 psS32 testEOCInit_Final(void)
    61 {
    62     return 0;
    6354}
    6455
    … …  
    171162    psFree(UT1time);
    172163    //Check return values from valid precession input
    173 
    174164    pmodel = psEOC_PrecessionModel(time);
    175165    if ( pmodel == NULL ) {
    … …  
    288278psS32 testEOCNutation(void)
    289279{
     280    psTime *in = psTimeAlloc(PS_TIME_TAI);
     281    in->sec = 1131579114;
     282    in->nsec = 498489000;
     283    in->leapsecond = false;
     284    psTime *empty = NULL;
     285    psEarthPole *nute = NULL;
     286
     287    //Return NULL for NULL input time.
     288    nute = psEOC_NutationCorr(empty);
     289    if (nute != NULL) {
     290        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     291                "psEOC_NutationCorr failed to return NULL for NULL input time.\n");
     292        return 1;
     293    }
     294    //Return NULL for UT1 time input
     295    psTime *UT1time = psTimeAlloc(PS_TIME_UT1);
     296    nute = psEOC_NutationCorr(UT1time);
     297    if (nute != NULL) {
     298        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     299                "psEOC_NutationCorr failed to return NULL for UT1 input time.\n");
     300        return 2;
     301    }
     302    psFree(UT1time);
     303    //Check return values from valid nutation time input
     304    nute = psEOC_PrecessionModel(in);
     305    if ( nute == NULL ) {
     306        psError(PS_ERR_BAD_PARAMETER_NULL, false,
     307                "psEOC_NutationCorr returned NULL for valid input.\n");
     308        return 3;
     309    } else {
     310        printf("Nutation Correction output = x,y,s = %.8g, %.8g, %.8g\n\n",
     311               nute->x, nute->y, nute->s);
     312    }
     313    psFree(nute);
     314    psFree(in);
     315
     316    if (!p_psEOCFinalize() ) {
     317        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n");
     318        return 12;
     319    }
    290320
    291321    return 0;

  • trunk/psLib/test/astro/verified/tst_psEarthOrientation.stderr

    r5522 r5533  
    1 /***************************** TESTPOINT ******************************************\
    2 *             TestFile: tst_psEarthOrientation.c                                   *
    3 *            TestPoint: psEarthOrientation{psEOCInit_Final}                        *
    4 *             TestType: Positive                                                   *
    5 \**********************************************************************************/
    6 
    7 
    8 ---> TESTPOINT PASSED (psEarthOrientation{psEOCInit_Final} | tst_psEarthOrientation.c)
    9 
    101/***************************** TESTPOINT ******************************************\
    112*             TestFile: tst_psEarthOrientation.c                                   *
    … …  
    7061\**********************************************************************************/
    7162
     63<DATE><TIME>|<HOST>|E|psEOC_NutationCorr (FILE:LINENO)
     64    Unallowable operation: time is NULL.
     65<DATE><TIME>|<HOST>|E|psEOC_NutationCorr (FILE:LINENO)
     66    Invalid time input.  Time cannot be of type UT1.
    7267
    7368---> TESTPOINT PASSED (psEarthOrientation{psEOCNutation()} | tst_psEarthOrientation.c)

  • trunk/psLib/test/astro/verified/tst_psEarthOrientation.stdout

    r5524 r5533  
    1010 -- Eop = x=-0.00045805944, y=0.00013368871, s=-2.986086e-05
    1111 -- PolarMotion = x=0.25980922, y=0.26040097, s=-2.986086e-05
     12Nutation Correction output = x,y,s = 1.1557751e+08, 8865375.4, -5.1231901e+14
     13
    1214
    1315The Value of T is = 1128530000.931
Note: See TracChangeset for help on using the changeset viewer.