Determining the orientation of the rotation axis of a flashing satellite

Doing what is described in the title has been an old dream for several years now. Knowing the direction of the rotation axis would allow a more exact determination of the rotation period (instead of the flash period we usually measure) and a better comparison of the characteristic times for exponential growth with theory. Following the evolution of the direction of the rotation axis is just as interesting as following the evolution of the rotation period, if not more. It would give insight in which forces are working on the rocket (interesting in case accelerations happen).

Mike McCants and I had a rather intensive discussion on Chapter 6 of the new brochure in English, which deals with the synodic effect. The end-result of the discussion was that Mike tracked 83- 69 J on January 1, from 02h00m59s till 02h10m40, i.e. from an elevation of 36 degrees down to 3.5 (!) degrees in the north, in order to catch as much as possible of the synodic effect. Here is what Mike has to say about the observation :

I have chosen period = 9.160 in order to fit the first 18 or 19 observations.
This gives about a 2 second error over the last 30 cycles or so, indicating
a change in period of about 0.07 seconds would be required to fit that.
That seems to me to be approximately the right size for a synodic effect.

 count         .00    1.00    2.00    3.00    6.00    7.00    9.00   11.00
 obs           .00    9.09   18.41   27.28   55.17   64.04   82.48  100.39
 pred          .00    9.16   18.32   27.48   54.96   64.12   82.44  100.76
 diff          .00    -.07     .09    -.20     .21    -.08     .04    -.37

 count       12.00   13.00   14.00   15.00   19.00   20.00   21.00   22.00
 obs        110.24  118.85  128.55  136.73  173.33  183.22  192.28  202.14
 pred       109.92  119.08  128.24  137.40  174.04  183.20  192.36  201.52
 diff          .32    -.23     .31    -.67    -.71     .02    -.08     .62

 count       23.00   28.00   29.00   32.00   33.00   34.00   35.00   37.00
 obs        210.41  256.20  265.89  293.74  303.59  312.44  322.00  339.88
 pred       210.68  256.48  265.64  293.12  302.28  311.44  320.60  338.92
 diff         -.27    -.28     .25     .62    1.31    1.00    1.40     .96

 count       39.00   41.00   42.00   43.00   44.00   45.00   47.00   48.00
 obs        357.88  376.89  386.74  394.48  404.80  413.93  432.19  441.34
 pred       357.24  375.56  384.72  393.88  403.04  412.20  430.52  439.68
 diff          .64    1.33    2.02     .60    1.76    1.73    1.67    1.66

 count       49.00   50.00   51.00   53.00   54.00   55.00   56.00   58.00
 obs        450.41  460.45  468.91  487.26  495.72  505.42  514.92  533.46
 pred       448.84  458.00  467.16  485.48  494.64  503.80  512.96  531.28
 diff         1.57    2.45    1.75    1.78    1.08    1.62    1.96    2.18

 count       59.00   60.00   61.00   62.00   63.00
 obs        542.30  551.65  561.34  570.69  580.23
 pred       540.44  549.60  558.76  567.92  577.08
 diff         1.86    2.05    2.58    2.77    3.15

After some serious handwaving, I decided I needed to write a computer program that calculates the direction of the rotation axis and hence also gives the rotation period. Since I didn't want to spend too much time on the problem, I used some of Mike's prediction software (thanks Mike!) and an easy approach to the reflection problem. I finished the software in 2 days and debugged it for another week...

Patrick Wils gave a good introduction of the reflection problem in his article in Flash 89 (last month). The theory behind my approach is very straightforward and allows a quick and 'dirty' calculation, which was what I was looking for.

