IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 19, 2005, 7:05:37 PM (21 years ago)
Author:
drobbin
Message:

Made several changes to EOC tests & fxns. Revamped GetPolarMotion to use interpolation.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r5793 r5814  
    88 *  @author Robert Daniel DeSonia, MHPCC
    99 *
    10  *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2005-12-15 02:37:48 $
     10 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2005-12-20 05:05:37 $
    1212 *
    1313 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    512512    S = SEC_TO_RAD(S * 1e-6);
    513513
    514 
    515514    // now calculate the non-poly portion from the tables
    516 
    517515    psF64* cols[17];
    518516    for (int lcv = 0; lcv < 17; lcv++) {
    … …  
    628626                out->x = cols[1][rowNum];
    629627                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;
     628                out->x = SEC_TO_RAD(out->x) * 1e-3;
     629                out->y = SEC_TO_RAD(out->y) * 1e-3;
    632630                rowNum = numRows;
    633631            } else {
    634632                out->x = cols[3][rowNum];
    635633                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;
     634                out->x = SEC_TO_RAD(out->x) * 1e-3;
     635                out->y = SEC_TO_RAD(out->y) * 1e-3;
    638636                rowNum = numRows;
    639637            }
    … …  
    762760    T += -2451545.0;
    763761    double theta = 2.0 * M_PI * (0.7790572732640 + 1.00273781191135448 * T);
    764     psSphereRot *out = psSphereRotAlloc(theta, 0.0, 0.0);
     762    psSphereRot *out = psSphereRotAlloc(-theta, 0.0, 0.0);
    765763    //    psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0);
    766764
    … …  
    773771{
    774772    PS_ASSERT_PTR_NON_NULL(time, NULL);
     773    PS_ASSERT_INT_WITHIN_RANGE(bulletin, PS_IERS_A, PS_IERS_B, NULL);
    775774
    776775    psEarthPole *out = psEarthPoleAlloc();
    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;
    783     //XXX:  Time must be converted to TAI before use?
    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     */
    791     //XXX: This may be the wrong idea... ADD says to use 3rd-order polys to interpolate
    792     //polar motion coordinates.
    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;
    803 
    804     //XXX: Apply polar tide correction here???
    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 
    814776    out->x = 0.0;
    815777    out->y = 0.0;
    … …  
    817779
    818780    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 
     781    /*
     782        if (MJD < 41684.0 || (MJD > 53434.0 && bulletin == PS_IERS_B) ||
     783                (MJD > 53858.0 && bulletin == PS_IERS_A) ) {
     784            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     785                    "Invalid time input.  Date, %lf, is out of range\n", MJD);
     786            return out;
     787        }
     788    */
    827789    // Check if EOC data loaded
    828790    if(! eocInitialized) {
    … …  
    841803    }
    842804    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;
     805    psVector *X = psVectorAlloc(numRows, PS_TYPE_F64);
     806    psVector *Y = psVectorAlloc(numRows, PS_TYPE_F64);
     807    psVector *S = psVectorAlloc(numRows, PS_TYPE_F64);
     808    psVector *T = psVectorAlloc(numRows, PS_TYPE_F64);
     809    if (bulletin == PS_IERS_A) {
     810        for (int rowNum = 0; rowNum < numRows; rowNum++) {
     811            T->data.F64[rowNum] = cols[0][rowNum];
     812            X->data.F64[rowNum] = cols[1][rowNum];
     813            Y->data.F64[rowNum] = cols[2][rowNum];
     814            S->data.F64[rowNum] = cols[3][rowNum];
     815        }
     816    } else {
     817        for (int rowNum = 0; rowNum < numRows; rowNum++) {
     818            T->data.F64[rowNum] = cols[0][rowNum];
     819            X->data.F64[rowNum] = cols[4][rowNum];
     820            Y->data.F64[rowNum] = cols[5][rowNum];
     821            S->data.F64[rowNum] = cols[6][rowNum];
     822        }
     823    }
     824
     825    double xOut = 0.0;
     826    double yOut = 0.0;
     827    double sOut = 0.0;
     828    double xTerm = 0.0;
     829    double yTerm = 0.0;
     830    double sTerm = 0.0;
     831    int k = 0;
     832    for (int i = 0; i < (numRows-1); i++) {
     833        if (MJD >= T->data.F64[i] && MJD < T->data.F64[i+1]) {
     834            k = i;
     835            if (k < 2) {
     836                k = 2;
    861837            }
    862         }
    863     }
    864 
     838            if (k > (numRows-2)) {
     839                k = numRows-2;
     840            }
     841            for (int m = k-1; m <= k+2; m++) {
     842                xTerm = X->data.F64[m];
     843                yTerm = Y->data.F64[m];
     844                sTerm = S->data.F64[m];
     845                for (int j = k-1; j <= k+2; j++) {
     846                    if ( m != j) {
     847                        double term = (MJD - T->data.F64[j])/(T->data.F64[m] - T->data.F64[j]);
     848                        xTerm *= term;
     849                        yTerm *= term;
     850                        sTerm *= term;
     851                    }
     852                }
     853                xOut += xTerm;
     854                yOut += yTerm;
     855                sOut += sTerm;
     856            }
     857            i = numRows-1;
     858        }
     859    }
     860    out->x = SEC_TO_RAD(xOut);
     861    out->y = SEC_TO_RAD(yOut);
     862    out->s = SEC_TO_RAD(sOut);
     863
     864
     865    /*
     866        for (int rowNum = 0; rowNum < numRows; rowNum++) {
     867            if ( (MJD - cols[0][rowNum]) < 1.0 ) {
     868                if (bulletin == PS_IERS_A) {
     869                    out->x = cols[1][rowNum];
     870                    out->y = cols[2][rowNum];
     871                    out->s = cols[3][rowNum];
     872                    out->x = SEC_TO_RAD(out->x);
     873                    out->y = SEC_TO_RAD(out->y);
     874                    out->s = SEC_TO_RAD(out->s);
     875                    rowNum = numRows;
     876                } else {
     877                    out->x = cols[4][rowNum];
     878                    out->y = cols[5][rowNum];
     879                    out->s = cols[6][rowNum];
     880                    out->x = SEC_TO_RAD(out->x);
     881                    out->y = SEC_TO_RAD(out->y);
     882                    out->s = SEC_TO_RAD(out->s);
     883                    rowNum = numRows;
     884                }
     885            }
     886        }
     887    */
     888    psFree(X);
     889    psFree(Y);
     890    psFree(S);
     891    psFree(T);
    865892    return out;
    866893}
    … …  
    10811108    y = motion->y;
    10821109    s = motion->s;
     1110    s = -s;
    10831111
    10841112
    … …  
    11301158
    11311159    /*
    1132             psSphereRot r,s,t;
     1160                psSphereRot r,p,t;
    11331161     
    1134             r.q0=sin(motion->y/2.0);
    1135             r.q1=0;
    1136             r.q2=0;
    1137             r.q3=cos(motion->y/2.0);
     1162                r.q0=sin(y/2.0);
     1163                r.q1=0;
     1164                r.q2=0;
     1165                r.q3=cos(y/2.0);
    11381166     
    1139             s.q0=0;
    1140             s.q1=sin(motion->x/2.0);
    1141             s.q2=0;
    1142             s.q3=cos(motion->x/2.0);
     1167                p.q0=0;
     1168                p.q1=sin(x/2.0);
     1169                p.q2=0;
     1170                p.q3=cos(x/2.0);
    11431171     
    1144             t.q0=0;
    1145             t.q1=0;
    1146             t.q2=sin(motion->s/2.0);
    1147             t.q3=cos(motion->s/2.0);
     1172                t.q0=0;
     1173                t.q1=0;
     1174                t.q2=sin(s/2.0);
     1175                t.q3=cos(s/2.0);
    11481176     
    1149             // calculate t*s*r.
    1150             psSphereRot* temp = psSphereRotCombine(NULL,&t,&s);
    1151             out = psSphereRotCombine(NULL, temp, &r);
    1152             psFree(temp);
     1177                // calculate t*s*r.
     1178                psSphereRot* temp = psSphereRotCombine(NULL,&t,&p);
     1179                out = psSphereRotCombine(NULL, temp, &r);
     1180                psFree(temp);
    11531181     
    1154             return out;
    1155     */
    1156 
    1157     s = -s;
     1182                return out;
     1183    */
     1184
    11581185    //Newest trial - mult. y matrix * x * z
     1186
    11591187    A[0][0] = cos(x)*cos(s);
    11601188    A[1][0] = cos(x)*sin(s);
    … …  
    11671195    A[2][2] = cos(x)*cos(y);
    11681196
    1169     /*    A[0][0] = cos(x)*cos(s);
    1170         A[0][1] = cos(x)*sin(s);
    1171         A[0][2] = -sin(x);
    1172         A[1][0] = sin(x)*sin(y)*cos(s) - cos(y)*sin(s);
    1173         A[1][1] = sin(x)*sin(y)*sin(s) + cos(y)*cos(s);
    1174         A[1][2] = cos(x)*sin(y);
    1175         A[2][0] = sin(x)*cos(y)*cos(s) + sin(y)*sin(s);
    1176         A[2][1] = sin(x)*cos(y)*sin(s) - sin(y)*cos(s);
    1177         A[2][2] = cos(x)*cos(y);
     1197    /*
     1198            A[0][0] = cos(x)*cos(s);
     1199            A[0][1] = cos(x)*sin(s);
     1200            A[0][2] = -sin(x);
     1201            A[1][0] = sin(x)*sin(y)*cos(s) - cos(y)*sin(s);
     1202            A[1][1] = sin(x)*sin(y)*sin(s) + cos(y)*cos(s);
     1203            A[1][2] = cos(x)*sin(y);
     1204            A[2][0] = sin(x)*cos(y)*cos(s) + sin(y)*sin(s);
     1205            A[2][1] = sin(x)*cos(y)*sin(s) - sin(y)*cos(s);
     1206            A[2][2] = cos(x)*cos(y);
    11781207    */
    11791208    //New trial - mult z * y * x
    … …  
    12031232    out = rotMatrix_To_Quat(A);
    12041233
    1205     //    out = psSphereRotAlloc(motion->y, motion->x, motion->s);
    1206 
    1207     //    out->q0 = SEC_TO_RAD(out->q0);
    1208     //    out->q1 = SEC_TO_RAD(out->q1);
    1209     //    out->q2 = SEC_TO_RAD(out->q2) + (out->q0+out->q1/2.0);
    1210     //    out->q3 = sqrt(1.0 - (out->q0*out->q0 + out->q1*out->q1 + out->q2*out->q2) );
    1211     /*    double diag_sum[3];
    1212         int maxi;
    1213         double recip;
    1214         diag_sum[0] = 1.0 + A[0][0] - A[1][1] - A[2][2];
    1215         diag_sum[1] = 1.0 - A[0][0] + A[1][1] - A[2][2];
    1216         diag_sum[2] = 1.0 - A[0][0] - A[1][1] + A[2][2];
    1217         diag_sum[3] = 1.0 + A[0][0] + A[1][1] + A[2][2];
    1218 
    1219         maxi = 0;
    1220         for (int i = 1; i < 4; ++i) {
    1221             if (diag_sum[i] > diag_sum[maxi]) {
    1222                 maxi = i;
    1223             }
    1224         }
    1225 
    1226         double p = 0.5 * sqrt(diag_sum[maxi]);
    1227         recip = 1.0 / (4.0 * p);
    1228 
    1229         if (maxi == 0) {
    1230             out->q0 = p;
    1231             out->q1 = recip * (A[0][1] + A[1][0]);
    1232             out->q2 = recip * (A[2][0] + A[0][2]);
    1233             out->q3 = recip * (A[1][2] - A[2][1]);
    1234         } else if (maxi == 1) {
    1235             out->q0 = recip * (A[0][1] + A[1][0]);
    1236             out->q1 = p;
    1237             out->q2 = recip * (A[1][2] + A[2][1]);
    1238             out->q3 = recip * (A[2][0] - A[0][2]);
    1239         } else if (maxi == 2) {
    1240             out->q0 = recip * (A[2][0] + A[0][2]);
    1241             out->q1 = recip * (A[1][2] + A[2][1]);
    1242             out->q2 = p;
    1243             out->q3 = recip * (A[0][1] - A[1][0]);
    1244         } else if (maxi == 3) {
    1245             out->q0 = recip * (A[1][2] - A[2][1]);
    1246             out->q1 = recip * (A[2][0] - A[0][2]);
    1247             out->q2 = recip * (A[0][1] - A[1][0]);
    1248             out->q3 = p;
    1249         }
    1250     */
    12511234    return out;
    12521235}
Note: See TracChangeset for help on using the changeset viewer.