####################################################################
#
# A class to integrate functions f(x) of one parameter x
#
#   I_f(x0,xN) = int_x0^x1 dx f(x)
#
# through quadratures.
#
# In general, a solution to such a problem is provided by
# discretisation in the integration variable, replacing
# the integration with a sum, 
#
# I_f(x0,xN) = sum_{i=1}^{N-1} [Deltax f(x_i)]
#
# i.e. by not taking the limit Deltax -> 0.
#
# Different methods will be supplemented by the students through 
# the homework exercise, in particular the simple Newton-Cotes,
# Trapezoidal and Simpson rules.
#
#
#
# In the initialisation of this class, through 
#
#   __init__(self,kernel,algorithm),
#
# the function f needs to be provided ("kernel"), and the 
# algorithm for the solution must be picked - the choices above 
# are selected through the names "Newton-Cotes", "Trapez", 
# "Simpson", etc.
# Note that in the implementation here, the kernel must provide 
# a method 
#  
#   value(x), 
#
# which depends on x only.
#
# In order to integrate the function with a FIXED number of steps
# N the method
#  
#  integrate_Nsteps(x0,x1,N)
#
# is invoked.  Here the arguments denote the boundary conditions
# x0 and x1 and the number of steps N such that the fixed stepsize
# delta_x = (x1-x0)/(N-1).
# The result of this integration is returned.
#


class Integrator:
    def __init__(self,kernel,algorithm):
        self.kernel  = kernel
        if algorithm == 'Newton-Cotes':
            self.algorithm = self.NCStep
        elif algorithm == 'Trapez':
            self.algorithm = self.TrapezStep
        elif algorithm == 'Simpson':
            self.algorithm = self.SimpsonStep
        else:
            raise Exception('ERROR: Unknown algorithm %s' % algorithm)

    def initialise(self,x0,x1,N):
        self.x0      = x0
        self.x1      = x1
        self.N       = N
        self.deltax  = (x1-x0)/float(N)

    def integrate_Nsteps(self,x0,x1,N):
        self.initialise(x0,x1,N)
        integral = self.algorithm()
        return integral

    def integrate_Accuracy(self,x0,x1,error):
        ##################################################
        # Fix this function (Q.2)
        ##################################################
        result = 0.123456789
        final_accuracy = 0
        return result, final_accuracy


    # Different algorithm choices, can be selected through the algorithm
    # argument in __init__

    def NCStep(self):
        ##################################################
        # Fix this function (Q.1a)
        ##################################################
        return 0.1

    def TrapezStep(self):
        ##################################################
        # Fix this function (Q.1b)
        ##################################################
        return 0.1

    def SimpsonStep(self):
        ##################################################
        # Fix this function (Q.1c)
        ##################################################
        return 0.1


if __name__ == "__main__":
    from Functions import xtoN
    from Functions import Gamma_From_Int
    from Functions import Gamma
    from Functions import BesselJ_From_Int
    from Functions import BesselJ


    print '###### Integrate x^4 from 0 to 1 (result = 0.2) ######'
    kernel     = xtoN(4)
    #integrator = Integrator(kernel,'Newton-Cotes')
    integrator = Integrator(kernel,'Simpson')
    #print integrator.integrate_Nsteps(0,1,9)
    result,accu = integrator.integrate_Accuracy(0,1,0.0000001)
    print 'Result = ',result,' +/- ',(accu*100.),'%.'



    from numpy import sqrt, pi
    print '###### Testing gamma function for x=0.5 (result = %s) ######' % sqrt(pi)
    argument    = 0.5
    accuracy    = 1.e-3
    kernel      = Gamma_From_Int(argument)
    integrator  = Integrator(kernel,'Simpson')
    result,accu = integrator.integrate_Accuracy(0.,1.,accuracy)
    seriesres   = Gamma(argument,accuracy)
    dev         = abs(2.*(seriesres-result)/(seriesres+result))
    if (result==0.) and (seriesres==0.): dev = 0.
    print 'Result = ',result,' +/- ',(accu*100.),'% vs ',seriesres,' deviation = ',dev



    from scipy.special import jn
    print '###### Testing Bessel J_2(5) (result = %s) ######' % jn(2,5)
    order       = 2
    argument    = 5
    accuracy    = 1.e-3
    kernel      = BesselJ_From_Int(order,argument)
    integrator  = Integrator(kernel,'Simpson')
    result,accu = integrator.integrate_Accuracy(0.,pi,accuracy)
    seriesres   = BesselJ(order,argument,accuracy)
    dev         = abs(2.*(seriesres-result)/(seriesres+result))
    if (result==0.) and (seriesres==0.): dev = 0.
    print 'Result = ',result,' +/- ',(accu*100.),'% vs ',seriesres,' deviation = ',dev