IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.