Annales Geophysicae High-precision measurement of satellite velocity using the EISCAT radar

This paper presents a method of measuring the velocity of a hard target using radar pulses reflected from the target flying through the radar beam. The method has two stages. First, the Doppler shifts of the echo pulses are calculated at a high accuracy with an algorithm which largely improves the accuracy given by the Fourier transform. The algorithm also calculates the standard deviations of the Doppler frequencies with Monte Carlo simulation. The second step is to fit the results from a sequence of radar pulses to a velocity model allowing linear variation of the second time derivative of target range. The achieved accuracies are demonstrated using radio pulses reflected by a satellite passing through the beam of the EISCAT UHF radar working at 930-MHz frequency. At high SNR levels, the standard deviations of the frequency from a single pulse reach typically down to 0.2 Hz. The best standard deviations of velocity fit are below 5 mm s −1 while those of the second time derivative of range are below 1 cm s −2.


Introduction
Since the beginning of the space ages in the 1950s, the number of man-made objects flying around the Earth has been continuously increasing.In addition to satellites operating for different scientific and economical purposes, these objects also consist of space debris.The space debris consists of defunct satellites, explosion fragments of rocket stages and smaller objects resulting from collisions between larger ob-jects.The number of space debris is continuously increasing and becoming more and more harmful to space activity.
Various methods exist for tracking satellites and other orbiting objects.They include telescopes, TV-cameras, radars as well as lasers (Wakker et al., 1991;Goldstein et al., 1998;Isobe and Japan Spaceguard Association, 1998;Janches et al., 2000;Foster et al., 2005).Exceptionally precise orbit determination is obtained by means of satellite-to-satellite ranging and synchronised clocks at the satellites (Vonbun et al., 1978).In principle, the satellite path is determined by six orbital elements, but continued observations are needed, since the satellite orbits are variable.The reasons of variations are e.g., air drag and the fine structure of the terrestrial gravitational field.
The EISCAT UHF radar (Folkestad et al., 1983) was designed for observing the ionosphere by means of the incoherent scatter method.From the very beginning of the measurements in the 1980s, short-lived signals were observed which obviously were not due to incoherent scatter, but were interpreted as echoes from satellites or meteors passing through the main lobe or the sidelobes of the radar beam.It has later been understood that these echoes, which are disturbances in the incoherent scatter work, contain useful information on these targets.The EISCAT UHF radar has actually been used for meteor studies (Pellinen-Wannberg and Wannberg, 1994;Janches et al., 2002).
In addition to satellites and meteors, also space debris produces echo signals in radar receiver.Due to its high transmitting power, the EISCAT UHF radar is a great tool for observing space debris.For this purpose no specific space debris measurement is necessarily needed, but observations Published by Copernicus Publications on behalf of the European Geosciences Union.

