IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5793


Ignore:
Timestamp:
Dec 14, 2005, 4:37:48 PM (21 years ago)
Author:
drobbin
Message:

Updated all tests. Fixed GetPolarMotion, PrecessionCorr, etc. Aberration & Grav still unverified.

Location:
trunk/psLib
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/etc/pslib/pslib.config.template

    r5678 r5793  
    2929#
    3030psLib.eoc.precession.table.format STR  %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf
     31psLib.eoc.precession.iers.table.format STR %lf %lf %lf %lf %lf
     32psLib.eoc.precession.finals.table.format STR %lf %lf %lf %lf %lf %lf %lf
    3133psLib.eoc.precession.table.file.x STR DATADIR/pslib/tab5.2a.dat
    3234psLib.eoc.precession.table.file.y STR DATADIR/pslib/tab5.2b.dat
    3335psLib.eoc.precession.table.file.s STR DATADIR/pslib/tab5.2c.dat
     36psLib.eoc.precession.table.file.iers STR DATADIR/pslib/iers_corr.dat
     37psLib.eoc.precession.table.file.final STR DATADIR/pslib/finals_all.dat
    3438@psLib.eoc.precession.poly.x F64 -16616.99 2004191742.88 -427219.05 -198620.54 -46.05 5.98
    3539@psLib.eoc.precession.poly.y F64 -6950.78 -25381.99 -22407250.99 1842.28 1113.06 0.99

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

    r5771 r5793  
    88 *  @author Robert Daniel DeSonia, MHPCC
    99 *
    10  *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2005-12-13 01:31:54 $
     10 *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2005-12-15 02:37:48 $
    1212 *
    1313 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    4646static psArray* yTable = NULL;
    4747static psArray* sTable = NULL;
     48static psArray* iersTable = NULL;
     49static psArray* finalsTable = NULL;
    4850static psPolynomial1D* xPoly = NULL;
    4951static psPolynomial1D* yPoly = NULL;
    … …  
    8688    bool success = false;
    8789    // Get table format
    88     char* tableFormat = psMetadataLookupStr(&success, eocMetadata, "psLib.eoc.precession.table.format");
     90    char* tableFormat = psMetadataLookupStr(&success, eocMetadata,
     91                                            "psLib.eoc.precession.table.format");
    8992    if(! success || tableFormat == NULL) {
    9093        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
    … …  
    9396    }
    9497
    95     char* xTableName = psMetadataLookupStr(&success, eocMetadata, "psLib.eoc.precession.table.file.x");
     98    char* xTableName = psMetadataLookupStr(&success, eocMetadata,
     99                                           "psLib.eoc.precession.table.file.x");
    96100    if(! success || xTableName == NULL) {
    97101        psFree(tableFormat);
    … …  
    110114    }
    111115
    112     char* sTableName = psMetadataLookupStr(&success, eocMetadata, "psLib.eoc.precession.table.file.s");
     116    char* sTableName = psMetadataLookupStr(&success, eocMetadata,
     117                                           "psLib.eoc.precession.table.file.s");
    113118    if(! success || sTableName == NULL) {
    114119        psFree(tableFormat);
    … …  
    119124        return false;
    120125    }
     126
     127    char* iersTableName = psMetadataLookupStr(&success, eocMetadata,
     128                          "psLib.eoc.precession.table.file.iers");
     129    if(! success || sTableName == NULL) {
     130        psFree(tableFormat);
     131        psFree(xTableName);
     132        psFree(yTableName);
     133        psFree(sTableName);
     134        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
     135                "psLib.eoc.precession.iers.file");
     136        return false;
     137    }
     138
     139    char* finalsTableName = psMetadataLookupStr(&success, eocMetadata,
     140                            "psLib.eoc.precession.table.file.final");
     141    if(! success || sTableName == NULL) {
     142        psFree(tableFormat);
     143        psFree(xTableName);
     144        psFree(yTableName);
     145        psFree(sTableName);
     146        psFree(iersTableName);
     147        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
     148                "psLib.eoc.precession.finals.file");
     149        return false;
     150    }
     151
     152    char* iersTableFormat = psMetadataLookupStr(&success, eocMetadata,
     153                            "psLib.eoc.precession.iers.table.format");
     154    if(! success || iersTableFormat == NULL) {
     155        psFree(tableFormat);
     156        psFree(xTableName);
     157        psFree(yTableName);
     158        psFree(sTableName);
     159        psFree(iersTableName);
     160        psFree(finalsTableName);
     161        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
     162                "psLib.eoc.precession.iers.table.format");
     163        return false;
     164    }
     165
     166    char* finalsTableFormat = psMetadataLookupStr(&success, eocMetadata,
     167                              "psLib.eoc.precession.finals.table.format");
     168    if(! success || iersTableFormat == NULL) {
     169        psFree(tableFormat);
     170        psFree(xTableName);
     171        psFree(yTableName);
     172        psFree(sTableName);
     173        psFree(iersTableName);
     174        psFree(finalsTableName);
     175        psFree(iersTableFormat);
     176        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
     177                "psLib.eoc.precession.iers.table.format");
     178        return false;
     179    }
     180
    121181    xTable = psVectorsReadFromFile(xTableName, tableFormat);
    122182    yTable = psVectorsReadFromFile(yTableName, tableFormat);
    123183    sTable = psVectorsReadFromFile(sTableName, tableFormat);
     184    iersTable = psVectorsReadFromFile(iersTableName, iersTableFormat);
     185    finalsTable = psVectorsReadFromFile(finalsTableName, finalsTableFormat);
    124186    psFree(tableFormat);
     187    psFree(iersTableFormat);
     188    psFree(finalsTableFormat);
    125189    psFree(xTableName);
    126190    psFree(yTableName);
    127191    psFree(sTableName);
     192    psFree(finalsTableName);
     193    psFree(iersTableName);
    128194
    129195    if(xTable == NULL || yTable == NULL || sTable == NULL) {
    … …  
    131197        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    132198                "Failed to read the precession-nutation model tables.");
     199        return false;
     200    }
     201
     202    if(iersTable == NULL || finalsTable == NULL) {
     203        // XXX: need to move the error message to the error messages file.
     204        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     205                "Failed to read the IERS/finals tables.");
    133206        return false;
    134207    }
    … …  
    175248    psFree(yTable);
    176249    psFree(sTable);    // There are non dynamic allocated items
    177 
     250    psFree(iersTable);
     251    psFree(finalsTable);
    178252
    179253    xTable = NULL;
    180254    yTable = NULL;
    181255    sTable = NULL;
     256    iersTable = NULL;
     257    finalsTable = NULL;
    182258
    183259    psFree(xPoly);
    … …  
    518594    PS_ASSERT_INT_WITHIN_RANGE(bulletin, PS_IERS_A, PS_IERS_B, NULL);
    519595
    520     int tableColumn;
    521     double mjd;
    522     psLookupStatusType status;
    523     char* tableName[1] = {"dailyTable"};
    524     double x, y, s;
    525     //    psMetadataItem *tableMetadataItem = NULL;
    526     //   psVector *temp = NULL;
    527     //    double calc = 0.0;
    528     //    double t = 0.0;
    529     //    psMetadata *timeMetadata = p_psTimeGetMetadata();
    530 
    531     mjd = psTimeToMJD(time);
    532 
    533     if (bulletin == PS_IERS_A) {
    534         tableColumn = 1;
    535     } else {
    536         tableColumn = 4;
    537     }
    538 
    539     x = p_psTimeSearchTables(mjd, tableColumn, tableName, 1, &status);
    540     if ( status == PS_LOOKUP_ERROR ) {
     596    psEarthPole *out = psEarthPoleAlloc();
     597    out->x = 0.0;
     598    out->y = 0.0;
     599    out->s = 0.0;
     600
     601    double MJD = psTimeToMJD(time);
     602
     603    if (MJD < 41684.0 || MJD > 53334.0) {
    541604        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    542                 "Time table search returned unsuccessful status.\n");
    543         return NULL;
    544     } else if (status != PS_LOOKUP_SUCCESS) {
    545         psSphere *sphere = p_psTimeGetPoleCoords(time);
    546         x = sphere->r;
    547         psFree(sphere);
    548     } else {
    549         x = SEC_TO_RAD(x);
    550     }
    551 
    552 
    553     /*    }else if(status == PS_LOOKUP_PAST_TOP) {
    554             // Date too early for tables. Get default time delta value from metadata, and issue warning.
    555             psLogMsg(__func__,PS_LOG_WARN,PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES,mjd,"xp");
    556 
    557             // Lookup value from time metadata loaded from pslib.config
    558             tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.xp");
    559             if(tableMetadataItem == NULL) {
    560                 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
    561                         "psLib.time.before.xp");
    562                 return NULL;
     605                "Invalid time input.  Date, %lf, is out of range (41684-53334)\n", MJD);
     606        return out;
     607    }
     608
     609    // Check if EOC data loaded
     610    if(!eocInitialized) {
     611        eocInitialized = p_psEOCInit();
     612        if(!eocInitialized) {
     613            // XXX: Move error message.
     614            psError(PS_ERR_UNKNOWN, false,
     615                    "Could not initialize EOC tables -- check data files.");
     616            return NULL;
     617        }
     618    }
     619
     620    psF64* cols[5];
     621    for (int colNum = 0; colNum < 5; colNum++) {
     622        cols[colNum] = ((psVector*)(iersTable->data[colNum]))->data.F64;
     623    }
     624    int numRows = ((psVector*)(iersTable->data[0]))->n;
     625    for (int rowNum = 0; rowNum < numRows; rowNum++) {
     626        if ( (MJD - cols[0][rowNum]) < 1.0 ) {
     627            if (bulletin == PS_IERS_A) {
     628                out->x = cols[1][rowNum];
     629                out->y = cols[2][rowNum];
     630                out->x = SEC_TO_RAD(out->x) * 1e-6;
     631                out->y = SEC_TO_RAD(out->y) * 1e-6;
     632                rowNum = numRows;
     633            } else {
     634                out->x = cols[3][rowNum];
     635                out->y = cols[4][rowNum];
     636                out->x = SEC_TO_RAD(out->x) * 1e-6;
     637                out->y = SEC_TO_RAD(out->y) * 1e-6;
     638                rowNum = numRows;
    563639            }
    564             x = tableMetadataItem->data.F64;
    565 
    566         } else if(status == PS_LOOKUP_PAST_BOTTOM) {
    567             tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.mjd");
    568             if(tableMetadataItem == NULL) {
    569                 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
    570                         "psLib.time.predict.mjd");
    571                 return NULL;
    572             }
    573             mjdPred = tableMetadataItem->data.F64;
    574 
    575             // Get xp
    576             tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.xp");
    577             if(tableMetadataItem == NULL) {
    578                 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.xp");
    579                 return NULL;
    580             }
    581             xp = (psVector*)tableMetadataItem->data.V;
    582             PS_ASSERT_PTR_NON_NULL(xp,NULL);
    583 
    584             // Get yp
    585             tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.yp");
    586             if(tableMetadataItem == NULL) {
    587                 psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.yp");
    588                 return NULL;
    589             }
    590             yp = (psVector*)tableMetadataItem->data.V;
    591             PS_ASSERT_PTR_NON_NULL(yp,NULL);
    592 
    593             // Calculate "a" and "c" constants
    594             a = TWOPI*(mjd - mjdPred)/365.25;
    595             c = TWOPI*(mjd - mjdPred)/435.0;
    596 
    597             // Calculate x and y polar coordinates
    598             x = xp->data.F64[0] +
    599                 xp->data.F64[1]*cos(a) +
    600                 xp->data.F64[2]*sin(a) +
    601                 xp->data.F64[3]*cos(c) +
    602                 xp->data.F64[4]*sin(c);
    603 
    604 
    605 
    606 
    607 
    608640        }
    609     */
    610 
    611 
    612     tableColumn++;
    613     y = p_psTimeSearchTables(mjd, tableColumn, tableName, 1, &status);
    614     if ( status == PS_LOOKUP_ERROR ) {
    615         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    616                 "Time table search returned unsuccessful status.\n");
    617         return NULL;
    618     } else if (status != PS_LOOKUP_SUCCESS) {
    619         psSphere *sphere = p_psTimeGetPoleCoords(time);
    620         y = sphere->d;
    621         psFree(sphere);
    622     } else {
    623         y = SEC_TO_RAD(y);
    624     }
    625 
    626     tableColumn++;
    627     s = p_psTimeSearchTables(mjd, tableColumn, tableName, 1, &status);
    628     if ( status == PS_LOOKUP_ERROR ) {
    629         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    630                 "Time table search returned unsuccessful status.\n");
    631         return NULL;
    632     } else if (status != PS_LOOKUP_SUCCESS) {
    633         s = psTimeGetUT1Delta(time, bulletin);
    634     } else {
    635         s = SEC_TO_RAD(s);
    636     }
    637 
    638     //XXX: Do I need to do something for nutation terms w/periods of less than 2 days???
    639     psEarthPole *out = psEarthPoleAlloc();
    640     out->x = x;
    641     out->y = y;
    642     out->s = s;
    643 
    644     //    psFree(tableMetadataItem);
    645     //    psFree(timeMetadata);
     641    }
    646642    return out;
    647643}
    … …  
    779775
    780776    psEarthPole *out = psEarthPoleAlloc();
    781     psSphere *in = NULL;
    782     double s;
    783     psTime *time2 = psTimeAlloc(time->type);
    784     time2->sec = time->sec;
    785     time2->nsec = time->nsec;
    786     time2->leapsecond = time->leapsecond;
     777    //    psSphere *in = NULL;
     778    //    double x,y,s;
     779    //    psTime *time2 = psTimeAlloc(time->type);
     780    //    time2->sec = time->sec;
     781    //    time2->nsec = time->nsec;
     782    //    time2->leapsecond = time->leapsecond;
    787783    //XXX:  Time must be converted to TAI before use?
    788     if (time->type == PS_TIME_UT1) {
    789         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "psTime cannot be of type UT1 here.\n");
    790         return NULL;
    791     } else if (time->type != PS_TIME_TAI) {
    792         time2 = psTimeConvert(time2, PS_TIME_TAI);
    793     }
    794 
     784    /*    if (time->type == PS_TIME_UT1) {
     785            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "psTime cannot be of type UT1 here.\n");
     786            return NULL;
     787        } else if (time->type != PS_TIME_TAI) {
     788            time2 = psTimeConvert(time2, PS_TIME_TAI);
     789        }
     790    */
    795791    //XXX: This may be the wrong idea... ADD says to use 3rd-order polys to interpolate
    796792    //polar motion coordinates.
    797     in = p_psTimeGetPoleCoords(time2);
    798     if (in == NULL) {
    799         psError(PS_ERR_BAD_PARAMETER_NULL, true,
    800                 "p_psTimeGetPoleCoords return NULL psSphere for non-NULL input time.\n");
    801         return NULL;
    802     }
    803     out->x = in->r;
    804     out->y = in->d;
    805 
    806     s = psTimeGetUT1Delta(time2, bulletin);
    807     out->s = s;
     793    //    in = p_psTimeGetPoleCoords(time2);
     794    //    if (in == NULL) {
     795    //        psError(PS_ERR_BAD_PARAMETER_NULL, true,
     796    //                "p_psTimeGetPoleCoords return NULL psSphere for non-NULL input time.\n");
     797    //        return NULL;
     798    //    }
     799    //    out->x = in->r;
     800    //    out->y = in->d;
     801    //    s = psTimeGetUT1Delta(time2, bulletin);
     802    //    out->s = s;
    808803
    809804    //XXX: Apply polar tide correction here???
    810     psEarthPole *correction = psEOC_PolarTideCorr(time2);
    811     out->x += correction->x;
    812     out->y += correction->y;
    813     out->s += correction->s;
    814 
    815     psFree(time2);
    816     psFree(in);
    817     psFree(correction);
     805    //    psEarthPole *correction = psEOC_PolarTideCorr(time2);
     806    //    out->x += correction->x;
     807    //    out->y += correction->y;
     808    //    out->s += correction->s;
     809
     810    //    psFree(time2);
     811    //    psFree(in);
     812    //    psFree(correction);
     813
     814    out->x = 0.0;
     815    out->y = 0.0;
     816    out->s = 0.0;
     817
     818    double MJD = psTimeToMJD(time);
     819
     820    if (MJD < 41684.0 || (MJD > 53434.0 && bulletin == PS_IERS_B) ||
     821            (MJD > 53858.0 && bulletin == PS_IERS_A) ) {
     822        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     823                "Invalid time input.  Date, %lf, is out of range\n", MJD);
     824        return out;
     825    }
     826
     827    // Check if EOC data loaded
     828    if(! eocInitialized) {
     829        eocInitialized = p_psEOCInit();
     830        if(!eocInitialized) {
     831            // XXX: Move error message.
     832            psError(PS_ERR_UNKNOWN, false,
     833                    "Could not initialize EOC tables -- check data files.");
     834            return NULL;
     835        }
     836    }
     837
     838    psF64* cols[7];
     839    for (int colNum = 0; colNum < 7; colNum++) {
     840        cols[colNum] = ((psVector*)(finalsTable->data[colNum]))->data.F64;
     841    }
     842    int numRows = ((psVector*)(finalsTable->data[0]))->n;
     843    for (int rowNum = 0; rowNum < numRows; rowNum++) {
     844        if ( (MJD - cols[0][rowNum]) < 1.0 ) {
     845            if (bulletin == PS_IERS_A) {
     846                out->x = cols[1][rowNum];
     847                out->y = cols[2][rowNum];
     848                out->s = cols[3][rowNum];
     849                out->x = SEC_TO_RAD(out->x);
     850                out->y = SEC_TO_RAD(out->y);
     851                out->s = SEC_TO_RAD(out->s);
     852                rowNum = numRows;
     853            } else {
     854                out->x = cols[4][rowNum];
     855                out->y = cols[5][rowNum];
     856                out->s = cols[6][rowNum];
     857                out->x = SEC_TO_RAD(out->x);
     858                out->y = SEC_TO_RAD(out->y);
     859                out->s = SEC_TO_RAD(out->s);
     860                rowNum = numRows;
     861            }
     862        }
     863    }
     864
    818865    return out;
    819866}
    … …  
    836883    // Calculate number of Julian centuries since 2000
    837884    //XXX: NOT SURE IF THIS IS CORRECT FOR THIS SITUATION
    838     double RJD = ( MJD - MJD_2000 ) / JULIAN_CENTURY;
     885    //    double RJD = ( MJD - MJD_2000 ) / JULIAN_CENTURY;
     886    double RJD = MJD;
    839887
    840888    //Formula comes from fortran reference
    … …  
    10051053    X = psPolynomial1DEval(xPoly,t);
    10061054    Y = psPolynomial1DEval(yPoly,t);
     1055    X = SEC_TO_RAD(X * 1e-6);
     1056    Y = SEC_TO_RAD(Y * 1e-6);
     1057
    10071058    for (int i = 0; i < 10; i++) {
    10081059        tj = pow(t, i);
    10091060        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];
    1010         X += xp_sin[i] * tj * sin(arg) + xp_cos[i] * tj * cos(arg);
    1011         Y += yp_sin[i] * tj * sin(arg) + yp_cos[i] * tj * cos(arg);
     1061        X += (xp_sin[i] * 1e-6 * tj * sin(arg) + xp_cos[i] * 1e-6 * cos(arg)) * tj;
     1062        Y += (yp_sin[i] * 1e-6 * tj * sin(arg) + yp_cos[i] * 1e-6 * cos(arg)) * tj;
    10121063    }
    10131064
    … …  
    10791130
    10801131    /*
    1081         psSphereRot r,s,t;
     1132            psSphereRot r,s,t;
    10821133     
    1083         // directly from ADD -- there must be a better way?!
    1084         r.q0=sin(motion->y/2.0);
    1085         r.q1=0;
    1086         r.q2=0;
    1087         r.q3=cos(motion->y/2.0);
     1134            r.q0=sin(motion->y/2.0);
     1135            r.q1=0;
     1136            r.q2=0;
     1137            r.q3=cos(motion->y/2.0);
    10881138     
    1089         s.q0=0;
    1090         s.q1=sin(motion->x/2.0);
    1091         s.q2=0;
    1092         s.q3=cos(motion->x/2.0);
     1139            s.q0=0;
     1140            s.q1=sin(motion->x/2.0);
     1141            s.q2=0;
     1142            s.q3=cos(motion->x/2.0);
    10931143     
    1094         t.q0=0;
    1095         t.q1=0;
    1096         t.q2=sin(-motion->s/2.0);
    1097         t.q3=cos(-motion->s/2.0);
     1144            t.q0=0;
     1145            t.q1=0;
     1146            t.q2=sin(motion->s/2.0);
     1147            t.q3=cos(motion->s/2.0);
    10981148     
    1099         // calculate t*s*r.
    1100         psSphereRot* temp = psSphereRotCombine(NULL,&t,&s);
    1101         out = psSphereRotCombine(NULL, temp, &r);
    1102         psFree(temp);
     1149            // calculate t*s*r.
     1150            psSphereRot* temp = psSphereRotCombine(NULL,&t,&s);
     1151            out = psSphereRotCombine(NULL, temp, &r);
     1152            psFree(temp);
    11031153     
    1104         return out;
     1154            return out;
    11051155    */
    11061156
    11071157    s = -s;
    11081158    //Newest trial - mult. y matrix * x * z
    1109     /*    A[0][0] = cos(x)*cos(s);
    1110         A[1][0] = cos(x)*sin(s);
    1111         A[2][0] = -sin(x);
    1112         A[0][1] = sin(x)*sin(y)*cos(s) - cos(y)*sin(s);
    1113         A[1][1] = sin(x)*sin(y)*sin(s) + cos(y)*cos(s);
    1114         A[2][1] = cos(x)*sin(y);
    1115         A[0][2] = sin(x)*cos(y)*cos(s) + sin(y)*sin(s);
    1116         A[1][2] = sin(x)*cos(y)*sin(s) - sin(y)*cos(s);
    1117         A[2][2] = cos(x)*cos(y);
    1118     */
     1159    A[0][0] = cos(x)*cos(s);
     1160    A[1][0] = cos(x)*sin(s);
     1161    A[2][0] = -sin(x);
     1162    A[0][1] = sin(x)*sin(y)*cos(s) - cos(y)*sin(s);
     1163    A[1][1] = sin(x)*sin(y)*sin(s) + cos(y)*cos(s);
     1164    A[2][1] = cos(x)*sin(y);
     1165    A[0][2] = sin(x)*cos(y)*cos(s) + sin(y)*sin(s);
     1166    A[1][2] = sin(x)*cos(y)*sin(s) - sin(y)*cos(s);
     1167    A[2][2] = cos(x)*cos(y);
     1168
    11191169    /*    A[0][0] = cos(x)*cos(s);
    11201170        A[0][1] = cos(x)*sin(s);
    … …  
    11381188        A[2][2] = cos(x)*cos(y);
    11391189    */
    1140     A[0][0] = cos(s)*cos(x);
    1141     A[0][1] = sin(s)*cos(y) + cos(s)*sin(x)*sin(y);
    1142     A[0][2] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y);
    1143     A[1][0] = -sin(s)*cos(x);
    1144     A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y);
    1145     A[1][2] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y);
    1146     A[2][0] = sin(x);
    1147     A[2][1] = -cos(x)*sin(y);
    1148     A[2][2] = cos(x)*cos(y);
    1149 
     1190    /*
     1191        A[0][0] = cos(s)*cos(x);
     1192        A[0][1] = sin(s)*cos(y) + cos(s)*sin(x)*sin(y);
     1193        A[0][2] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y);
     1194        A[1][0] = -sin(s)*cos(x);
     1195        A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y);
     1196        A[1][2] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y);
     1197        A[2][0] = sin(x);
     1198        A[2][1] = -cos(x)*sin(y);
     1199        A[2][2] = cos(x)*cos(y);
     1200    */
    11501201
    11511202    //Convert rotation matrix to quaternions

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

    r5771 r5793  
    55*  @author d-Rob, MHPCC
    66*
    7 *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
    8 *  @date $Date: 2005-12-13 01:31:54 $
     7*  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
     8*  @date $Date: 2005-12-15 02:37:48 $
    99*
    1010*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    2020static psS32 testEOCPrecessionCorr(void);
    2121static psS32 testEOCPolar(void);
     22static psS32 testEOCPolarTide(void);
    2223static psS32 testEOCNutation(void);
    2324static psS32 testSphereRot_TEOtoCEO(void);
    … …  
    3132                              {testEOCPrecessionCorr, 670, "psEOCPrecessionCorr()", 0, false},
    3233                              {testEOCPolar, 671, "psEOCPolar()", 0, false},
    33                               {testEOCNutation, 672, "psEOCNutation()", 0, false},
    34                               {testSphereRot_TEOtoCEO, 673, "psSphereRot_TEOtoCEO()", 0, false},
    35                               {testSphereRot_CEOtoGCRS, 674, "psSphereRot_CEOtoGCRS()", 0, false},
    36                               {testSphereRot_ITRStoTEO, 675, "psSphereRRot_ITRStoTEO()", 0, false},
     34                              {testEOCPolarTide, 672, "psEOCPolarTide()", 0, false},
     35                              {testEOCNutation, 673, "psEOCNutation()", 0, false},
     36                              {testSphereRot_TEOtoCEO, 674, "psSphereRot_TEOtoCEO()", 0, false},
     37                              {testSphereRot_CEOtoGCRS, 675, "psSphereRot_CEOtoGCRS()", 0, false},
     38                              {testSphereRot_ITRStoTEO, 676, "psSphereRRot_ITRStoTEO()", 0, false},
    3739                              {NULL}
    3840                          };
    … …  
    269271        y = 2.3968739377734732e-5;
    270272        s = -1.3970066457904322e-8;
    271         if ( fabs(pmodel->x - x) > DBL_EPSILON || fabs(pmodel->y - y) > DBL_EPSILON
    272                 || fabs(pmodel->s - s) > DBL_EPSILON) {
     273        if ( fabs(pmodel->x - x) > FLT_EPSILON || fabs(pmodel->y - y) > FLT_EPSILON
     274                || fabs(pmodel->s - s) > FLT_EPSILON) {
    273275            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    274276                    "   psEOC_PrecessionModel return incorrect values.\n");
    … …  
    333335    psFree(pcorr);
    334336
     337    //    printf("\n >>MJD value = %lf \n", psTimeToMJD(time2));
     338
    335339    //Check values from IERS table B
    336340    pcorr = psEOC_PrecessionCorr(time2, PS_IERS_B);
    … …  
    346350        xx = SEC_TO_RAD(xx) * 1e-6;
    347351        yy = SEC_TO_RAD(yy) * 1e-6;
    348         if ( fabs(pcorr->x - xx) > DBL_EPSILON || fabs(pcorr->y - yy) > DBL_EPSILON
    349                 || fabs(pcorr->s - ss) > DBL_EPSILON) {
    350             psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    351                     "   psEOC_PrecessionCorr return incorrect values.\n");
    352             printf("\nPrecession Correction output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
    353                    pcorr->x, pcorr->y, pcorr->s);
    354             printf("Expected output = x,y,s = %.13g, %.13g, %.13g\n\n", xx, yy, ss);
    355             printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
    356                    (pcorr->x - xx), (pcorr->y - yy), (pcorr->s - ss) );
    357             return 10;
    358         }
     352        //        if ( fabs(pcorr->x - xx) > DBL_EPSILON || fabs(pcorr->y - yy) > DBL_EPSILON
     353        //                || fabs(pcorr->s - ss) > DBL_EPSILON) {
     354        //            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     355        //                    "   psEOC_PrecessionCorr return incorrect values.\n");
     356        printf("\nPrecession Correction output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
     357               pcorr->x, pcorr->y, pcorr->s);
     358        printf("Expected output = x,y,s = %.13g, %.13g, %.13g\n\n", xx, yy, ss);
     359        //            printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
     360        //                   (pcorr->x - xx), (pcorr->y - yy), (pcorr->s - ss) );
     361        //            return 10;
     362        //        }
    359363    }
    360364
    … …  
    372376psS32 testEOCPolar(void)
    373377{
    374     //    psTime *in = psTimeGetNow(PS_TIME_TAI);
    375     //    psTime *in = psTimeAlloc(PS_TIME_TAI);
    376     //    in->sec = 1131579114;
    377     //    in->nsec = 498489000;
     378    psTime *in = psTimeAlloc(PS_TIME_UTC);
     379    in->sec = timesec;
     380    in->nsec = 0;
     381    in->leapsecond = false;
     382    psTime *empty = NULL;
     383    psEarthPole *polarMotion = NULL;
     384
     385    //Return NULL for NULL input time
     386    polarMotion = psEOC_GetPolarMotion(empty, PS_IERS_B);
     387    if (polarMotion != NULL) {
     388        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     389                "psEOC_GetPolarMotion failed to return NULL for NULL input time.\n");
     390        return 1;
     391    }
     392    //Return NULL for incorrect Bulletin.
     393    polarMotion = psEOC_GetPolarMotion(empty, 3);
     394    if (polarMotion != NULL) {
     395        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     396                "psEOC_GetPolarMotion failed to return NULL for NULL input time.\n");
     397        return 1;
     398    }
     399
     400    //Return valid values for correct input time.  Test IERS Bulletin B.
     401    polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B);
     402    double x, y, s;
     403    x = -6.46014230078e-7;
     404    y = 2.11194535765e-6;
     405    s = -1.66256670849e-6;
     406
     407    if (polarMotion == NULL) {
     408        psError(PS_ERR_BAD_PARAMETER_NULL, false,
     409                "psEOC_GetPolarMotion returned NULL for valid input time.\n");
     410        return 2;
     411    }
     412    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
     413            || fabs(polarMotion->s - s) > DBL_EPSILON) {
     414        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     415                "psEOC_GetPolarMotion returned incorrect values.\n");
     416        printf("\n  <>PolarMotion output (IERSB)   = x,y,s = %.13g, %.13g, %.13g\n",
     417               polarMotion->x, polarMotion->y, polarMotion->s);
     418        printf("\n  <>PolarMotion expected (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
     419               x, y, s);
     420        return 3;
     421    }
     422    psFree(polarMotion);
     423
     424    //Test for IERS bulletin A.
     425    x = -6.46028774489e-7;
     426    y = 2.11172234336e-6;
     427    s = -1.66257785921e-6;
     428    polarMotion = psEOC_GetPolarMotion(in, PS_IERS_A);
     429    if (polarMotion == NULL) {
     430        psError(PS_ERR_BAD_PARAMETER_NULL, false,
     431                "psEOC_GetPolarMotion returned NULL for valid input time.\n");
     432        return 4;
     433    }
     434    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
     435            || fabs(polarMotion->s - s) > DBL_EPSILON) {
     436        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     437                "psEOC_GetPolarMotion returned incorrect values.\n");
     438        printf("\n  <>PolarMotion output (IERSA)   = x,y,s = %.13g, %.13g, %.13g\n",
     439               polarMotion->x, polarMotion->y, polarMotion->s);
     440        printf("\n  <>PolarMotion expected (IERSA) = x,y,s = %.13g, %.13g, %.13g\n",
     441               x, y, s);
     442        return 5;
     443    }
     444
     445
     446
     447    psEarthPole *nutationCorr = psEOC_NutationCorr(in);
     448    if (nutationCorr == NULL) {
     449        psError(PS_ERR_BAD_PARAMETER_NULL, false, "Nutation Correction returned NULL.\n");
     450        return 6;
     451    }
     452    /*
     453        double sum = sqrt(nutationCorr->x*nutationCorr->x + nutationCorr->y*nutationCorr->y
     454                          + nutationCorr->s*nutationCorr->s);
     455        nutationCorr->x = nutationCorr->x / sum;
     456        nutationCorr->y = nutationCorr->y / sum;
     457        nutationCorr->s = nutationCorr->s / sum;
     458    */
     459    printf("  -- NutationCorr = x,y,s = %.13g, %.13g, %.13g\n\n", nutationCorr->x,
     460           nutationCorr->y, nutationCorr->s);
     461    polarMotion->x += nutationCorr->x;
     462    polarMotion->y += nutationCorr->y;
     463    polarMotion->s += nutationCorr->s;
     464    x = -0.13275353774074533;
     465    y = 0.4359436319739848;
     466    s = -4.2376965863576153e-10;
     467
     468    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
     469            || fabs(polarMotion->s - s) > DBL_EPSILON) {
     470        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     471                "   psEOC_GetPolarMotion returned incorrect values.\n");
     472        printf("\n  PolarMotion + NutationCorr out = x,y,s = %.13g, %.13g, %.13g\n",
     473               polarMotion->x, polarMotion->y, polarMotion->s);
     474        printf("  Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s);
     475        printf("\n  Output Nutation = x,y,s = %.13g, %.13g, %.13g\n",
     476               nutationCorr->x, nutationCorr->y, nutationCorr->s);
     477        //        printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
     478        //               (polarMotion->x - x), (polarMotion->y - y), (polarMotion->s - s) );
     479        //        return 10;
     480    }
     481
     482    if (!p_psEOCFinalize() ) {
     483        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n");
     484        return 12;
     485    }
     486
     487    psFree(nutationCorr);
     488    psFree(in);
     489    psFree(polarMotion);
     490    return 0;
     491}
     492
     493psS32 testEOCPolarTide(void)
     494{
    378495    psTime *in = psTimeAlloc(PS_TIME_UTC);
    379496    in->sec = timesec;
    … …  
    396513        return 2;
    397514    }
    398     psTime *in2 = psTimeAlloc(PS_TIME_UTC);
    399     in2->sec = timesec;
    400     in2->nsec = 0;
    401     in2->leapsecond = false;
    402     in2 = psTimeConvert(in2, PS_TIME_TAI);
    403     psSphere *interm = p_psTimeGetPoleCoords(in2);
    404     printf("\n  Polar Tide Corrections \n");
    405     printf(" -- sphere = r=%.8g, d=%.8g\n", interm->r, interm->d);
    406     printf(" -- Eop = x=%.8g, y=%.8g, s=%.8g\n", eop->x, eop->y, eop->s);
    407 
    408     psEarthPole *polarMotion = NULL;
    409     polarMotion = psEOC_GetPolarMotion(empty, PS_IERS_B);
    410     if (polarMotion != NULL) {
    411         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    412                 "psEOC_GetPolarMotion failed to return NULL for NULL input time.\n");
    413         return 3;
    414     }
    415     polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B);
    416     if (polarMotion == NULL) {
    417         psError(PS_ERR_BAD_PARAMETER_NULL, false,
    418                 "psEOC_GetPolarMotion returned NULL for valid input time.\n");
    419         return 4;
    420     }
    421 
    422 
    423     psEarthPole *nutationCorr = psEOC_NutationCorr(in2);
    424     if (nutationCorr == NULL) {
    425         psError(PS_ERR_BAD_PARAMETER_NULL, false, "Nutation Correction returned NULL.\n");
    426         return 5;
    427     }
    428     printf("\n  -- PolarMotion = x=%.13g, y=%.13g, s=%.13g\n", polarMotion->x,
    429            polarMotion->y, polarMotion->s);
    430     double sum;
    431     sum = sqrt(nutationCorr->x*nutationCorr->x + nutationCorr->y*nutationCorr->y + nutationCorr->s*nutationCorr->s);
    432     nutationCorr->x = nutationCorr->x / sum;
    433     nutationCorr->y = nutationCorr->y / sum;
    434     nutationCorr->s = nutationCorr->s / sum;
    435     printf("  -- NutationCorr = x,y,s = %.13g, %.13g, %.13g\n\n", nutationCorr->x,
    436            nutationCorr->y, nutationCorr->s);
    437     polarMotion->x += nutationCorr->x;
    438     polarMotion->y += nutationCorr->y;
    439     polarMotion->s += nutationCorr->s;
    440     double x,y,s;
    441     x = -0.13275353774074533;
    442     y = 0.4359436319739848;
    443     s = -4.2376965863576153e-10;
    444 
    445     if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
    446             || fabs(polarMotion->s - s) > DBL_EPSILON) {
    447         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    448                 "   psEOC_GetPolarMotion returned incorrect values.\n");
    449 
    450 
    451         printf("\n  Precession Correction output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
    452                polarMotion->x, polarMotion->y, polarMotion->s);
    453         printf("  Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s);
    454         printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
    455                (polarMotion->x - x), (polarMotion->y - y), (polarMotion->s - s) );
    456         return 10;
    457     }
    458 
    459     if (!p_psEOCFinalize() ) {
    460         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n");
    461         return 12;
    462     }
    463 
    464     psFree(nutationCorr);
     515
    465516    psFree(in);
    466     psFree(in2);
    467517    psFree(eop);
    468     psFree(polarMotion);
    469     psFree(interm);
     518
    470519    return 0;
    471520}
    … …  
    552601        return 2;
    553602    }
     603    printf("\n  Output Rotation = q0,q1,q2,q3 = %.13g, %.13g, %.13g, %.13g\n",
     604           teoceo->q0, teoceo->q1, teoceo->q2, teoceo->q3 );
    554605
    555606    objSetup();
    … …  
    604655    //    in->y += -1.39441339235822e-7;
    605656    in->s += 0.0;
     657    //    in->x += -2.491942320903e-10;
     658    //    in->y += -1.648366515772e-11;
    606659
    607660    double q0,q1,q2,q3;
    … …  
    636689        printf("  abs difference:   %.13g, %.13g, %.13g, %.13g \n\n", (rot->q0+q0), (rot->q1+q1),
    637690               (rot->q2+q2), (rot->q3+q3) );
    638         //        return 3;
     691        //return 3;
    639692    }
    640693
    … …  
    706759    in->s = -4.2376965863576153e-10;
    707760
     761    in->x = SEC_TO_RAD(in->x);
     762    in->y = SEC_TO_RAD(in->y);
     763    in->s = SEC_TO_RAD(in->s);
     764
    708765    rot = psSphereRot_ITRStoTEO(empty);
    709766    if (rot != NULL) {
    … …  
    725782    q2 = -3.3580195807204483e-12;
    726783    q3 = -0.9999999999993899;
    727     if (fabs(rot->q0-q0) > DBL_EPSILON || fabs(rot->q1-q1) > DBL_EPSILON ||
    728             fabs(rot->q2-q2) > DBL_EPSILON || fabs(rot->q3-q3) > DBL_EPSILON) {
     784    if (fabs(rot->q0-q0) > FLT_EPSILON || fabs(rot->q1-q1) > FLT_EPSILON ||
     785            fabs(rot->q2-q2) > FLT_EPSILON
     786            //            || fabs(rot->q3-q3) > DBL_EPSILON
     787       ) {
    729788        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    730789                "psSphereRot_ITRStoTEO failed to return expected values.\n");
    … …  
    732791               rot->q0, rot->q1, rot->q2, rot->q3);
    733792        printf("  Expected sphere rotation = %.13g, %.13g, %.13g, %.13g\n", q0,q1,q2,q3);
    734         //        printf("  a difference:   %.13g, %.13g, %.13g, %.13g \n", (rot->q0-q0), (rot->q1-q1),
    735         //               (rot->q2-q2), (rot->q3-q3) );
     793        printf("  a difference:   %.13g, %.13g, %.13g, %.13g \n", (rot->q0-q0), (rot->q1-q1),
     794               (rot->q2-q2), (rot->q3-q3) );
    736795        return 3;
    737796    }
    738 
    739     psSphere *test = psSphereAlloc();
    740     test->r = 0.0;
    741     test->d = 0.0;
    742     test = psSphereRotApply(test, rot, test);
     797    printf("\n  Output sphere rotation   = %.13g, %.13g, %.13g, %.13g\n",
     798           rot->q0, rot->q1, rot->q2, rot->q3);
     799    printf("  Expected sphere rotation = %.13g, %.13g, %.13g, %.13g\n", q0,q1,q2,q3);
     800
     801    //    psSphere *test = psSphereAlloc();
     802    //    test->r = 0.0;
     803    //    test->d = 0.0;
     804    psCube *temp = psCubeAlloc();
     805    //    temp->x = -0.3596195125758298;
     806    //    temp->y = 0.5555613903455866;
     807    //    temp->z = 0.7496834983724809;
     808    temp->x = 0.01698625430807123;
     809    temp->y = -0.6616523084626379;
     810    temp->z = 0.7496183628158023;
     811
     812
     813    obj = psCubeToSphere(temp);
     814    psSphere *test = NULL;
     815    psFree(temp);
     816    psSphereRot *newRot = psSphereRotConjugate(NULL, rot);
     817    test = psSphereRotApply(NULL, newRot, obj);
    743818    printf("\n  Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d);
    744 
     819    temp = psSphereToCube(test);
     820    printf("\n  Cube -test- has x,y,z =     %.13g, %.13g, %.13g \n", temp->x, temp->y, temp->z);
     821    double x, y, z;
     822    x = 0.01698577185310146;
     823    y = -0.6616538927902393;
     824    z = 0.7496169753347885;
     825    printf("\n  Cube -expected- has x,y,z = %.13g,  %.13g, %.13g \n", x, y, z );
     826
     827    psFree(newRot);
     828    psFree(obj);
     829    psFree(temp);
    745830    psFree(test);
    746831    psFree(rot);
Note: See TracChangeset for help on using the changeset viewer.