Patrick Wallace, Nicole Capitaine, and Dennis McCarthy,
Thank you very much for your reply in the forum. I apologize that this has taken me so long to reply, but I wanted to get most the questions ironed out. I agree with your position to use the new techniques and I’ve been working to make that happen, at least within AGI. Your answers helped greatly in achieving working code for the several different options, and the example problem helped get those last few bugs out. I have corresponded with some other folks in the field, and in particular, Scott Polorelc has helped with some additional insight and test cases. It turns out that the example problem I have for my text was pretty close, but I’ve made one or two minor corrections to be more in line with results I’m now obtaining (this is the LEO case I show later). I have looked at the example problem and matched nearly all the matrices, for both approaches, etc. I’m still compiling results on that, but within the accuracy of MatLab, I match the test case. The paper you suggested last time (AA 412) was also helpful.
1. As a process to refine my codes, I first matched CIP coordinate values to sample values I had received from Scott. I refer to 3 approaches - the complete iau00xys with the long series summations, and the iau00pna and iau00pnb methods. These are shown in the attached pdf (ch03opt.pdf) as Alg 12, 14, and 15, and the last two use the iau2000a and iau2000b respectively. Although the TN32 (latest download from 13 Jul 04) describes how to obtain the CIP coordinates from the dpsi and deps values (eq 18-24) (Alg 13 in ch03opt.pdf), the discussion of approaches to form the transformation matrices (Section 5.10 and method 1, 2a and 2b) does not seem to use this approach. I have tried to summarize the pertinent equations from these methods in the pdf attachment, and in my Matlab scripts. The SOFA and IERS routines seem to follow this same convention. You can see the approximation of the IAU2000b theory in the comparisons. All units are arcsecs except for the first case which is the IERS test case provided in the forum. The second set of values are in radians.
> testxys
======
00xys x 68.058443366770 y 5.636327037859 s -0.002900374048 a 0.000138888893
x 0.00032995664459 y 0.00002732568459
00pna x 68.058443081208 y 5.636326168085 s -0.002900373900 a 0.000138888893
x 0.00032995664321 y 0.00002732568037
00pnb x 68.058064390369 y 5.636322053613 s -0.002900716853 a 0.000138888893
x 0.00032995480726 y 0.00002732566043
======
00xys x 206.610841582264 y 2.585385622140 s 0.001595200421 a 0.000138888924
00pna x 206.610841858982 y 2.585384376360 s 0.001595201043 a 0.000138888924
00pnb x 206.610905554775 y 2.585129796599 s 0.001594803986 a 0.000138888924
======
00xys x 393.886209774124 y -2.794390800668 s 0.000371446625 a 0.000138889016
00pna x 393.886210272449 y -2.794391474989 s 0.000371447273 a 0.000138889016
00pnb x 393.886586470054 y -2.794170217647 s 0.000370788614 a 0.000138889016
======
00xys x 607.721390016045 y -3.305193180132 s 0.006690153552 a 0.000138889190
00pna x 607.721390406484 y -3.305191752409 s 0.006690151452 a 0.000138889190
00pnb x 607.721312850785 y -3.305153289299 s 0.006689731104 a 0.000138889190
======
00xys x 794.647293098830 y -1.071208974025 s -0.003411330324 a 0.000138889404
00pna x 794.647294065373 y -1.071209230503 s -0.003411329828 a 0.000138889404
00pnb x 794.647083199480 y -1.071281449283 s -0.003411456217 a 0.000138889404
======
00xys x 1006.970440352566 y -10.953839627805 s 0.021699015926 a 0.000138889717
00pna x 1006.970440551235 y -10.953840238664 s 0.021699017423 a 0.000138889717
00pnb x 1006.970439779244 y -10.954314442631 s 0.021699986129 a 0.000138889717
======
00xys x 1196.602882840200 y -2.167484619281 s -0.008771623367 a 0.000138890057
00pna x 1196.602883092569 y -2.167485947154 s -0.008771619514 a 0.000138890057
00pnb x 1196.603070252649 y -2.167596684449 s -0.008771367689 a 0.000138890057
======
00xys x 1405.105826206628 y -19.523092453299 s 0.045532389922 a 0.000138890500
00pna x 1405.105827953005 y -19.523092233876 s 0.045532389257 a 0.000138890500
00pnb x 1405.105863827744 y -19.523109446006 s 0.045532395894 a 0.000138890500
======
00xys x 1599.804510166085 y -6.213280729265 s -0.010652329042 a 0.000138890978
00pna x 1599.804510268629 y -6.213280055773 s -0.010652331652 a 0.000138890978
00pnb x 1599.804833695014 y -6.213619795053 s -0.010650927002 a 0.000138890978
======
00xys x 1802.300060269284 y -27.887115008647 s 0.072719272201 a 0.000138891541
00pna x 1802.300062141819 y -27.887115252719 s 0.072719273394 a 0.000138891541
00pnb x 1802.299642385323 y -27.887228395478 s 0.072719747298 a 0.000138891541
======
00xys x 2002.931741313923 y -13.781836749762 s -0.001407445031 a 0.000138892163
00pna x 2002.931740667369 y -13.781837074850 s -0.001407443474 a 0.000138892163
00pnb x 2002.932019636299 y -13.782315154329 s -0.001404970604 a 0.000138892163
======
00xys x 80.531879792773 y 7.273921786482 s -0.003026567667 a 0.000138888894
00pna x 80.531879525285 y 7.273920773562 s -0.003026567465 a 0.000138888894
00pnb x 80.531879241826 y 7.273838229790 s -0.003026716378 a 0.000138888894
>
2. The next step of testing was to examine some test cases. For these, I have used the iau00xys approach shown above. I selected two test cases, a LEO and a GEO satellite. The attached following spreadsheet values show the results. I setup the original excel file so that one could place any position and velocity value in the “baseline” line, and the differences would be calculated for all the remaining values. Currently, the iers ceo-based method is set as the baseline. I may not have the “sofa Eq based” inputs set correct (dpsi, deps).
3. I have generated a few more questions, but at least I believe now that I am running something that is close to the desired technical results.
a. The dspi frame bias appears as two different values in the TN 32 pg 40 eq 17 (0.016617”), and in the IERS CPBN2000 routine (0.041775). These values show up in the AA 412 paper as well. The example problem seemed to use the latter value and I’m confused on which to use.
b. You asked about my previous question concerning the plot of the dx dy values and I’ve included the graph from the excel data. The dx and dy values are from the bulletin b values of dx and dy. The 91 and 99 data look well outside the “norm”.
c. The TN13 shows the 1980IAU theory of nutation fundamental arguments with only up to t^3 values, while the TN21 and TN32 show up to t^4. I presume that one should use the later expanded expressions. Also, TN21 shows the second term for F5, or omega, as 6962890.2665”, while the others have it as 6962890.5431” (difference is the fractional portion). Again, I presume the later versions are correct for implementations using FK5 and those using the IAU2000, but I wanted to be sure.
d. In the attached pdf (ch03opt.pdf) file, I show the several intermediate coordinate systems with the old FK5 theory. In the new IAU2000 theory, there also appears to be a “pseudo-earth fixed” pef system, but there could also be an intermediate reference epoch (ire). I think it’s good to reduce the number of available options as it should promote less confusion. However, I wonder if the ire will be the same if one uses the CEO-based or equinox based approaches. Specifically, they both use the same polar motion matrix. However, since the CEO uses the ERA angle and the equinox based uses the GST angle, the resulting vector will be [slightly] different between the two approaches. Thus, the ire could potentially differ depending on which approach one uses.
e. The question I had asked last time about the “FK5” approach referred to the Alg 15 (ch03opt.pdf). Essentially, I wondered if the only change for IAU2000 was to add the delta delta psi and eps corrections from finals.data?
Thank you again for your consideration.
dav
David A Vallado
Technical Program Manager, CSSI/AGI
719-573-2600, 719-573-9079 FAX
E-Mail:
Web: http://www.centerforspace.com