T. Nygrén et al.: Satellite velocity
can be made in connection with the normal radar experiments planned for ionospheric studies.Such investigations have actually been made and a large number of space debris particles have been observed (Markkanen et al., 2002(Markkanen et al., , 2005(Markkanen et al., , 2009)).It has turned out that the radar is capable of observing particles with diameters down to 2 cm at a distance of 1000 km.
For tracking the paths of individual particles, their orbital elements should be known.A radar is able to measure the distance and Doppler shift of the target.The latter means that the line-of-sight velocity can be determined.If a tristatic radar is available and the target flies through the common scattering volume, all six orbital elements can be obtained.This is applicable to man-made objects like satellites and space debris as well as to natural objects like meteorite heads.The latter one is a scientific geophysical application, whereas the former one has mainly a practical use.Still, high-accuracy determination of satellite or space debris orbits may even have an importance in space weather studies, since the atmospheric density is affected by space weather.
The present paper introduces a method of determining the line-of-sight velocity of a radar target at a very high accuracy.A method of very precise determination of the satellite range will be presented in a later paper.When put together, the methods in these two works can be used for determining the orbital parameters of the target with an extremely high precision.The methods will be demonstrated using observations of satellite echoes made by the EISCAT UHF radar in November and December 2010.The power of the two methods is described by the standard deviations of the results; the best standard deviations of the velocity are of the order of a few millimetres per second and those of range are of the order of a few tens of centimetres.

Measurements
The data analysed in this paper is part of a larger dataset taken with the EISCAT UHF radar (Folkestad et al., 1983) during a measurement campaign conducted under the European SSA Preparatory Phase CO-VI activity in late November and early December 2010.The radar target in the present dataset is EUMETSAT's METOP-A satellite, which has a body size of about 6.3 m × 2.5 m × 2.5 m and large solar panels.The satellite was on a sun-synchronous polar orbit with an inclination of 98.7 • and an altitude of about 820 km.The radar transmitted pulses with lengths of 1920 µs at intervals of 20 ms.When the beam is directed to a point on the satellite orbit, several hundreds of echo pulses are received during a time interval of several seconds while the satellite passes through the radar beam.The signal is downconverted to the baseband and sampled at intervals of 1 µs.Data flow is continuous and it also includes the downconverted and sampled transmission signal.
The data flow is divided into IPPs (inter pulse periods) so that each IPP starts with the transmission pulse and, af- ter some delay, an echo from the satellite is observed.The transmission contains a phase modulation with a 32-bit alternating code (Lehtinen and Häggström, 1981).The downconversion shifts the transmission frequency to zero and, therefore, only the transmission envelope is visible within this part of the data flow.The echo pulse, on the other hand, lies at a non-zero frequency due to the Doppler shift caused by the satellite motion.
The complex modulation envelope of the transmitted pulse with a length T is (t) ∝ exp[iφ m (t)], where φ m (t) is the time variation of the modulation phase.Then the signal x(t) within the echo pulse is also proportional to (t), i.e., where ω is the angular Doppler frequency and φ is the phase.Thus, the echo signal can be demodulated by multiplying it with the complex conjugate of the modulation envelope.The resulting signal ] is a monochromatic signal multiplied by a boxcar function with a length of the transmitted pulse.The task is to determine the Doppler frequency of this pulse, which gives the line-of-sight target velocity.
Since the target velocity is not constant during the time when the pulse is reflected from the target, also the frequency of the echo pulse is not exactly constant.The consequences of this fact will be discussed later.
Figure 1 demonstrates the application of demodulation to a single IPP.The top panel shows the real part of the complex modulation envelope, consisting of several 180 • phase shifts of the alternating code.The second panel shows the front end of the real part of the radar echo.The echo starts at sample number 11323 and only 680 points from the beginning of the echo are shown (the grey area in the top panel corresponds to this part of the echo pulse).The phase shifts are visible in the echo pattern, and they are indicated by vertical dashed lines in the figure.The demodulated echo is plotted in the bottom panel.This shows that the demodulation removes the phase shifts.The result is a sinusoidal variation, except for clear peaks at the points of the phase shifts.These errors can simply be removed by interpolation.The imaginary parts are not portrayed here, but they show a similar behaviour.
The top panel of Fig. 1 indicates that the amplitude of the modulation envelope increases throughout the transmitting pulse.The amplitude of the imaginary part would have a corresponding decrease.This behaviour mainly results from a small drift of the signal phase during the transmission.It is worth pointing out that this phase drift is actually also phase modulation and the applied demodulation removes it together with the phase modulation of the alternating code.Hence, even if the radar would transmit only simple pulses, they would contain some phase drifts, which could be removed by the present method.
A more detailed study of the amplitude behaviour indicates that the radar power does not remain constant during the transmission of a single pulse, but it is slightly reduced.This means that even the amplitude of the echo pulse slightly changes.It would be possible to make a power correction to the echo pulses.However, test were made which indicated that the effects of the power drift on the accuracy of the determined Doppler shift is so small that it could be neglected.Therefore, no power correction is made in the present paper.

