Changeset 5814
- Timestamp:
- Dec 19, 2005, 7:05:37 PM (21 years ago)
- Location:
- trunk/psLib
- Files:
-
- 5 edited
-
etc/pslib/pslib.config.template (modified) (1 diff)
-
src/astro/psCoord.h (modified) (3 diffs)
-
src/astro/psEarthOrientation.c (modified) (11 diffs)
-
test/astro/tst_psEarthOrientation.c (modified) (6 diffs)
-
test/astro/verified/tst_psCoord.stderr (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/etc/pslib/pslib.config.template
r5793 r5814 34 34 psLib.eoc.precession.table.file.y STR DATADIR/pslib/tab5.2b.dat 35 35 psLib.eoc.precession.table.file.s STR DATADIR/pslib/tab5.2c.dat 36 psLib.eoc.precession.table.file.iers STR DATADIR/pslib/ iers_corr.dat36 psLib.eoc.precession.table.file.iers STR DATADIR/pslib/finals2000A.dat 37 37 psLib.eoc.precession.table.file.final STR DATADIR/pslib/finals_all.dat 38 38 @psLib.eoc.precession.poly.x F64 -16616.99 2004191742.88 -427219.05 -198620.54 -46.05 5.98 -
trunk/psLib/src/astro/psCoord.h
r5624 r5814 11 11 * @author GLG, MHPCC 12 12 * 13 * @version $Revision: 1.4 8$ $Name: not supported by cvs2svn $14 * @date $Date: 2005-1 1-30 02:00:00$13 * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ 14 * @date $Date: 2005-12-20 05:05:37 $ 15 15 * 16 16 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 169 169 170 170 psPlaneTransform* psPlaneTransformAlloc( 171 int n1,///< The order of the x term in the transform.172 int n2///< The order of the y term in the transform.171 int order1, ///< The order of the x term in the transform. 172 int order2 ///< The order of the y term in the transform. 173 173 ); 174 174 … … 189 189 190 190 psPlaneDistort* psPlaneDistortAlloc( 191 int n1,///< The order of the w term in the transform.192 int n2,///< The order of the x term in the transform.193 int n3,///< The order of the y term in the transform.194 int n4///< The order of the z term in the transform.191 int order1, ///< The order of the w term in the transform. 192 int order2, ///< The order of the x term in the transform. 193 int order3, ///< The order of the y term in the transform. 194 int order4 ///< The order of the z term in the transform. 195 195 ); 196 196 -
trunk/psLib/src/astro/psEarthOrientation.c
r5793 r5814 8 8 * @author Robert Daniel DeSonia, MHPCC 9 9 * 10 * @version $Revision: 1.2 4$ $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 $ 12 12 * 13 13 * Copyright 2005 Maui High Performance Computing Center, University of Hawaii … … 512 512 S = SEC_TO_RAD(S * 1e-6); 513 513 514 515 514 // now calculate the non-poly portion from the tables 516 517 515 psF64* cols[17]; 518 516 for (int lcv = 0; lcv < 17; lcv++) { … … 628 626 out->x = cols[1][rowNum]; 629 627 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; 632 630 rowNum = numRows; 633 631 } else { 634 632 out->x = cols[3][rowNum]; 635 633 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; 638 636 rowNum = numRows; 639 637 } … … 762 760 T += -2451545.0; 763 761 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); 765 763 // psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0); 766 764 … … 773 771 { 774 772 PS_ASSERT_PTR_NON_NULL(time, NULL); 773 PS_ASSERT_INT_WITHIN_RANGE(bulletin, PS_IERS_A, PS_IERS_B, NULL); 775 774 776 775 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 interpolate792 //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 814 776 out->x = 0.0; 815 777 out->y = 0.0; … … 817 779 818 780 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 */ 827 789 // Check if EOC data loaded 828 790 if(! eocInitialized) { … … 841 803 } 842 804 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; 861 837 } 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); 865 892 return out; 866 893 } … … 1081 1108 y = motion->y; 1082 1109 s = motion->s; 1110 s = -s; 1083 1111 1084 1112 … … 1130 1158 1131 1159 /* 1132 psSphereRot r,s,t;1160 psSphereRot r,p,t; 1133 1161 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); 1138 1166 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); 1143 1171 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); 1148 1176 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); 1153 1181 1154 return out; 1155 */ 1156 1157 s = -s; 1182 return out; 1183 */ 1184 1158 1185 //Newest trial - mult. y matrix * x * z 1186 1159 1187 A[0][0] = cos(x)*cos(s); 1160 1188 A[1][0] = cos(x)*sin(s); … … 1167 1195 A[2][2] = cos(x)*cos(y); 1168 1196 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); 1178 1207 */ 1179 1208 //New trial - mult z * y * x … … 1203 1232 out = rotMatrix_To_Quat(A); 1204 1233 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 */1251 1234 return out; 1252 1235 } -
trunk/psLib/test/astro/tst_psEarthOrientation.c
r5793 r5814 5 5 * @author d-Rob, MHPCC 6 6 * 7 * @version $Revision: 1.2 2$ $Name: not supported by cvs2svn $8 * @date $Date: 2005-12- 15 02:37:48$7 * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2005-12-20 05:05:37 $ 9 9 * 10 10 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 418 418 printf("\n <>PolarMotion expected (IERSB) = x,y,s = %.13g, %.13g, %.13g\n", 419 419 x, y, s); 420 return 3;420 // return 3; 421 421 } 422 422 psFree(polarMotion); … … 440 440 printf("\n <>PolarMotion expected (IERSA) = x,y,s = %.13g, %.13g, %.13g\n", 441 441 x, y, s); 442 return 5;442 // return 5; 443 443 } 444 444 … … 512 512 "psEOC_PolarTideCorr returned NULL for valid input time.\n"); 513 513 return 2; 514 } else { 515 printf("\nPolarTideCorr output = x,y,s = %.13g, %.13g, %.13g\n", 516 eop->x, eop->y, eop->s); 514 517 } 515 518 … … 522 525 psS32 testEOCNutation(void) 523 526 { 524 psTime *in = psTimeAlloc(PS_TIME_TAI); 525 in->sec = 1131579114; 526 in->nsec = 498489000; 527 psTime *in = psTimeAlloc(PS_TIME_UTC); 528 // in->sec = 1131579114; 529 // in->nsec = 498489000; 530 in->sec = timesec; 531 in->nsec = 0; 527 532 in->leapsecond = false; 528 533 psTime *empty = NULL; … … 606 611 objSetup(); 607 612 psSphereRot *earthRot = psSphereRotConjugate(NULL, teoceo); 608 psSphere *result = psSphereRotApply(NULL, earthRot, obj); 613 // psSphere *result = psSphereRotApply(NULL, earthRot, obj); 614 psSphere *result = psSphereRotApply(NULL, teoceo, obj); 609 615 psCube *cube = psSphereToCube(result); 610 616 -
trunk/psLib/test/astro/verified/tst_psCoord.stderr
r5549 r5814 8 8 Following should generate an error message: negative x order 9 9 <HOST>|E|psPlaneTransformAlloc (FILE:LINENO) 10 Error: n1 is less than 0.10 Error: order1 is less than 0. 11 11 <HOST>|I|testPlaneTransformAlloc 12 12 Following should generate an error message: negative y order 13 13 <HOST>|E|psPlaneTransformAlloc (FILE:LINENO) 14 Error: n2 is less than 0.14 Error: order2 is less than 0. 15 15 16 16 ---> TESTPOINT PASSED (psCoord{psPlaneTransformAlloc()} | tst_psCoord.c) … … 25 25 Following should generate an error message for negative x order 26 26 <HOST>|E|psPlaneDistortAlloc (FILE:LINENO) 27 Error: n1 is less than 0.27 Error: order1 is less than 0. 28 28 <HOST>|I|testPlaneDistortAlloc 29 29 Following should generate an error message for negative y order 30 30 <HOST>|E|psPlaneDistortAlloc (FILE:LINENO) 31 Error: n2 is less than 0.31 Error: order2 is less than 0. 32 32 <HOST>|I|testPlaneDistortAlloc 33 33 Following should generate an error message for negative z order 34 34 <HOST>|E|psPlaneDistortAlloc (FILE:LINENO) 35 Error: n3 is less than 0.35 Error: order3 is less than 0. 36 36 <HOST>|I|testPlaneDistortAlloc 37 37 Following should generate an error message for negative w order 38 38 <HOST>|E|psPlaneDistortAlloc (FILE:LINENO) 39 Error: n4 is less than 0.39 Error: order4 is less than 0. 40 40 41 41 ---> TESTPOINT PASSED (psCoord{psPlaneDistortAlloc()} | tst_psCoord.c)
Note:
See TracChangeset
for help on using the changeset viewer.
