Calculate the Fourier series with the trigonometry

I try to implement the Fourier series function according to the following formulas:

...where...

...and...

Here is my approach to the problem:

import numpy as np
import pylab as py

# Define "x" range.
x = np.linspace(0, 10, 1000)

# Define "T", i.e functions' period.
T = 2
L = T / 2

# "f(x)" function definition.
def f(x): 
    return np.sin(np.pi * 1000 * x)

# "a" coefficient calculation.
def a(n, L, accuracy = 1000):
    a, b = -L, L
    dx = (b - a) / accuracy
    integration = 0
    for i in np.linspace(a, b, accuracy):
        x = a + i * dx
        integration += f(x) * np.cos((n * np.pi * x) / L)
    integration *= dx
    return (1 / L) * integration

# "b" coefficient calculation.
def b(n, L, accuracy = 1000):
    a, b = -L, L
    dx = (b - a) / accuracy
    integration = 0
    for i in np.linspace(a, b, accuracy):
        x = a + i * dx
        integration += f(x) * np.sin((n * np.pi * x) / L)
    integration *= dx
    return (1 / L) * integration

# Fourier series.   
def Sf(x, L, n = 10):
    a0 = a(0, L)
    sum = 0
    for i in np.arange(1, n + 1):
        sum += ((a(i, L) * np.cos(n * np.pi * x)) + (b(i, L) * np.sin(n * np.pi * x)))
    return (a0 / 2) + sum    

# x axis.
py.plot(x, np.zeros(np.size(x)), color = 'black')

# y axis.
py.plot(np.zeros(np.size(x)), x, color = 'black')

# Original signal.
py.plot(x, f(x), linewidth = 1.5, label = 'Signal')

# Approximation signal (Fourier series coefficients).
py.plot(x, Sf(x, L), color = 'red', linewidth = 1.5, label = 'Fourier series')

# Specify x and y axes limits.
py.xlim([0, 10])
py.ylim([-2, 2])

py.legend(loc = 'upper right', fontsize = '10')

py.show()

...and here is what I get after plotting the result:

I've read the How to calculate a Fourier series in Numpy? and I've implemented this approach already. It works great, but it use the expotential method, where I want to focus on trigonometry functions and the rectangular method in case of calculating the integraions for a_{n} and b_{n} coefficients.

Thank you in advance.

UPDATE (SOLVED)

Finally, here is a working example of the code. However, I'll spend more time on it, so if there is anything that can be improved, it will be done.

from __future__ import division
import numpy as np
import pylab as py

# Define "x" range.
x = np.linspace(0, 10, 1000)

# Define "T", i.e functions' period.
T = 2
L = T / 2

# "f(x)" function definition.
def f(x): 
    return np.sin((np.pi) * x) + np.sin((2 * np.pi) * x) + np.sin((5 * np.pi) * x)

# "a" coefficient calculation.
def a(n, L, accuracy = 1000):
    a, b = -L, L
    dx = (b - a) / accuracy
    integration = 0
    for x in np.linspace(a, b, accuracy):
        integration += f(x) * np.cos((n * np.pi * x) / L)
    integration *= dx
    return (1 / L) * integration

# "b" coefficient calculation.
def b(n, L, accuracy = 1000):
    a, b = -L, L
    dx = (b - a) / accuracy
    integration = 0
    for x in np.linspace(a, b, accuracy):
        integration += f(x) * np.sin((n * np.pi * x) / L)
    integration *= dx
    return (1 / L) * integration

# Fourier series.   
def Sf(x, L, n = 10):
    a0 = a(0, L)
    sum = np.zeros(np.size(x))
    for i in np.arange(1, n + 1):
        sum += ((a(i, L) * np.cos((i * np.pi * x) / L)) + (b(i, L) * np.sin((i * np.pi * x) / L)))
    return (a0 / 2) + sum   

# x axis.
py.plot(x, np.zeros(np.size(x)), color = 'black')

# y axis.
py.plot(np.zeros(np.size(x)), x, color = 'black')

# Original signal.
py.plot(x, f(x), linewidth = 1.5, label = 'Signal')

# Approximation signal (Fourier series coefficients).
py.plot(x, Sf(x, L), '.', color = 'red', linewidth = 1.5, label = 'Fourier series')

# Specify x and y axes limits.
py.xlim([0, 5])
py.ylim([-2.2, 2.2])

py.legend(loc = 'upper right', fontsize = '10')

py.show()

Consider developing your code in a different way, block by block. You should be surprised if a code like this would work at the first try. Debugging is one option, as @tom10 said. The other option is rapid prototyping the code step by step in the interpreter, even better with ipython.

Above, you are expecting that b_1000 is non-zero, since the input f(x) is a sinusoid with a 1000 in it. You're also expecting that all other coefficients are zero right?

Then you should focus on the function b(n, L, accuracy = 1000) only. Looking at it, 3 things are going wrong. Here are some hints.

the multiplication of dx is within the loop. Sure about that?
in the loop, i is supposed to be an integer right? Is it really an integer? by prototyping or debugging you would discover this
be careful whenever you write (1/L) or a similar expression. If you're using python2.7, you're doing likely wrong. If not, at least use a from __future__ import division at the top of your source. Read this PEP if you don't know what I am talking about.

If you address these 3 points, b() will work. Then think of a in a similar fashion.