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:
and:
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.