I assume the timings to be very correct (no residuals) and precise. At each of the timings I calculate the direction of the sun and that of the observer, both viewed from the satellite. I know that a flash will occur whenever the reflecting surface of the rocket (a cylinder's body) is perpendicular to the bisectrix between the vector satellite-sun and the vector satellite-observer. I assume that the rocket is tumbling end-over-end, which makes that the rotation axis is always perpendicular to the rocket body's long axis. So, I have two vectors (rot. axis and bisectrix) that have to be perpendicular to the body axis. If I take the vector product of the rot. axis and the bisectrix I get the body axis at the time of the flash.

Obviously, I don't know the direction of the rotation axis. So, I try out all possible directions and do the above calculation for each of these directions. Assume I take a certain direction of the rotation axis (right ascension and declination ). For each of the 45 timings that Mike made, I calculate the direction of the body axis of the rocket. I know the time difference between two flashes, and I also know the direction of the body axis at both these times. I calculate the angle between the two directions of the body axis and I know that . Here is the time difference between the two flashes. From this formula, I can calculate . I now proceed to calculate for each couple of successive flashes. So, for each couple of and , I end up with 44 values of .

We know that in reality the rocket has one single rotation period that remains as good as constant during one pass. So, we need to find that couple of and for which the 44 values of are one and the same. That would then be the direction of the rotation axis and the one value is the rotation period.

Naturally, our measurements are never infinitely accurate. The observer suffers from a variable reaction time. So, we can't expect all measurements to be spot on the flash. That's why we can't expect to find a direction of the rotation axis for which all 44 values of coincide. What we can expect is to find directions for which the values of are grouped more tightly around their average. One measure for this 'tightness' is the variance or the standard deviation . The variance is the mean squared-deviation of each value for with respect to the average rotation period. Wherever is minimal, we can hope to find the direction of the rotation axis.

What does this look like in practice? Figure 1 shows a twodimensional contour plot of for Mike's 83-69 J timings. On the x-axis we find the right ascension in degrees (from 0 to 360), the y-axis shows the declination in degrees (from -90 to 90) of the rotation axis. The lines connect areas where the standard deviation is the same. The step between each contour line is 0.005 seconds, only areas where the standard deviation is lower than 0.47 seconds are shown. The white areas at 10 to 140 , 50 to 90 on the one hand, and 220 to 320, -80 to -50 on the other hand are areas where the standard deviation is substantially higher than in the rest of the sky. These areas can most probably be excluded as possible directions of the rotation axis.

Figure 1 : Contour plot of the standard deviation as a function of the right ascensions and declination , for lower than 0.47 seconds. Step size of consequent contour lines is 0.005 seconds. Satellite is 83- 69 J.

Figure 2 shows a similar contour-plot. This time I have excluded all regions in the sky where the standard deviation is higher than 0.42 seconds, with a step of 0.0025 seconds between two contour lines. This time only a few specific regions come into play. The two isolated regions on the right actually have the lowest values of , down to 0.36 seconds. So, it would seem these two isolated regions are the most probable directions of the rotation axis. But wait a minute, why only these two regions? Why not the whole huge region as shown on figure 2? It's no good limiting the regions artificially to a certain threshold value, we should try to prove the two isolated regions are significantly better choices for the direction of the rotation axis.

Figure 2 : Contour plot of the standard deviation as a function of the right ascensions and declination , for lower than 0.42 seconds. Step size of consequent contour lines is 0.0025 seconds. Satellite is 83- 69 J.

Let's look at one region in detail, the one around . Figure 3 is a three-dimensional plot. On the x-axis (values around 300) is the right ascension, on the y-axis the declination (values around -55) and finally on the z-axis we find . So the higher the curve gets, the better the fit is, the more chance the direction in question really was the rotation axis. You can clearly see a bump upwards at around 305, -50. Its radius is a few degrees wide. The average rotation period is 18.27. The flash periods Mike observed were of the order 9.16 to 9.19, which makes the rotation period not all too far away from that (9.14). But the rotation period isn't really all that accurately known from my analysis, since there is quite some dispersion (given by the standard deviation of about 0.36 seconds) in the values in that region. For a more accurate determination of I would need to do a more rigourous analysis where each flash is calculated starting from a certain value of . This is exactly the kind of program I tried to avoid, since it is more complex (involving iterations -see Sat 4.75) and its results are much more difficult to interprete. You would need to do a multiple of the calculations I need to do. You need to fix the four unknowns -see Patrick Wils' article last month- whereas I only have two unknowns (alpha and delta)! This, in fact, is the elegance of this method.

Figure 3 : as a function of right ascension and declination . The bump upwards is a region on the sky which represents a better fit of the data. Satellite is 83- 69 J.

But I still need to prove that our regions are statistically significant. This will be a bit more technical. Firstly it turns out that the 44 values of are for almost all directions of the rot. axis normally distributed around the average (using a test). Only in the regions that are blank in figure 1 (they are in fact the regions where the rot. axis is almost parallel to the resp. bisectrix) are there significant deviations from the normal distribution. Most of the sky has the same (about 0.42). The two regions I think might be the rotation axis have values of about 0.37. If I now perform an F-test to see whether these two variances are significantly different (taking into account that the average value is observed much more often in the parameter space), I find that this is indeed the case. The chance that the real variance is the same in the whole sky (despite the 0.37 -0.42 discrepancy) is less than 1%. This makes my statement that those two regions are better fits statistically significant.

In fact the statistical calculations in the above have been made with 33 observations, not all 45. I excluded the observations that were obviously suffering from large residuals (larger than 0.8 seconds). The interesting thing is that the two regions which I found to have significantly different variances, were already (though less pronounced) present in the case for where all 45 observations were included. The two regions weren't statistically significant though. It is quite normal that the observations with high residuals influence the significance so clearly, since the calculation of the variance gives more weight to observations with high residuals. The 12 observations that were excluded were spread quite uniformly over the whole pass.

Some interesting tidbits that may help show that our calculations are not unrealistic are shown in figure 4 and 5.

Figure 4 : Histogram of rotation period as determined for each couple and . The x-axis has the rotation period (or rather half of it), the y-axis the relative occurence frequency of that rotation period. Satellite is 83- 69 J.

Figure 4 shows a histogram of the half periods obtained with my method for each couple of and . The double peak is intriguing. One is higher than the flash period, the other is a bit less than the observed flash period. The exact position of the peaks on the x-axis could well be off by several hundredths of a second. Figure 5 shows which regions in the sky produce a higher rotation period, and which a lower period. Ultimately it is the sense of rotation that is behind these graphs. The synodic effect will be different depending on the sense of rotation.

Figure 5 : Rotation Period as a function of right ascension and declination . Periods range between 18.4 and 17.8 seconds. Satellite is 83- 69 J.

Finally, I performed this method of analysis to the observations of Walter Nissen of 91-29 B, on which Patrick Wils already reported last month. Figure 6 shows the right ascension on x, declination on y and . This time the peak is better defined than in the example above. This is due to the extreme geometry Patrick alluded to. Our solution is found at 16h20m RA and 31.9 declination. Patrick found a solution of 16h28m and 29.9 declination. The difference in declination may be due to the lower accuracy of the Sat calculations, or alternatively due to the way I define the 'rotation axis'.

Figure 6 : as a function of right ascension and declination . The bump upwards is a region on the sky which represents a better fit of the data. Satellite is 91- 29 B.

It is a common adage that one can prove everything (false or true) using statistics. I'm pretty sure that the rotation axis of 91- 29 B can be reliably determined with my method. I do feel however that the 'significance' of the determined rotation axis for 83- 69 J may not mean all that much. It is rather suprising that the best-of-fit regions are very close to worst-of-fit regions, as was also the case for 91- 29 B. However, Mike did not see anything really out of the ordinary, unlike Walter who definitely saw a very weird pass. Mike just saw a gradual synodic effect, not an extreme one! So, it is possible that the determined directions of the rotation axis are just a spurious result. I am however very confident that the rotation axis is to be found somewhere in the vast region as shown in figure 2. This article is partly meant to provoke reactions, comments and opinions, since I myself am not sure how accurately the rotation axis can be determined with my method.

One other thing I am pretty confident about is that the method will be able to 'pinpoint' a rotation axis given simultaneous observations of a satellite. Though we have not been able to perform simultaneous observations as yet (for one because Mike lives in Austin, Texas and I live in Munich...), we (coincidentially) do have three observations of 82- 40 J during 24 hours. Mike made 31 timings on January 28, 2h21m UT, I made 20 timings on January 28, 18h07m UT and finally Mike observed 39 periods on January 29, 29, 1h52m UT. It is probably a safe assumption that the rotation axis did not change its direction a lot during those 24 hours. We can reasonably argue that the usual torques (graviational and magnetic) are too small to change the direction over one day. Furthermore there were no signs of accelerations (in rotation period), during which higher torques (due to fuel leaks) might happen.

So, we can try to use all three sets of measurements to find the best fit for all of them. In this way we can probably circumvent the possibly spurious solutions (i.e. the ones close to where the bisectrix was pointing to) since for most of these measurements those 'bisectrix' regions are different. Figure 7 shows the regions where the bisectrix between the vectors sat-sun and sat-obs are. There are six regions in total for three strings of measurements. Each string of measurements has two bisectrix regions in opposite directions. The two longer strings are Mike's observations. They nearly coincide because he observed the satellite from the same place and during similar passes (e.g. both are in the north, fictitious example). So for future determinations it may be preferrable to observe passes that have different sat-obs vectors. The bisectrix region of my observation is quite distinct from the those that Mike did, simply because I observed it from a different place under different geometrical conditions. The sun's position of course did not change, which is why all of these bisectrix regions are still about in the same region of the sky.

Figure 7 : X-axis is right ascension , y-axis is declination . Shown here is the location in the sky of the bisectrix between the vectors satellite-sun and satellite-observer, for each of the three measurements used in this analsysis of 82- 40 J data.

If we make contour plots of the standard deviation (as we did for 83- 69 J) for each of these measurements, we are again not sure which solutions are 'real' and which are 'spurious', i.e. in the vicinity of the bisectrix regions. If we now plot the sum of each we get figure 8 We have three of those values in each point of the sky, as determined for both of Mike's and my observations.

Figure 8 : The sum of as a function of right ascension and declination (or rather ). To improve visibility a skirt is drawn around this surface plot. It can be clearly seen that the fit gets better closer to the north pole region. The large dropouts visible below the x and y-axis are due to the perspective.

Again the x-axis is and the y-axis is , while the z-axis is the sum as described above. We now have a measure of how good the fit is for the three measurements combined. Where z is maximal we find the best fit. It turns out that this is the region around , i.e. the north pole. Whereas the solutions around and are the best fits as determined from the individual contour plots, they are not very good fits if we combine the three measurements! So, as I had hoped, combining quasi-simultaneous observations does get rid of those solutions that are adjacent to bisectrix regions. The one region that is left is the north polar region.

Figure 9 : Contour plot of the sum of the standard deviations (as determined for each of Mike's and my measurements) as a function of the right ascension and declination , for lower than 4.0 seconds. Step size of consequent contour lines is 0.025 seconds. Satellite is 82- 40 J.

Figure 9 shows the contour plot of the sum of . Only values below 4.0 seconds (i.e. the average standard deviation is thus 1.33 seconds) are plotted. Again the polar region sticks out. One could argue that the polar region is also adjacent to the bisectrix regions around and , but the angular distance is over 30 degrees here, while for the other regions (around ) that angular distance is more something like 5 to 10 degrees.

To conclude this long article, it looks like a coordinated effort to observe a rocket quasi-simultaneously under as diverse geometrical conditions as possible, might well lead to a determination of the rotation axis of that rocket. The software to analyze observations of this kind has been written, we should now come up with the observations. It would probably be best to follow one (or several) rockets over a longer time period, so the evolution of the rotation axis direction can be followed. I would love to hear your comments and reactions!

Bart De Pontieu

Many thanks to Mike McCants for supplying large pieces of the prediction software.