#!/usr/bin/python3

import numpy as np
import sys

# c TB: taken from https://en.wikipedia.org/wiki/Simpson%27s_rule

'''
program simpsonI_non_cum.py
	- "running" integration of user provided data (2-column format) using Simpson 
	method. 

      	- example of usage:

	python3 simpsonI_non_cum.py grad.dat 

	Here, the integration variable and integrand are stored in file grad.dat, 
	result is written on standard output
'''

def read_file(name):
    old = open(name, 'r')
    xvect = []
    yvect = []
    dyvect = []
    for line in old.readlines():
        prefield = line.split()
        xvect.append(float(prefield[0]))
        yvect.append(float(prefield[1]))
    old.close()
    return np.array(xvect), np.array(yvect)


def integration_err(x, f, df, step):
    err = 0.
    for i in range(len(f)):
        f_ = 1*f
        f_[i] += step
        arreaP = simpson_nonuniform(x, f_)
        f_ = 1*f
        f_[i] -= step
        arreaM = simpson_nonuniform(x, f_)
        dA = (arreaP-arreaM)/(2*step)
        err += (dA*df[i])**2
    err = err**0.5
    return err


def simpson_nonuniform(x, f):
    """
    Simpson rule for irregularly spaced data.

        Parameters
        ----------
        x : list or np.array of floats
                Sampling points for the function values
        f : list or np.array of floats
                Function values at the sampling points

        Returns
        -------
        float : approximation for the integral
    """
    N = len(x) - 1
    h = np.diff(x)

    result = 0.0
    for i in range(1, N, 2):
        hph = h[i] + h[i - 1]
        result += f[i] * (h[i]**3 + h[i - 1]**3
                          + 3. * h[i] * h[i - 1] * hph)\
            / (6 * h[i] * h[i - 1])
        result += f[i - 1] * (2. * h[i - 1]**3 - h[i]**3
                              + 3. * h[i] * h[i - 1]**2)\
            / (6 * h[i - 1] * hph)
        result += f[i + 1] * (2. * h[i]**3 - h[i - 1]**3
                              + 3. * h[i - 1] * h[i]**2)\
            / (6 * h[i] * hph)

    if (N + 1) % 2 == 0:
        result += f[N] * (2 * h[N - 1]**2
                          + 3. * h[N - 2] * h[N - 1])\
            / (6 * (h[N - 2] + h[N - 1]))
        result += f[N - 1] * (h[N - 1]**2
                              + 3*h[N - 1] * h[N - 2])\
            / (6 * h[N - 2])
        result -= f[N - 2] * h[N - 1]**3\
            / (6 * h[N - 2] * (h[N - 2] + h[N - 1]))
    return result


name = sys.argv[1]
x, f = read_file(name)

step = 1e-3

dA = np.zeros(len(f), float)

arrea = simpson_nonuniform(x, f)
# print arrea
# print ""

print(x[0], 0.)
for i in range(1, len(f)):
    f_ = f[:i+1]
    x_ = x[:i+1]
    arrea_ = simpson_nonuniform(x_, f_)
    print(x[i], arrea_)
