IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5771


Ignore:
Timestamp:
Dec 12, 2005, 3:31:54 PM (20 years ago)
Author:
drobbin
Message:

Updated tests and debugged code. PrecessionModel now working.

Location:
trunk/psLib
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r5749 r5771  
    88 *  @author Robert Daniel DeSonia, MHPCC
    99 *
    10  *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2005-12-08 02:49:04 $
     10 *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2005-12-13 01:31:54 $
    1212 *
    1313 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    432432    double S = psPolynomial1DEval(sPoly,t);
    433433
     434    X = SEC_TO_RAD(X * 1e-6);
     435    Y = SEC_TO_RAD(Y * 1e-6);
     436    S = SEC_TO_RAD(S * 1e-6);
     437
     438
    434439    // now calculate the non-poly portion from the tables
    435440
    … …  
    448453        double as = cols[1][lcv];
    449454        double ac = cols[2][lcv];
    450 
    451         X += as*tj*sin(arg) + ac*tj*cos(arg);
     455        //        printf("as-x, ac-x, = %.13g, %.13g\n", as, ac);
     456        as = SEC_TO_RAD(as) * 1e-6;
     457        ac = SEC_TO_RAD(ac) * 1e-6;
     458        //        X += as*tj*sin(arg) + ac*tj*cos(arg);
     459        X += (as*tj*sin(arg) + ac*cos(arg)) * tj;
    452460    }
    453461
    … …  
    465473        double as = cols[1][lcv];
    466474        double ac = cols[2][lcv];
    467 
    468         Y += as*tj*sin(arg) + ac*tj*cos(arg);
     475        as = SEC_TO_RAD(as) * 1e-6;
     476        ac = SEC_TO_RAD(ac) * 1e-6;
     477
     478        //        Y += as*tj*sin(arg) + ac*tj*cos(arg);
     479        Y += (as*tj*sin(arg) + ac*cos(arg)) * tj;
    469480    }
    470481
    … …  
    482493        double as = cols[1][lcv];
    483494        double ac = cols[2][lcv];
    484 
    485         S += as*tj*sin(arg) + ac*tj*cos(arg);
     495        as = SEC_TO_RAD(as) * 1e-6;
     496        ac = SEC_TO_RAD(ac) * 1e-6;
     497
     498        //        S += as*tj*sin(arg) + ac*tj*cos(arg);
     499        S += (as*tj*sin(arg) + ac*cos(arg)) * tj;
    486500    }
    487501
    … …  
    700714        A[2][2] = 1.0 - a*(pole->x*pole->x + pole->y*pole->y);
    701715    */
     716
    702717    A[0][0] = (1.0 - a*pole->x*pole->x)*cos(pole->s) - a*pole->x*pole->y*sin(pole->s);
    703718    A[1][0] = -a*pole->x*pole->y*cos(pole->s) + (1.0 - a*pole->y*pole->y)*sin(pole->s);
    … …  
    709724    A[1][2] = -pole->y;
    710725    A[2][2] = 1.0 - a*(pole->x*pole->x + pole->y*pole->y);
     726
    711727    double x, y, s;
    712728    x = pole->x;
    … …  
    10101026    psSphereRot *out = NULL;
    10111027
     1028    double x,y,s;
     1029    x = motion->x;
     1030    y = motion->y;
     1031    s = motion->s;
     1032
     1033
    10121034    //Setup Rotation Matrix for transformation (x,y,z rotation)
    10131035    //XXX: May need to be (z,y,x as in Mathworld?)
    1014     A[0][0] = cos(motion->x)*cos(-motion->s);
    1015     A[0][1] = sin(motion->y)*sin(motion->x)*cos(-motion->s) + cos(motion->y)*sin(-motion->s);
    1016     A[0][2] = -cos(motion->y)*sin(motion->x)*cos(-motion->s) + sin(motion->y)*sin(-motion->s);
    1017     A[1][0] = -cos(motion->x)*sin(-motion->s);
    1018     A[1][1] = -sin(motion->y)*sin(motion->x)*sin(-motion->s) + cos(motion->y)*cos(-motion->s);
    1019     A[1][2] = cos(motion->y)*sin(motion->x)*sin(-motion->s) + sin(motion->y)*cos(-motion->s);
    1020     A[2][0] = sin(motion->x);
    1021     A[2][1] = -sin(motion->y)*cos(motion->x);
    1022     A[2][2] = cos(motion->y)*cos(motion->x);
     1036    /*   A[0][0] = cos(motion->x)*cos(-motion->s);
     1037        A[0][1] = sin(motion->y)*sin(motion->x)*cos(-motion->s) + cos(motion->y)*sin(-motion->s);
     1038        A[0][2] = -cos(motion->y)*sin(motion->x)*cos(-motion->s) + sin(motion->y)*sin(-motion->s);
     1039        A[1][0] = -cos(motion->x)*sin(-motion->s);
     1040        A[1][1] = -sin(motion->y)*sin(motion->x)*sin(-motion->s) + cos(motion->y)*cos(-motion->s);
     1041        A[1][2] = cos(motion->y)*sin(motion->x)*sin(-motion->s) + sin(motion->y)*cos(-motion->s);
     1042        A[2][0] = sin(motion->x);
     1043        A[2][1] = -sin(motion->y)*cos(motion->x);
     1044        A[2][2] = cos(motion->y)*cos(motion->x);
     1045    */
     1046    /*
     1047        A[0][0] = cos(motion->x)*cos(-motion->s);
     1048        A[1][0] = sin(motion->y)*sin(motion->x)*cos(-motion->s) + cos(motion->y)*sin(-motion->s);
     1049        A[2][0] = -cos(motion->y)*sin(motion->x)*cos(-motion->s) + sin(motion->y)*sin(-motion->s);
     1050        A[0][1] = -cos(motion->x)*sin(-motion->s);
     1051        A[1][1] = -sin(motion->y)*sin(motion->x)*sin(-motion->s) + cos(motion->y)*cos(-motion->s);
     1052        A[2][1] = cos(motion->y)*sin(motion->x)*sin(-motion->s) + sin(motion->y)*cos(-motion->s);
     1053        A[0][2] = sin(motion->x);
     1054        A[1][2] = -sin(motion->y)*cos(motion->x);
     1055        A[2][2] = cos(motion->y)*cos(motion->x);
     1056    */
     1057    /*
     1058        A[0][0] = cos(s)*cos(x);
     1059        A[0][1] = sin(s)*sin(y) + cos(s)*sin(x)*sin(y);
     1060        A[0][2] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y);
     1061        A[1][0] = -sin(s)*cos(x);
     1062        A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y);
     1063        A[1][2] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y);
     1064        A[2][0] = sin(x);
     1065        A[2][1] = -cos(x)*sin(y);
     1066        A[2][2] = cos(x)*cos(y);
     1067    */
     1068    /*
     1069        A[0][0] = cos(s)*cos(x);
     1070        A[1][0] = sin(s)*sin(y) + cos(s)*sin(x)*sin(y);
     1071        A[2][0] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y);
     1072        A[0][1] = -sin(s)*cos(x);
     1073        A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y);
     1074        A[2][1] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y);
     1075        A[0][2] = sin(x);
     1076        A[1][2] = -cos(x)*sin(y);
     1077        A[2][2] = cos(x)*cos(y);
     1078    */
     1079
     1080    /*
     1081        psSphereRot r,s,t;
     1082     
     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);
     1088     
     1089        s.q0=0;
     1090        s.q1=sin(motion->x/2.0);
     1091        s.q2=0;
     1092        s.q3=cos(motion->x/2.0);
     1093     
     1094        t.q0=0;
     1095        t.q1=0;
     1096        t.q2=sin(-motion->s/2.0);
     1097        t.q3=cos(-motion->s/2.0);
     1098     
     1099        // calculate t*s*r.
     1100        psSphereRot* temp = psSphereRotCombine(NULL,&t,&s);
     1101        out = psSphereRotCombine(NULL, temp, &r);
     1102        psFree(temp);
     1103     
     1104        return out;
     1105    */
     1106
     1107    s = -s;
     1108    //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    */
     1119    /*    A[0][0] = cos(x)*cos(s);
     1120        A[0][1] = cos(x)*sin(s);
     1121        A[0][2] = -sin(x);
     1122        A[1][0] = sin(x)*sin(y)*cos(s) - cos(y)*sin(s);
     1123        A[1][1] = sin(x)*sin(y)*sin(s) + cos(y)*cos(s);
     1124        A[1][2] = cos(x)*sin(y);
     1125        A[2][0] = sin(x)*cos(y)*cos(s) + sin(y)*sin(s);
     1126        A[2][1] = sin(x)*cos(y)*sin(s) - sin(y)*cos(s);
     1127        A[2][2] = cos(x)*cos(y);
     1128    */
     1129    //New trial - mult z * y * x
     1130    /*    A[0][0] = cos(s)*cos(x);
     1131        A[1][0] = sin(s)*cos(y) + cos(s)*sin(x)*sin(y);
     1132        A[2][0] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y);
     1133        A[0][1] = -sin(s)*cos(x);
     1134        A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y);
     1135        A[2][1] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y);
     1136        A[0][2] = sin(x);
     1137        A[1][2] = -cos(x)*sin(y);
     1138        A[2][2] = cos(x)*cos(y);
     1139    */
     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
    10231150
    10241151    //Convert rotation matrix to quaternions
    10251152    out = rotMatrix_To_Quat(A);
     1153
     1154    //    out = psSphereRotAlloc(motion->y, motion->x, motion->s);
     1155
     1156    //    out->q0 = SEC_TO_RAD(out->q0);
     1157    //    out->q1 = SEC_TO_RAD(out->q1);
     1158    //    out->q2 = SEC_TO_RAD(out->q2) + (out->q0+out->q1/2.0);
     1159    //    out->q3 = sqrt(1.0 - (out->q0*out->q0 + out->q1*out->q1 + out->q2*out->q2) );
    10261160    /*    double diag_sum[3];
    10271161        int maxi;

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

    r5749 r5771  
    55*  @author d-Rob, MHPCC
    66*
    7 *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
    8 *  @date $Date: 2005-12-08 02:49:04 $
     7*  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
     8*  @date $Date: 2005-12-13 01:31:54 $
    99*
    1010*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    344344        yy = -0.0287618408203125;
    345345        ss = 0.0;
     346        xx = SEC_TO_RAD(xx) * 1e-6;
     347        yy = SEC_TO_RAD(yy) * 1e-6;
    346348        if ( fabs(pcorr->x - xx) > DBL_EPSILON || fabs(pcorr->y - yy) > DBL_EPSILON
    347349                || fabs(pcorr->s - ss) > DBL_EPSILON) {
    … …  
    599601    //    in->y += -0.0287618408203125;
    600602    //    in->s += 0.0;
     603    //    in->x += 3.05224300720406e-7;
     604    //    in->y += -1.39441339235822e-7;
     605    in->s += 0.0;
     606
    601607    double q0,q1,q2,q3;
    602608    q0 = -1.1984522406756289e-5;
    … …  
    619625    }
    620626
    621     //    psCube *tempCube = psCubeAlloc();
    622     //    tempCube->x = -0.35963388069046304;
    623     //    tempCube->y = 0.5555192509816625;
    624     //    tempCube->z = 0.749078321908413;
    625     //    obj = psCubeToSphere(tempCube);
    626     //    psFree(tempCube);
    627     obj = psSphereAlloc();
    628     obj->r = objR;
    629     obj->d = objD;
     627    if (fabs(rot->q0-q0) > DBL_EPSILON || fabs(rot->q1-q1) > DBL_EPSILON ||
     628            fabs(rot->q2-q2) > DBL_EPSILON || fabs(rot->q3-q3) > DBL_EPSILON) {
     629        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     630                "psSphereRot_CEOtoGCRS failed to return expected values.\n");
     631        printf("\n  Output sphere rotation   = %.13g,%.13g,%.13g,%.13g\n",
     632               rot->q0, rot->q1, rot->q2, rot->q3);
     633        printf("  Expected sphere rotation = %.13g,%.13g,%.13g,%.13g\n", q0,q1,q2,q3);
     634        printf("  a difference:   %.13g, %.13g, %.13g, %.13g \n", (rot->q0-q0), (rot->q1-q1),
     635               (rot->q2-q2), (rot->q3-q3) );
     636        printf("  abs difference:   %.13g, %.13g, %.13g, %.13g \n\n", (rot->q0+q0), (rot->q1+q1),
     637               (rot->q2+q2), (rot->q3+q3) );
     638        //        return 3;
     639    }
     640
     641
     642    psCube *tempCube = psCubeAlloc();
     643    tempCube->x = -0.35963388069046304;
     644    tempCube->y = 0.5555192509816625;
     645    tempCube->z = 0.7497078321908413;
     646    obj = psCubeToSphere(tempCube);
     647    psFree(tempCube);
     648    //    obj = psSphereAlloc();
     649    //    obj->r = objR;
     650    //    obj->d = objD;
    630651    psSphereRot *precessionNutation = psSphereRotConjugate(NULL, rot);
    631652    psSphere *result = psSphereRotApply(NULL, precessionNutation, obj);
    … …  
    639660    printf("A difference of:   %.13g, %.13g, %.13g\n\n",
    640661           (x-cube->x), (y-cube->y), (z-cube->z));
     662
     663    //great sphere difference between result and expected
     664    //Re = 6.38x10^6m
     665    psCube *expect = psCubeAlloc();
     666    expect->x = x;
     667    expect->y = y;
     668    expect->z = z;
     669    psSphere *expected = psCubeToSphere(expect);
     670    psFree(expect);
     671    double d = (6.38e6) * acos(cos(result->r)*cos(expected->r)*cos(result->d-expected->d) + sin(result->r)*sin(expected->r));
     672    printf("\n\nGreat circle difference of:  %.13g \n", d);
     673
     674    psFree(expected);
    641675    psFree(obj);
    642676
    643677
    644     if (fabs(rot->q0-q0) > DBL_EPSILON || fabs(rot->q1-q1) > DBL_EPSILON ||
    645             fabs(rot->q2-q2) > DBL_EPSILON || fabs(rot->q3-q3) > DBL_EPSILON) {
    646         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    647                 "psSphereRot_CEOtoGCRS failed to return expected values.\n");
    648         printf("\n  Output sphere rotation   = %.13g,%.13g,%.13g,%.13g\n",
    649                rot->q0, rot->q1, rot->q2, rot->q3);
    650         printf("  Expected sphere rotation = %.13g,%.13g,%.13g,%.13g\n", q0,q1,q2,q3);
    651         printf("  a difference:   %.13g, %.13g, %.13g, %.13g \n", (rot->q0-q0), (rot->q1-q1),
    652                (rot->q2-q2), (rot->q3-q3) );
    653         printf("  abs difference:   %.13g, %.13g, %.13g, %.13g \n\n", (rot->q0+q0), (rot->q1+q1),
    654                (rot->q2+q2), (rot->q3+q3) );
    655         return 3;
    656     }
    657678
    658679    psSphere *test = psSphereAlloc();
    … …  
    708729        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    709730                "psSphereRot_ITRStoTEO failed to return expected values.\n");
    710         printf("\n  Output sphere rotation   = %.13g,%.13g,%.13g,%.13g\n",
     731        printf("\n  Output sphere rotation   = %.13g, %.13g, %.13g, %.13g\n",
    711732               rot->q0, rot->q1, rot->q2, rot->q3);
    712         printf("  Expected sphere rotation = %.13g,%.13g,%.13g,%.13g\n", q0,q1,q2,q3);
    713         printf("  a difference:   %.13g, %.13g, %.13g, %.13g \n", (rot->q0-q0), (rot->q1-q1),
    714                (rot->q2-q2), (rot->q3-q3) );
     733        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) );
    715736        return 3;
    716737    }
Note: See TracChangeset for help on using the changeset viewer.