#!/usr/bin/python

# ================================================================
# John Kerl
# kerl.john.r@gmail.com
# 2006-03-03
#
# This is an example of truncated Fourier series in one variable.
# ================================================================
# This software is released under the terms of the GNU GPL.
# Please see LICENSE.txt in the same directory as this file.
# ================================================================

# ----------------------------------------------------------------
from __future__ import division # 1/2 = 0.5, not 0.
from math import *
import sys

# ----------------------------------------------------------------
# Numerical quadrature for computing Fourier coefficients.

def fourcoeff(func, k, xlo, xhi, nx, dx):
	sum = 0.0
	x = xlo
	for i in range(0, nx):
		expnikx = cos(k*x) - 1j*sin(k*x)
		sum += expnikx * func(x) * dx
		x   += dx
	return sum / 2.0 / pi

# ----------------------------------------------------------------
# Function to be approximated via truncated Fourier series

def f(x):
	#return x
	#return pi - abs(x)
	#return cos(x)
	#return x^2
	#return pi**2 - x**2
	#return cos(1.7*x)
	#return cos(1.7*x+0.4)
	#return x**3
	#return cos(4*x) + cos(5*x) + cos(6*x) + 1.2 * cos(7*x)
	#return exp(x)
	#return 1 / (1 + x**2)
	#return cosh(x)
	#return sinh(x)
	#return exp(-pi*x**2)

	#a = pi/5
	a = pi/2
	if (abs(x) <= a):
		return 1.0
	else:
		return 0.0

# ================================================================
# Start of main program

xlo = -pi;
xhi =  pi;
nx  = 400;
dx  = (xhi - xlo) / (1.0 * nx);

kmax = 20;

# ----------------------------------------------------------------
# First, compute the Fourier coefficients.
# Index the ck array from 0 to 2*kmax+1.

ck = range(0, 2*kmax+1)
for k in range (-kmax, kmax+1):
	ck[k+kmax] = fourcoeff(f, k, xlo, xhi, nx, dx)

if 1:
	print "Fourier coefficients:"
	for k in range (-kmax, kmax+1):
		print "%3d: %11.7f %11.7f" % (k, ck[k+kmax].real, ck[k+kmax].imag)
	print

# ----------------------------------------------------------------
# Next, tabulate the approximations.

x = xlo
for i in range (0, nx):
	print "%11.7f" % (x),
	print "%11.7f" % (f(x)),
	y = ck[kmax+0]
	print "%11.7f" % (y.real),
	for k in range (1, kmax):

		cposk = ck[kmax+k]
		expikx = cos(k*x) + 1j*sin(k*x)
		y += cposk * expikx

		cnegk = ck[kmax-k]
		expnikx = cos(k*x) - 1j*sin(k*x)
		y += cnegk * expnikx

		print "%11.7f" % (y.real),
	print
	x += dx