Frequency determination
The data from a single satellite pass contains the transmission envelopes as well as the sampled echo pulses and their arrival times.The first step in data analysis is to demodulate the pulses as described in Sect. 2.
The second step is to determine the Doppler frequency ν of each pulse.Then the the line-of-sight target velocity, i.e., the time derivative of range, is given by where ν r is the radar frequency and c is the speed of light.
Here positive velocity means away from the radar.This equation also defines the required frequency accuracy for a given velocity accuracy.In the present case, ν r = 929.6MHz.Then a velocity accuracy of 1 m s −1 implies a frequency accuracy of 6.2 Hz, but in order to go down to 10 cm s −1 , 1 cm s −1 or 1 mm s −1 , respective accuracies of 0.62 Hz, 0.062 Hz or 0.0062 Hz would be needed.
The data vector from an echo pulse consists of 1920 points at 1-µs intervals.A plain discrete Fourier transform of this data vector leads to a frequency step of 1/(1920×10 −6 ) Hz ≈ 520 Hz.One can improve the frequency accuracy by padding the data vector with zeros.However, in order to achieve a frequency step of 6.2 Hz (1-m s −1 velocity accuracy), the length of the padded data vector should be about 1.6×10 5 points (to pass this limit, a 2 18 -point fast Fourier transform should be used).Decreasing the frequency step by three orders of magnitude, to reach the velocity accuracy of 1 mm s −1 , would mean a data vector of 1.6 × 10 8 points.In order to pass this accuracy, a fast Fourier transform of 2 28 points is needed.
A multitude of methods for determining the frequency of a monochromatic pulse have been developed (Pisarenko, 1973;Chan et al., 1981;Kay and Marple, 1981;McMahon and Barett, 1986;Quinn, 1994Quinn, , 1997;;MacLeod, 1998;Aboutanios and Mulgrew, 2005;Provencher, 2010;Candan, 2011).The method used in this paper is explained in Appendix A. It gives both the frequency and phase of the echo signal at a high accuracy.It first calculates fast Fourier transform (2 16point transform is used in this paper) and then finds the maximum point of the amplitude spectrum at an essentially higher resolution.The true accuracy, however, depends on the SNR of the echo signal.The effect of noise on the frequency accuracy is studied in Appendix B, which shows a linear dependence between the logarithm of frequency error and the logarithm of the SNR.
The frequency analysis itself does not give error limits.Therefore, the standard deviation of frequency is calculated for a 1920-µs pulse using a Monte Carlo simulation as explained in Appendix B. Due to the linear dependence, it is sufficient to do the simulation for two SNR values only.When the SNR value of the echo is known, the standard deviation of the determined frequency is obtained from the linear variation.The Monte Carlo simulation could be done separately for the frequency of each pulse.However, since the frequency variation is rather small during a single satellite pass, the SNR simulation is here made only for a single frequency for each satellite pass.

