'''This file contains gaussian quadrature formula of order 2, 4 and 6
These functions are called by the script python13_1.py.

Input:
    x... grid points of equidistant partition of interval.
    fcn... the function to be integrated

Output: 
    approx... approximation to the integral by means of quadrature formula
'''

import numpy as np

def gauss_ord2(x, fcn):
    n = len(x)-1  # number of subintervals of partition
    h = x[1] - x[0]  # length of subintervals (stepsize)
    approx = 0
    
    c1 = .5
                  
    for j in range(n):
        Y1 = fcn(x[j] + h*c1)
        
        approx += h*Y1
        
    return approx


def gauss_ord4(x, fcn):
    n = len(x)-1  # number of subintervals of partition
    h = x[1] - x[0]  # length of subintervals (stepsize)
    approx = 0
    
    c1 = .5 - np.sqrt(3)/6
    c2 = .5 + np.sqrt(3)/6
                  
    for j in range(n):
        Y1 = fcn(x[j] + h*c1)
        Y2 = fcn(x[j] + h*c2)
        
        approx += 0.5*h*(Y1 + Y2)
        
    return approx


def gauss_ord6(x, fcn):
    n = len(x)-1  # number of subintervals of partition
    h = x[1] - x[0]  # length of subintervals (stepsize)
    approx = 0
    
    c1 = .5 - np.sqrt(15)/10
    c2 = .5
    c3 = .5 + np.sqrt(15)/10
                  
    for j in range(n):
        Y1 = fcn(x[j] + h*c1)
        Y2 = fcn(x[j] + h*c2)
        Y3 = fcn(x[j] + h*c3)
        
        approx += h/18*(5*Y1 + 8*Y2 + 5*Y3)
        
    return approx