We consider specular reflections off a cylindrical reflective body, tumbling end-over-end, i.e. the rotation axis is perpendicular to the main body axis (see Figure 1). Note that in this case the rotation axis coincides with the angular momentum vector and the instantaneous rotation vector (assuming external torques are zero or small). We take a reference frame whose origin is at the center of mass of the satellite. The x, y and z axes are the usual astronomical directions. This means that the x-axis points towards the vernal equinox point, the z-axis towards the celestial pole and the y-axis completes this right-handed frame. We call the unit vector satellite-sun , the unit vector satellite-observer and the unit vector along the rotation axis . Note that the rotation about the rotation axis is defined in a right-handed way, and that the angular coordinates of are () so that , and . The cylindrical body axis is assumed to be perpendicular to at all times (i.e. 'end-over-end tumbling'), and rotating about with the rotation period .
Figure 1: The reference frame used in our calculations.
From Snell's law we know that the observer will only see specular flashes off the rocket if the bisectrix between and is perpendicular to . Let us call the unit vector along the bisectrix between and the vector . So, the reflection condition becomes:
Whereas we also assume an end-over-end tumble:
From (1) and (2) it is straightforward to deduce that:
It should be noted that for a typical pass of an artificial satellite, the vector is time-dependent, since the satellite moves with respect to the observer. The sun's position in our coordinate system is constant for all practical purposes, given that a satellite pass lasts for 5 to 10 minutes and the sun's (apparent) celestial motion is about 1 degree per day. We will assume that both the rotation axis and the rotation period remain constant during the pass. Both of these assumptions are pretty safe as we will justify below.
Let us now consider the problem at hand, the determination of the rotation axis . The input parameters consist of an array of n times () at which the observer recorded a flash (visually) and also of an accompanying array of n indices , expressing the index number of the flash observed. The latter takes into account that not all flashes timed are consecutive. Usually though, the flash index array will be identical to an array containing the consecutive numbers . If the observer misses (say) the second flash, the index array would be .
Further input parameters are the observer's geographical coordinates and the satellite's USSPACECOM orbital elements with epoch around the time of observation. Using these input parameters, the satellite ephemeris software program (SGP4) can calculate the vector at each time accurately enough for the present analysis. The position of the sun in the satellite's rest frame remains the same during one satellite pass (which typically lasts 5 to 10 minutes) for all practical purposes. It can be calculated with popular astronomical software programs.
The remaining unknown parameters of the problem are and the rotation period . If these were known, it would be possible to calculate at each flash, using the time of that flash. From the it is possible to calculate the flash periods :
in which is the angle between and . The observed flash periods are given by .
A laborious and not very successful way to determine and consists of calculating for all possible directions of , and all possible values of . One then looks for the and values that minimize the differences . To improve the minimization, one can also introduce a fourth unknown, the exact time of the first flash. This approach is very laborious because it implies searching for four unknowns. It is not very successful because the accuracy of the visual observations is not high enough to reliably determine , the time of the first flash and the direction of at the same time.
Our new method avoids both problems by getting rid of the third unknown and the fourth unknown (the time of the first flash).
We treat as the only unknowns. We calculate and , and from those vectors also . Using , we can calculate for each couple by using (3). It is then straightforward to calculate . So far, we have described the approach described in the above, which we labelled as laborious and unsuccessful.
The crucial difference is that we do not treat as an unknown, but calculate it from (4) by postulating:
where is given by the observations.
From (4) we know that (5) defines , i.e. we can calculate for each and each flash period a corresponding rotation period. So, instead of treating as an unknown, we calculate for each pair of successive timings i and i+1 and for each direction of the rotation axis.
For any random with coordinates , these n-1 rotation period values are different from one another. They show a standard deviation around an average where:
Only for the real direction of the rotation axis are the n-1 rotation periods identical to one another. Only for that couple is the standard deviation equal to 0.
This then is how one can determine the rotation axis in a much less laborious fashion. One only has to scan for a zero value. One thus scans through two unknown parameters ( and ) instead of four. In our method, the third unknown () is a side-result from the search for the other two, whereas the fourth unknown doesn't play a role, since we assume that the timings made by the observer do not suffer from inaccurate measuring, i.e. we assume that the times observed are exactly the times at which the flashes occurred (so the time of the first flash is now defined too).
Of course, the above assumes that the timings made by the observer are perfect, i.e. that they have no residuals. This is not a very realistic assumption, given the human reaction time and other factors, which is also why in all practical cases on cannot hope for a specific direction for which all n-1 values are identical. One can however hope for a direction (,) for which becomes minimal.
Our method thus consists of looking for a direction of with coordinates for which is minimal. If the minimum is statistically significant, the direction of the rotation axis has been found, and an estimate for the rotation period can be calculated by taking . Naturally, the accuracy of will only be as good as the original accuracy of the timings. The power of our new method lies in the fact that the rotation axis can be found without accurate knowledge of . This is demonstrated in the next section.