Velocity fit
The duration of a satellite passing through the radar beam is a few seconds.During this time, the line-of-sight velocity of the satellite varies and it turns out that neither the second time derivative of range can be taken as constant.A sequence of velocity measurements v i = v(t i ), i = 1, 2, . . ., n with standard deviations σ i , i = 1, 2, . . ., n is obtained from a single satellite pass.The result of each Doppler measurement corresponds to the time when the pulse meets the target.Therefore, due to the varying distance between the radar and the satellite, the separations of the times t i are not constant.
The exact timing of the observations is quite important, since both the first and second time derivatives of target range change quite rapidly.The times of transmitting the front end of a pulse and of arrival of the front end of the echo are observed.By putting the start of transmission to zero, the echo T. Nygrén et al.: Satellite velocity from the front end will be received at some instant t = t f and from the rear end at t = t r .Then the front and rear ends will be reflected from the target at times t f /2 and t r /2, respectively.Since the rear end is transmitted at the time t = T , where T is the length of the transmitted pulse, and the target is moving at a velocity v, the length of the reflected pulse is T r = t r − t f = cT /(c − v).Using a pulse length T = 1920 µs and v = 5000 m s −1 , this gives T r − T ≈ 0.03 µs.Therefore, the change of pulse length is so small that it should not affect the frequency determination and there is no need to take it into account.
The next question is timing of an individual Doppler observation, when the pulse length during the impact can be taken as constant.An order-of-magnitude value of the second time derivative of range (see Sect. 5) could be 25 m s −2 , for instance.Then the variation of velocity within a single pulse of 1920 µs is 4.8 cm s −1 .The results in Sect. 5 show that this is clearly larger than the smallest standard deviation of velocity.Therefore, the accuracy of timing should be better than the pulse length.Each pulse is actually a chirp rather than being monochromatic.The Fourier transform of a chirp pulse is a function containing Dawson's integrals and, therefore, it is too complicated to be studied analytically.Hence, the problem must be investigated numerically.One would expect that the present method, which assumes a monochromatic pulse, gives the frequency at the centre of a chirp pulse.This is indeed what happens with chirp rates expected from satellite echoes.(If the chirp rate is too high, the resulting spectrum may be double-peaked and then the method does not work.)Hence, the correct timing is made by putting each t i to correspond to the time when the centre of the pulse hits the target.
The observed line-of-sight velocity is the first time derivative of range.This is given in the frame of reference fixed to the radar.One should notice that, although the second time derivative of range is the time derivative of this velocity, it is not the line-of-sight acceleration of the target, but only an acceleration describing the temporal variation of range.In what follows, it is still called acceleration for simplicity.
Assuming that the acceleration has a linear variation during the small time interval of the satellite pass, the values of the acceleration at the times of observation can be written as where a 0 and ȧ are constants.Then the velocities are The velocities can be collected to a single n-dimensional column vector where T indicates transpose.Correspondingly, the standard deviations of the velocities make a single vector The unknown parameters can be collected to an unknown vector Then the relation between the unknowns and measurements is where is a column vector containing the true errors of the components of V .The linear relation between the unknowns and the velocities is given by a matrix where 1 is a column vector with each component equal to 1 and the time vector T is Equation ( 7) presents a linear inversion problem which can be solved in a standard manner (see e.g., Nygrén et al., 1997).
The solution is where is the diagonal covariance matrix of the velocity vector V .
The error covariance matrix of the unknown vector is In conclusion, the most probable values of the parameters v 0 , a 0 and ȧ are the three components of x and the diagonal components of x are their variances.These results can be used to calculate the most probable value of the velocity at each instant of time.The fitted velocity vector is and its covariance matrix is The standard deviations of the velocities are given by the square roots of the diagonal elements of this matrix.
Knowing the best values and covariances of a 0 and ȧ, the best fit of a i can be calculated.By defining a vector xa = (a 0 , ȧ) T and a matrix the fitted values of a i are, according to Eq. ( 2), given by a vector Their covariance matrix is where the 2 × 2 matrix xa is the part of x corresponding to the components xa of x.

Sample results
Figure 2 shows an example of the SNR and the calculated Doppler shift from a single satellite pass.Here the horizontal axis is time from the start of the measurement.The top panel shows clearly the main and first sidelobes of the antenna pattern.The maximum SNR reaches to a maximum value of about 640.
The Doppler shift in the bottom panel of Fig. 2 is calculated with the method in Appendix A, starting with a fast Fourier transform of a length of 2 16  = 65 536 points.The positive Doppler shift indicates velocity towards the radar and the velocity is decreasing with time.The result reveals large inaccuracies between the lobes of the antenna pattern, where the SNR goes down to 10 −2 .The results to be shown later are always calculated from pulses with SNR > 2. These parts are plotted in red in Fig. 2.
Examples of frequencies determined from four beam passes are shown in Fig. 3.Here again, zero time refers to the the start of the measurement, and the analysis is limited to values SNR > 2. The frequency variations are very small and they take place around a trend similar to that in the bottom panel of Fig. 2. Thus, the trend is first calculated by making a linear weighted fit ν L = a • t + b to the determined Doppler frequencies.The results after subtracting the trend are plotted by dots in Fig. 3, where the continuous lines indicate 3σ limits.The values of the fitting parameters a and b are shown in each panel.
The results show that the statistical errors in the centre parts of the plots (regions of high SNR) are clearly below 1 Hz.Still, the four panels show quite different deviations from the linear trend.The top panel indicates a smooth systematic behaviour, while the second panel contains a variation, which looks sinusoidal with a period of the order of 2 s.The third panel contains a faster oscillation with a period of about 0.2 s and, finally, the bottom panel displays a more irregular variation.Considering the indicated 3σ limits, these variations cannot be due to statistics, but they must be real.A plausible explanation is that they are caused by satellite rotation.During different passes, the satellite is seen from different directions and, therefore, also the Doppler shifts from the rotating satellite may vary.This must depend on the geometrical structure of the satellite.However, it is also possible that some of the oscillations or irregularities are caused by scintillation due to instabilities in the ionospheric plasma.
When the Doppler velocities and their standard deviations from individual pulses are known, the velocities and their standard deviations can be calculated using Eq.(1).These are shown (after detrending) in limits by the thin lines.The standard deviations of velocities from individual pulses (not shown) lie between 2.5 cm s −1 at the centres and 50 cm s −1 at the edges of the plots.The behaviour of the fitting results is mainly determined by observations with small standard deviations and, therefore, the curves at the edges do not always follow the cloud of dots near the edges of the panels.
Fitting of the model drastically reduces the error of the velocity estimate.The standard deviations obtained from the four beam passes in Fig. 4 are plotted in Fig. 5.The results indicate that the smallest standard deviations from the middle part of each satellite pass is always smaller than 5 mm s −1 (the smallest minimum value is 3.4 mm s −1 on curve 1).Almost tenfold standard deviations are encountered on the edges.
The linearly varying accelerations a i with standard deviations (notice that this acceleration is the second derivative of range, not the target acceleration) are calculated in the four cases following Eqs.( 16) and ( 17).The results are plotted in Fig. 6.The results indicate that, indeed, the variation of acceleration is important within the interval of a few seconds Table 1.Results from 12 beam passes on 1 December 2010, after 16:00 UT.The beam passes are numbered from 1 to 12.The next two columns show the beam azimuths and elevations, t min and t s give the start times of the IPPs in minutes and seconds, t gives the times after the start of the IPPs, and the last four columns contain the line-of-sight velocities and accelerations with their standard deviations.The examples in Fig. 4 are from beam passes 1, 12, 10 and 8.  during the satellite pass.The top and bottom panels indicate increasing and the two centre panels decreasing acceleration (remembering that positive velocity is away from the radar).

