Choosing Sampling time and number of data points

The previous simulations were made with a sampling time of 1000 yrs.  This is not enough to correctly retrieve the frequency associated with the precession of \omega, which for S/2000S5 has a period of 676.2 yr (the rules of having at least 10 data points per period apply).  We need to choose a smaller sampling time.  On the other hand, a sampling time too small would include in the frequency power spectra frequencies too large, so reducing the resolution.   Another problem to consider is the aliasing of frequencies larger than the Nyquist frequency into the sampled frequencies.   The problem is to find a sampling time small enough to accurately sample the \omega precession frequency, but large enough to have a good resolution.

To solve this problem I first performed an integration with a very small sampling time, 0.1 yr, which corresponds to a Nyquist frequency of 6.48*106 "/yr.   Using 32768 data points, the resolution was of 395.5 "/yr.  Results of the FFT of e*exp(i\omega) are in the following figures.  The largest frequency was at  1464.8 "/yr, which corresponds to a period of 885 yrs.  This is larger than the observed period of 676 yrs, and the difference is due to the poor resolution of this FFT.   Definetely, 0.1 yr is too small a sampling time.  However, this simulations show us that most of the power spectrum is between +/-1.6*104 "/yr, with the third largest frequency at way out at 8.4*104 "/yr.  1.6*104 "/yr would be therefore an excellent Nyquist frequency.   This corresponds to a sampling period of 40 yrs.
 
 

We therefore performed another simulation with this sampling time.   Using again 32768 data points, the resolution of the FFT would be of ~1"/yr.   Results of the FFT and the time evolution of e and \omega are in the following figures; this time the highest frequency was of 1916 "/yr, in excellent agreement with the 676 yr period of \omega precession.
 
 

The integration lasted 4 Myr, so it was possible to FFT the osculating elements of S/2000S5 three times using 32768 data points.   If, as in Robutel and Laskar 2001, we define

(4)                 \sigma = 1-f(2)/f(1)

(where f(2) is the frequency computed on the second interval, and f(1)is the frequency computed on the first) as a way to measure chaotic diffusion, the value found for S\2000S5 was of 2*10-3.