'''Experiment 17.2: randomly chosen Riemann sums.

The program computes Riemann sums for the function f(x,y) = x**2+y**2 on
the rectangle [0,1]x[0,1], plots a grid with random intermediate points
and the approximation of the graph by the top faces of the cuboids.

Input:
    n... Number of summands (n <= 50 is recommended)

Output:
    sum... Approximation of the integral of f(x,y) = x**2+y**2 on [0,1]x[0,1].

Typical call of program:
    from python17_1 import riemann_sum
    riemann_sum(10)      # riemann_sum(n)

'''

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import mpl_toolkits.mplot3d as a3

def riemann_sum(n = 10):
    
    if n > 50:
        print('To limit computing time, n <= 50 is recommended - or '+'else                 comment Figure(2)')
    
    t = np.arange(0, 1, 1./n)
    x, y = np.meshgrid(t,t)
    
    '''Here we create random intermediate points (xi, yi) such that
    x[i,j] <= xi[i,j] <= x[i,j+1] and y[i,j] <= yi[i,j] <= y[i+1,j]
    for all i = 0...n-1, j = 0...n-1.
    '''
    xi = x + np.random.uniform(0,1./n,x.shape)
    yi = y + np.random.uniform(0,1./n,y.shape)
    
    fx = xi**2 + yi**2
    r_sum = sum(sum(fx))/n**2
    
    '''We plot the function x**2 + y**2 and the upper faces of the cuboids representing the summands of the Riemann sum
    '''  
    mpl.rcParams['xtick.labelsize'] = 8
    mpl.rcParams['lines.linewidth'] = 0.8
    mpl.rcParams['lines.markersize'] = 4

    plt.figure(1)
    plt.grid()
    plt.plot(xi,yi,'r.')
    plt.xticks(np.linspace(0,1,11))
    plt.yticks(np.linspace(0,1,11))
    plt.title('Top faces of cuboids')
    
    '''3D-image of top faces of the cuboids''' 
    t = np.linspace(0, 1, n+1)
    ax = a3.Axes3D(plt.figure(2))
    for i in range(n):
        for j in range(n):
            vtx = np.array(
                   [[t[i], t[j], fx[i,j]],
                   [t[i+1], t[j], fx[i,j]],
                   [t[i+1], t[j+1], fx[i,j]],
                   [t[i], t[j+1], fx[i,j]],
                   [t[i], t[j], fx[i,j]],])
            boxes = a3.art3d.Poly3DCollection([vtx])
            ax.add_collection3d(boxes)
    ax.set_zlim([0,2])
    plt.xticks(np.linspace(0,1,11))
    plt.yticks(np.linspace(0,1,11))
    ax.set_zticks([0,1,2])
    ax.view_init(20,220)
    
    '''Boundary of the surface z = x^2 + y^2'''
    t = np.linspace(0, 1, 100)
    ax.plot(np.zeros(t.shape), t, t**2, 'k')
    ax.plot(t, np.zeros(t.shape), t**2, 'k')
    ax.plot(np.ones(t.shape), t, 1 + t**2, 'k')
    ax.plot(t, np.ones(t.shape), 1 + t**2, 'k')
    
    plt.show()
    return r_sum