No. az [
Figure 7 shows the time variations of the standard deviations of accelerations a i during the four beam passes.The smallest standard deviations in the centre of all passes is always smaller than 5 mm s −2 and reaches 15-23 mm s −2 at the edges.The smallest minimum is 3.7 mm s −2 and it is found on curve 1.
Since the results give the minimum points of the two standard deviations, it is possible to calculate the most accurate values of velocity and acceleration.They are obtained from Eqs. ( 13) and ( 16) using times of the minimum standard deviations.However, the minima of the two standard deviations do not occur quite at the same time.Since giving only simultaneous velocity and acceleration is feasible, one has to choose between giving either the most accurate velocity or the most accurate acceleration.In what follows, both velocity and acceleration are given at times when the standard deviation of velocity gets its minimum.
Table 1 contains results from a single satellite path on 1 December 2010 after 16:00 UT, obtained by this method.The radar beam is pointed sequentially at 12 different positions along the satellite trajectory waiting for the satellite to pass the beam.The beam directions are defined by the azimuth and elevation angles in columns 2 and 3 (the real pointing accuracy of the beam may be weaker than indicated and most likely the satellite does not pass through the beam centre).Positive and negative azimuth values mean eastwards and westwards from geographic north, respectively.Columns 4 and 5 contain start times (t min gives minutes and t s gives seconds) of the IPPs giving the smallest standard deviations.In column 6, t gives the times of the measured smallest standard deviations of velocity with respect to the start times.These should be added to the start times to get the times of the results.The correct timing also depends on the synchronisation of the EISCAT clock with UTC and on possible delays in the radar electronics, but these are not considered here.The next two columns contain the velocities and their standard deviations and the last two columns the accelerations and their standard deviations.
The azimuth behaviour shows that the satellite is flying westwards at latitudes to the north of the radar position.The maximum elevation angle is below 30 • so that the trajectory lies far from the radar zenith.The velocities are first negative, in accordance with the trajectory and then they turn to positive.The maximum velocities at the two ends of the trajectory exceed 4 km s −1 .Except for a single case of low SNR, the standard deviations of the velocity are always T. Nygrén et al.: Satellite velocity below 1 cm s −1 .The acceleration is always positive and the values lie below 40 m s −2 with standard deviations which are usually only a few mm s −2 .Again, one should remember that the acceleration here is the second time derivative of range, not the acceleration of the target.

Discussion and conclusions
This paper introduces a method of measuring the beam aligned satellite velocity from radar observations.The method consists of two parts: determining the Doppler shift from individual pulses and fitting a physical model to the results.In order to achieve the best possible accuracy of velocity measurements, the Doppler shift must be determined at an accuracy of the order of a few mHz.This is achieved by a method explained in Appendix A. After obtaining Doppler shifts of a sequence of echo pulses, a physical velocity model is fitted to the results.Model fitting is made by means of stochastic inversion, which implies that standard deviations of individual Doppler shifts must be known.The standard deviations depend of the SNR of the radar pulses.Although the pulses themselves do not give the standard deviations directly, the standard deviation for any frequency, pulse length and SNR value are obtained using the Monte Carlo method.This is explained in Appendix B and it gives the tools needed in model fitting.
The model takes into account the fact that the time derivative of the target range, i.e., the line-of-sight velocity, changes while the target is passing through the radar beam.In addition, it also allows that the second time derivative has a linear time variation.Although the satellite passes the radar beam in a short time interval of the order of 5-10 s, the latter variation must be taken into account when aiming at the best possible accuracies.
The effect of model fitting is that all applied measurements affect the results at each time of observation.Therefore, even Doppler shifts with large standard deviations reduce the standard deviations at times of the most accurate Doppler shifts.The result is that the model fitting reduces the standard deviations of the most accurate velocities by an order of magnitude.
The achieved velocity accuracy is so high that a very precise timing of measurements becomes important.For instance, if we consider that the second derivative of range is 25 m s −2 , which is a typical value in the above results, the velocity change during a single pulse of 1920 µs is 4.8 cm s −1 .This is larger than the best standard deviations of velocity by an order of magnitude.Hence, timing must be better than the pulse length.A more precise timing was obtained when it was noticed that, within the Doppler chirp rates encountered in practice, the frequency analysis always gives the frequency at the centre of a chirp pulse.
It turns out that the results contain features which are not completely in agreement with the applied model of satellite motion.These variations are larger than the statistical standard deviations.In some cases they seem irregular or in other cases they contain regular structures with periods of a couple of seconds or fractions of a second.It is possible that these structures are caused by the rotation of the satellite which has large solar panels and periodic structures on its surface.When the satellite is observed at different times on its orbit, the radar probes it from different directions with respect to the rotation axis and, therefore, the time variation of the Doppler shift may vary from one satellite pass to another.A second possible reason which could cause the variations is that scintillation due to ionospheric irregularities might affect variations in the signal path which show up as changes in the Doppler shift.In both cases, the model applied in fitting is not exactly valid and this causes errors which are not of statistical origin.There is no means of getting rid of such errors and the best one can hope for is that they more or less cancel out in model fitting.
The radar experiment applied in the present paper was designed for investigating the capabilities of the EISCAT radars in probing the paths of satellites and space debris particles.The resulting statistical accuracy of satellite velocity is of the order of a few millimetres per second, and that of range (to be presented in a second paper) is a few tens of centimetres.It seems clear that these accuracies are insurmountable.Hence, these methods have a permanent practical value.Due to the lower SNR values, however, the velocity accuracies of small space debris particles are not expected to be as high as those shown in the present paper.
In addition to practical purposes, the methods are also expected to have importance in geophysical investigations.The most obvious topic is velocities and positions of meteor heads with improved accuracy.Since the atmospheric density at high altitudes is affected by space weather, satellite orbits must vary accordingly.This offers a view on studies of geophysical conditions on the variations of orbital parameters.Since f (0) = 1, the solution of Eq. ( A6) is obtained by solving numerically the equation 1 24 The correct root x can be easily identified and it gives the unknown frequency The signal phase can be calculated form Eq. (A4) by choosing the phase angle (ν m ).This gives φ = (ν m ) + π ν m n 0 − π ν 0 (n 0 + 2). (A17) Although this approach already gives quite accurate results, the accuracy can still be improved in a simple way using the following error correction method.One can create a pulse with the frequency and phase obtained by the above method.When the analysis is repeated with this pulse, slightly different results are obtained.Since the original frequency and phase are known in this case, the errors of the second analysis are also known.Furthermore, since the results of the first analysis are already quite close to the true values, these errors are also quite close to the errors of the first analysis.Hence, these errors can be used to correct the results of the first analysis.Numerical tests indicate that this correction may increase the accuracy, for example, by two or three orders of magnitude.
Figure A1 shows the errors from frequency and phase analysis of pulses with frequencies ranging from zero to 0.5, where the frequency unit is the reciprocal of the sampling interval.The pulse length is 1920 samples and the data sequences have been padded with zeros up to a length of 2 13  = 8192 points.In this case, the frequency error is of the order of 10 −11 , while the frequency step given by the Fourier transform is only 1/8192 ≈ 1.22 × 10 −4 .The phase error is of the order of 10 −7 rad.
Even more accurate results are obtained if longer Fourier transforms are used.If the length of the Fourier transform is 2 16  = 65 536 points, the frequency error is of the order of 10 −12 before error correction and goes down to 10 −16 (the computer eps) after correction.The phase error after correction is of the order of 10 −12 rad.These accuracies are achieved, although the frequency step of the Fourier transform is only 1/2 16 ≈ 1.5 × 10 −5 .

Appendix B
The effect of noise The pulses used for calculating the results in Appendix A were free of noise.It is obvious that the presence of noise affects the signal spectrum even at frequencies close to the main peak.This necessarily reduces the accuracy of frequency and phase determination.Since noise is always present, it is important to study the effect of noise to the accuracy of frequency determination.An easy way of investigating this is to use Monte Carlo simulation with noise added to the monochromatic signal.
As an example, a pulse of 1920 data samples with unit amplitude, a frequency of 0.01 and a phase π/3 is used here.Gaussian random noise is added to the pulse.If the standard deviation of a basic random Gaussian signal is unity, a noise signal leading to a signal-to-noise value SNR is obtained by multiplying the basic random signal by an amplitude The noisy pulse is padded with zeros up to a length of 2 13 points, and the frequency is determined using the method in Appendix A. This is repeated 10 000 times for a set of SNR values.The standard deviations of the resulting frequency distributions give the error estimates in each case.These are shown in Fig. B1 for SNR values extending from 1 to 900.The results indicate that the errors are always better than 5 × 10 −6 but, at SNR values greater than about 10, the errors with the chosen frequency vary from 10 −6 to 10 −7 .
In spite of the noise, the errors are always much better than 1/2 13 = 1.22 × 10 −4 , which is the frequency resolution provided by the Fourier transform.In order to obtain a frequency resolution of 10 −6 , for instance, a 10 6 -point Fourier transform would be needed.

Fig. 1 .
Fig. 1.Example of data from a single IPP.Top: Real part of the modulation envelope.Middle: Real part of the echo pulse (only 680 first samples are shown).Bottom: Same as middle panel after demodulation.Grey shading in the top panel indicates the part of transmission corresponding to the part of the echo pulse in the other two panels.The vertical dashed lines in the middle and bottom panel indicate the positions of the 180 • phase shifts caused by the modulation envelope.

Fig. 2 .
Fig. 2. Example of SNR (top) and Doppler shift (bottom) from a single satellite pass.Curves corresponding to SNR > 2 are plotted in red.

Fig. 3 .
Fig. 3. Detrended Doppler frequencies from four beam passes (dots), together with 3σ limits (continuous lines).Time is measured from the start of each measurement.The subtracted trend is ν L = a • t + b with the numbers given in each panel.

Fig. 4 .
Fig. 4. Detrended velocities calculated from the Doppler frequencies shown in Fig. 3 (open circles), together with velocities fitted to the model v i = v 0 + a 0 t i + ȧt 2 i /2 (heavy line).The 3σ limits of the fitted velocity are shown by the thin lines.The subtracted trend is v L = a • t + b with the numbers given in each panel.

Fig. 5 .Fig. 6 .
Fig. 5. Standard deviations of fitted velocities from the four beam passes.Numbers on the curves refer to the panels in Fig. 4.

Fig. 7 .
Fig. 7. Standard deviations of fitted accelerations from the four beam passes.Numbers on the curves refer to the panels in Fig. 4.
Fig. A1.Errors of frequency (top) and phase (bottom) determination for pulses with different frequencies.The pulse length is 1920 samples and the length of Fourier transform is 2 13 points.A single phase value of π/3 is used in all cases and the error correction method is applied.
Fig.B1.The errors (standard deviations) of frequency determination as a function of SNR.The pulse frequency is 0.01 and the pulse length is 1920 samples before padding with zeros up to a length of 2 13 points.