#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Examples of plotting with pylab to illustrate Riemann left and right rectangles, and Trapezoidal methods for computing definite integrals of regular functions (of R -> R).

@date: Sun Apr 05 15:00:27 2015.
@author: Lilian Besson for CS101 course at Mahindra Ecole Centrale 2015.
@licence: MIT Licence (http://lbesson.mit-license.org).
"""

from pylab import *

print "Plotting and computing the integral of cos(x) from -4 to 4."
print "Formally, the value is about -1.5136049 (= 2*sin(4))."


# %% Riemann left rectangles rule

def plot_riemann_left(f, xmin, xmax, name_f, n=10):
    """ Plot the function f from xmin to xmax, and n Riemann left rectangles."""
    figure()
    xvalues = linspace(xmin, xmax, 1000)
    yvalues = [f(x) for x in xvalues]
    plot(xvalues, yvalues, 'r-')  # plots the function as a red line
    # Now plot the rectangles
    h = (xmax - xmin) / float(n-1)
    xi = [xmin + i*h for i in xrange(0, n-1)]
    yi = [f(x) for x in xi]  # left rectangles!
    for x, y in zip(xi, yi):
        plot([x, x, x+h, x+h, x], [0, y, y, 0, 0], 'b-')
    # Finally, computing the area
    area = round(h * sum(yi), 4)
    title("Riemman left rectangles for " + name_f + " with " + str(n) + " rectangles. Area is " + str(area))
    savefig("riemannleft.png")
    return area


# One example
print "With left rectangles:", plot_riemann_left(cos, -4, 4, "cos", n=15)


# %% Riemann right rectangles rule

def plot_riemann_right(f, xmin, xmax, name_f, n=10):
    """ Plot the function f from xmin to xmax, and n Riemann right rectangles."""
    figure()
    xvalues = linspace(xmin, xmax, 1000)
    yvalues = [f(x) for x in xvalues]
    plot(xvalues, yvalues, 'r-')  # plots the function as a red line
    # Now plot the rectangles
    h = (xmax - xmin) / float(n-1)
    xi = [xmin + i*h for i in xrange(0, n-1)]
    yi = [f(x+h) for x in xi]  # right rectangles!
    for x, y in zip(xi, yi):
        plot([x, x, x+h, x+h, x], [0, y, y, 0, 0], 'b-')
    # Finally, computing the area
    area = round(h * sum(yi), 4)
    title("Riemman right rectangles for " + name_f + " with " + str(n) + " rectangles. Area is " + str(area))
    savefig("riemannright.png")
    return area


# One example
print "With right rectangles:", plot_riemann_right(cos, -4, 4, "cos", n=15)


# %% Riemann center rectangles rule

def plot_riemann_center(f, xmin, xmax, name_f, n=10):
    """ Plot the function f from xmin to xmax, and n Riemann center rectangles."""
    figure()
    xvalues = linspace(xmin, xmax, 1000)
    yvalues = [f(x) for x in xvalues]
    plot(xvalues, yvalues, 'r-')  # plots the function as a red line
    # Now plot the rectangles
    h = (xmax - xmin) / float(n-1)
    xi = [xmin + i*h for i in xrange(0, n-1)]
    yi = [f(x + h/2.0) for x in xi]  # center rectangles!
    for x, y in zip(xi, yi):
        plot([x, x, x+h, x+h, x], [0, y, y, 0, 0], 'b-')
    # Finally, computing the area
    area = round(h * sum(yi), 4)
    title("Riemman center rectangles for " + name_f + " with " + str(n) + " rectangles. Area is " + str(area))
    savefig("riemanncenter.png")
    return area


# One example
print "With center rectangles:", plot_riemann_center(cos, -4, 4, "cos", n=15)


# %% Trapezoidal rule

def plot_trapez(f, xmin, xmax, name_f, n=10):
    """ Plot the function f from xmin to xmax, and n trapezoides (as used for the Trapezoidal rule)."""
    figure()
    xvalues = linspace(xmin, xmax, 1000)
    yvalues = [f(x) for x in xvalues]
    plot(xvalues, yvalues, 'r-')  # plots the function as a red line
    # Now plot the trapezoides
    h = (xmax - xmin) / float(n-1)
    xi = [xmin + i*h for i in xrange(0, n-1)]
    # trapezoides !
    for x in xi:
        plot([x, x, x+h, x+h, x], [0, f(x), f(x+h), 0, 0], 'b-')  # trapez
        plot([x, x, x+h, x+h, x], [0, 0.5 * (f(x) + f(x+h)), 0.5 * (f(x) + f(x+h)), 0, 0], 'g:')  # equivalent rectangle
    # Finally, computing the area
    yi = [0.5 * (f(x) + f(x+h)) for x in xi]
    area = round(h * sum(yi), 4)
    title("Trapezoidal rule for " + name_f + " with " + str(n) + " trapezoides. Area is " + str(area))
    savefig("trapezoides.png")
    return area


# One example
print "With trapez:", plot_trapez(cos, -4, 4, "cos", n=15)