def double_integral(func, limits, res=1000):
    t = time.clock()
    t1 = time.clock()
    t2 = time.clock()
    s = 0
    a, b = limits[0], limits[1]
    outer_values = np.linspace(a, b, res)
    c_is_func = callable(limits[2])
    d_is_func = callable(limits[3])
    for y in outer_values:
        if c_is_func:
            c = limits[2](y)
        else:
            c = limits[2]
        if d_is_func:
            d = limits[3](y)
        else:
            d = limits[3]
        dA = ((b - a) / res) * ((d - c) / res)
        inner_values = np.linspace(c, d, res)
        for x in inner_values:
            t2 = time.clock() - t2
            s += func(x, y) * dA
        t1 = time.clock() - t1
     t = time.clock() - t
     return s, t, t1 / res, t2 / res**2

def f(x, y):
    if (4 - y**2 - x**2) < 0:
        return 0              #This is to avoid taking the root of negarive #'s
    return np.sqrt(4 - y**2 - x**2)

def c(y):
    return np.sqrt(2 * y - y**2)

def d(y):
    return np.sqrt(4 - y**2)
# b d
# S S f(x,y) dx dy
# a c
a, b, = 0, 2
print(double_integral(f, [a, b, c, d]))