# Calculating the period of the sunspot cycle

Two weeks ago, I used the

sunspot number data provided by the Solar Physics Group at NASA's Marshall Space Flight Center to demonstrate positioning plots in window. This week, I'd like to show how to calculate the period of the sunspot cycle.
If you haven't already done so, download the sunspot numbers file and place it in your IDL path. Read it with the astrolib

READCOL procedure:

file = **file_which**('spot_num.txt', /include)
**readcol**, file, year, month, sunspots

Next, transform the sunspot series to the frequency domain and compute magnitude and power spectra:

mspec = **abs**(**fft**(sunspots))
pspec = mspec^**2**

(Aside:

FFT: it's all you need.)
I'd like to display the power spectrum as a function of frequency. This requires a few statements to set up a frequency vector based on the time data from the sunspot numbers file:

sampling_interval = **1**/**12.0** ; years, from file
freq_nyquist = **0.5**/sampling_interval ; years^{-1}
n_sunspots = **n_elements**(sunspots)
freq = **findgen**(n_sunspots/**2**)/(n_sunspots/**2**-**1**)*freq_nyquist

The sampling interval is one month. This allows me to determine the

Nyquist frequency (the maximum resolvable frequency), which along with the number of sunspots in the file, gives me the information to define a frequency vector. Note that I'm defining only positive frequencies up to the Nyquist frequency; this is because I'm going to display a one-sided power spectrum. Fold the negative Fourier modes into the spectrum, omitting the zero (DC) mode:

freq_nodc = freq[**1**:n_sunspots/**2**-**1**]
pspec_nodc = **2***pspec[**1**:n_sunspots/**2**-**1**]

and display the result on logarithmic axes:

p = **plot**(freq_nodc, pspec_nodc, /xlog, /ylog, $
xtickunits='exponent', $ ; IDL 8.2.1
xtitle='Frequency ($yr^{-1}$)', $
ytitle='Spectral Density', $
title='Power Spectrum of Sunspot Numbers, 1749-2010')

Note that I can set the style of tick units on the axes: 'exponent' and 'scientific' are new in IDL 8.2.1.
The power spectrum has a peak near 0.1 yr^{-1}. Locate the peak frequency with

MAX and mark it on the plot with

POLYLINE:

pspec_nodc_peak = **max**(pspec_nodc, i_peak)
freq_nodc_peak = freq_nodc[i_peak]
!null = **polyline**([**1.0**,**1.0**]*freq_nodc_peak, p.*yrange*, color='orange', /data)

The inverse of the peak frequency gives the period of the sunspot cycle: approximately 11 years. Mark this on the plot with

TEXT:

sunspot_peak_period = **1.0**/freq_nodc_peak
str = '$T_{peak}$ = ' + **strtrim**(sunspot_peak_period,**2**) + ' years'
!null = **text**(**1e-1**, **1e-4**, str, /data)

Here's my result:

Finally, just because it's interesting, test

Parseval's theorem (energy is conserved in the time- and frequency-domain representations of the signal) in several forms:

IDL> **print**, **total**(pspec)*n_sunspots, **total**(sunspots^**2**)
1.46437e+007 1.46437e+007
IDL> **print**, mspec[**0**], **mean**(sunspots)
52.0133 52.0134
IDL> **print**, **total**(pspec[**1**:*]), **stddev**(sunspots)^**2***(n_sunspots-**1**)/n_sunspots ; convert from sample variance
1965.67 1965.66

Conservation of greatness (of FFT)!