#!/usr/bin/env python3

# ================================================================
# 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 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 = [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), end=' ')
	print("%11.7f" % (f(x)), end=' ')
	y = ck[kmax+0]
	print("%11.7f" % (y.real), end=' ')
	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), end=' ')
	print()
	x += dx