import math
import numpy

class DomainError(Exception):
    """Exception class for signalling exceptions in these functions."""
    pass


class xtoN:
    def __init__(self,N):
        self.N = N

    def __call__(self,x):
        return pow(x,self.N)


class Gamma_From_Int:
    """Represents the kernel of the integral form of Gamma"""
    def __init__(self,arg):
        self.arg = arg

    def __call__(self,x):
        if x < 0 or x > 1:
            raise DomainError('Trying to integrate Gamma_From_Int outside [0,1]')
        if x == 0:
            return 0.
        if self.arg < 1. and x == 1.: 
            return 0.
        res = math.pow( numpy.log(1./x), self.arg - 1. ) 
        # res  = math.pow(x,self.arg-1.)*numpy.exp(-x)
        return res



def Gamma(x,accuracy=None):
    result  = 1.
    if x == int(x):
        if x > 0:
            i = 1;
            while i < x:
                result *= i
                i      += 1
        else:
            raise DomainError('Gamma function for negative integer argument is undefined')
    else:
        result = 1./x
        factor = 2 * accuracy
        n      = 1
        while abs(1.-factor) > accuracy:
            factor  = math.pow( 1.+1./n, x ) / ( 1. + x/n )
            result *= factor
            n      += 1
    return result

class BesselJ_From_Int:
    def __init__(self,order,arg):
        self.order = int(order)
        self.arg   = arg
        print 'Initialised integrand for BesselJ[',self.order,', ',self.arg,'].'

    def __call__(self,x):
        #################################################################
        # Fix this function with the kernel of the integral form (Q3.b)
        # Use 'raise DomainError' to signal undefined parameter points
        #################################################################
        return 0.1

def BesselJ(order,x,accuracy=None):
    #################################################################
    # Fix this function with the series form (Q.3a)
    # Use 'raise DomainError' to signal undefined parameter points
    #################################################################
    return 0.1



if __name__ == "__main__":
    import scipy.special as sp
    k = 0
    for i in range(0,4): 
        BesselJk = BesselJ(k,i,0.00001)
        print 'J[',k,',',i,'] = ',BesselJk,'(',sp.jn(k,i),')'