#********************************
def int_1d_tpz(f,a,b,rel0):
#********************************
#This subrountine calculates one dimentional
#integral with fixed limits a and b in trapezoidal
#approximation int_1d_tpz = \int_a^b  f(x)dx

#f     -  integrand - user defined function f(x), input
#a,b   - limits of integration - real, input
#rel0  - reaquired relative error - real imput

#************ EXAMPLE ****************
#def fun(x):
#    f = exp(-x*x)
#    return f
#rel = 0.0001
#a = 0.0
#b = 2.0
#sint = int_1d_tpz(fun,a,b,rel)
#print(sint)
#**************************************
#V.1
#Misak Sargsian
#FIU, 22-Feb-2017

    N1 = 10
    rel = 2.0*rel0
    h1 = (b-a)/N1
    s1 = 0.5*(f(a) + f(b))
    for k in range(1,N1):
        s1 += f(a+k*h1)
    val1 = h1*s1



    while(rel>rel0):    
        N2 = N1*2
        h2 = h1/2.0
        s2 = 0.0
        for k in range(1,N2,2):
            s2 += f(a+k*h2)
        val2  = 0.5*val1 + h2*s2
        rel = abs((val2-val1)/3.0)/abs(val2)
        val1 = val2
        N1 = N2
        h1 = h2
    return val2  
        
#********************************
def int_1d_simps(f,a,b,rel0):
#********************************
#This subrountine calculates one dimentional
#integral with fixed limits a and b in Simpson
#approximation int_1d_tpz = \int_a^b  f(x)dx

#f     -  integrand - user defined function f(x), input
#a,b   - limits of integration - real, input
#rel0  - reaquired relative error - real imput

#************ EXAMPLE ****************
#def fun(x):
#    f = exp(-x*x)
#    return f
#rel = 0.0001
#a = 0.0
#b = 2.0
#sint = int_1d_tpz(fun,a,b,rel)
#print(sint)
#**************************************
#V.1
#Misak Sargsian
#FIU, 22-Feb-2017

    N1 = 10
    rel = 2.0*rel0
    h1 = (b-a)/N1
    s1 = 1.0/3.0* (f(a) + f(b))
    t1 = 2.0/3.0*f(a+(N1-1)*h1)
    for k in range(2,N1,2):
        s1 += 2.0/3.0*f(a+k*h1)
        t1 += 2.0/3.0*f(a+(k-1)*h1)
    val1 = h1*(s1 + 2*t1)
    
    while(rel>rel0):    
        N2 = N1*2
        h2 = h1/2.0
        s2 = s1+t1
        t2 = 0.0
        for k in range(1,N2,2):
            t2 += 2.0/3.0*f(a+k*h2)
        val2  = h2*(s2 + 2.0*t2)
        rel = abs(val2-val1)/abs(val2)
        t1 = t2
        s1 = s2        
        val1 = val2
        N1 = N2
        h1 = h2       
    return val2               
  

#********************************
def int_2d_tpz(fun,a,b,fly,fuy,rel0):
#********************************
#This subrountine calculates two dimentional
#integral using trapezoidal methos
#with fixed limits a and b for the outher variable x
#and fl(x) and fu(x) for inner variable y
#\int_a^b(\int_fl(x)^fu(x) fun(x,y)dy)dx    
#f(x,y) -  integrand - user defined function f(x,y), input
#a,b   - limits of integration  by x- real, input
#fly(x), fuy(y) - limits of integration by y, (functions of x), input 
#rel0  - reaquired relative error - real imput

#************ EXAMPLE ****************
#def fun(x,y):
#    f = x*y
#    return f
#rel = 0.0001
#a = 0.0
#b = 2.0
#def fl(x)
#    ymin = 2.0*x
#    return ymin
#def fu(x)
#    ymax = 4.0*x
#    return ymax
#rel0 = 0.0000001
#a = 0.0
#b = 2.0
#
#r = int_2d_tpz(fun2,a,b,fl,fu,rel0)
#print (r)
#**************************************
#V.1
#Misak Sargsian
#FIU, 27-Feb-2017
    def funx(x):
        N1y = 10
        rely = 2.0*rel0
        ay = fly(x)
        by = fuy(x)
        h1y = (by - ay)/N1y
        s1y = 0.5*(fun(x,ay) + fun(x,by))
        for k in range(1,N1y):
            s1y += fun(x,ay+k*h1y)
        val1y = h1y*s1y

        while(rely>rel0):    
            N2y = N1y*2
            h2y = h1y/2.0
            s2y = 0.0
            for k in range(1,N2y,2):
                s2y += fun(x,ay+k*h2y)
            val2y  = 0.5*val1y + h2y*s2y

            if val2y==0.0:
               rely = abs((val2y-val1y)/3.0)
            else:   
               rely = abs((val2y-val1y)/3.0)/abs(val2y)
            
            val1y = val2y
            N1y = N2y
            h1y = h2y
        return val1y  

    sint = int_1d_tpz(funx,a,b,rel0)
    return sint






#********************************
def int_2d_simps(fun,a,b,fly,fuy,rel0):
#********************************
#This subrountine calculates two dimentional
#integral with Simpson method
#with fixed limits a and b for the outher variable x
#and fl(x) and fu(x) for inner variable y
#\int_a^b(\int_fl(x)^fu(x) fun(x,y)dy)dx    
#f(x,y) -  integrand - user defined function f(x,y), input
#a,b   - limits of integration  by x- real, input
#fly(x), fuy(y) - limits of integration by y, (functions of x), input 
#rel0  - reaquired relative error - real imput

#************ EXAMPLE ****************
#def fun(x,y):
#    f = x*y
#    return f
#rel = 0.0001
#a = 0.0
#b = 2.0
#def fl(x)
#    ymin = 2.0*x
#    return ymin
#def fu(x)
#    ymax = 4.0*x
#    return ymax
#rel0 = 0.0000001
#a = 0.0
#b = 2.0
#
#r = int_2d_simps(fun2,a,b,fl,fu,rel0)
#print (r)
#**************************************
#V.1
#Misak Sargsian
#FIU, 27-Feb-2017
    def funx(x):
        N1y = 10
        rely = 2.0*rel0
        ay = fly(x)
        by = fuy(x)
        h1y = (by - ay)/N1y
        s1y = 1.0/3.0* (fun(x,ay) + fun(x,by))
        t1y = 2.0/3.0*fun(x,ay+(N1y-1)*h1y)
        for k in range(2,N1y,2):
            s1y += 2.0/3.0*fun(x,ay+k*h1y)
            t1y += 2.0/3.0*fun(x,ay+(k-1)*h1y)
        val1y = h1y*(s1y + 2*t1y)

        while(rely>rel0):    
            N2y = N1y*2
            h2y = h1y/2.0
            s2y = s1y+t1y
            t2y = 0.0
            for k in range(1,N2y,2):
                t2y += 2.0/3.0*fun(x,ay+k*h2y)
            val2y  = h2y*(s2y + 2.0*t2y)
            if val2y==0.0:
                rely = abs((val2y-val1y)/15.0)
            else:
                rely = abs((val2y-val1y)/15.0)/abs(val2y)
            
        t1y = t2y
        s1y = s2y        
        val1y = val2y
#        print(val1y,'val1')
        N1y = N2y
        h1y = h2y        
#    print(N2,rel)     
        return val1y               
  
    sint = int_1d_simps(funx,a,b,rel0)
    return sint