#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Sep 17 14:04:45 2024

@author: sam
"""

from scipy.special import p_roots

def quadz(f,a,b,tol): # functions from last week
    Iold = 0
    for i in range(1,100):
        x,w = p_roots(i)
        I = (b-a)/2*sum(w*f((b-a)/2*x+(a+b)/2))
        if i > 3 and abs(I - Iold) < tol:
            return I
        Iold = I
    return I

def I(f,a,b,n): # functions from last week
    I = 0
    x,w = p_roots(n)
    for i in range(len(x)):
        I += w[i]*f((b-a)/2*x[i]+(a+b)/2)
    I *= (b-a)/2
    return I

def f(x,y):
    return (x**2)*y



def squad2d(f,a,b,c,d):
    def yint(f,x,c,d): # interior of integral
        def fy(y): # single valued function of f, where x is constant
            return f(x,y)
        return I(fy,c,d,13) # integrates f(y) between c and d
    def fx(x): # exterior of integral, function in terms of only f
        return yint(f,x,c,d)
    return I(fx,a,b,13) # integrates f(x) between a and b



print(squad2d(f,0,1,0,2))