Need to compute sunspot frequencies ? Here's some code for you !

#! /usr/bin/python
# $Id: sunspots.py,v 1.1 2009/05/31 12:08:08 george Exp george $
#
# Plot sunspot data.
#
# Adapted from
#
#   http://linuxgazette.net/115/andreasen.html
#
# this was the original program.  Dataset from
#
#   http://linuxgazette.net/115/misc/andreasen/sunspots.dat
#
# Nice examples of using the newer plotting libraries can be found here:
#
#   http://www.daniweb.com/code/snippet691.html
#
# History:
#   Created: George Jones , 5/30/09
#
#   $Log: sunspots.py,v $
#   Revision 1.1  2009/05/31 12:08:08  george
#   Initial revision
#
#

import math
import pylab  # matplotlib
from scipy import *
import scipy.io.array_import

# Generate x,y datasets (year,wolfer)
tempdata = scipy.io.array_import.read_array('sunspots.dat')
year=tempdata[:,0]
wolfer=tempdata[:,1]

#
# Plot year vs number of sunspots
#

pylab.xlabel("Year")
pylab.ylabel("Wolfer number")

pylab.plot(year, wolfer, 'b')

# save the plot as a PNG image file (optional)
pylab.savefig('sunspots_time.png')

# show the pylab plot window
# you can zoom the graph, drag the graph, change the margins, save the graph
pylab.show()


#
# Take FFT of #s of sunspots, generating real and imaginary
#

Y=fft(wolfer)

pylab.xlabel("real(FFT)")
pylab.ylabel("img(FFT)")
pylab.title("Meas")
pylab.plot(Y.real, Y.imag, 'ro')
pylab.xlim(-4000,2000)
pylab.savefig('sunspots_FFT.png')
pylab.xlim(-4000,2000)
pylab.show()

#
# Compute the frequency of sunspots
#

n=len(Y)
power = abs(Y[1:(n/2)])**2
nyquist=1./2
freq=array(range(n/2))/(n/2.0)*nyquist

pylab.xlabel("Frequency [1/year]")
pylab.ylabel("|FFT|**2")
pylab.title("Meas")
pylab.plot(freq[1:len(freq)], power, 'b')
pylab.xlim(0,0.20)
pylab.savefig('sunspots_freq.png')
pylab.xlim(0,0.20)
pylab.show()

#
# Given the frequency, compute the period
#


period=1./freq

pylab.xlabel("Period [year]")
pylab.ylabel("|FFT|**2")
pylab.title("Meas")
pylab.plot(period[1:len(period)], power, 'b')
pylab.xlim(0,40)
pylab.savefig('sunspots_period.png')
pylab.xlim(0,40)
pylab.show()