Here is the code used:

import numpy as np
GaussTable = [[[0], [2]], [[-1/np.sqrt(3), 1/np.sqrt(3)], [1, 1]], [[-np.sqrt(3/5), 0, np.sqrt(3/5)], [5/9, 8/9, 5/9]], [[-0.861136, -0.339981, 0.339981, 0.861136], [0.347855, 0.652145, 0.652145, 0.347855]], [[-0.90618, -0.538469, 0, 0.538469, 0.90618], [0.236927, 0.478629, 0.568889, 0.478629, 0.236927]], [[-0.93247, -0.661209, -0.238619, 0.238619, 0.661209, 0.93247], [0.171324, 0.360762, 0.467914, 0.467914, 0.360762,  0.171324]]]



def f(r,th):
    return r**2*np.cos(th)
def IGAL(f, n1, r1, r2, n2,th1,th2):
  n1 = int(n1)
  n2=int(n2)
  return sum([sum([(r2 - r1)/2*GaussTable[n1 - 1][1][i]*(th2 - th1)/2*GaussTable[n2 - 1][1][j]*f((r2 - r1)/2*(GaussTable[n1 - 1][0][i] + 1) + r1,(th2 - th1)/2*(GaussTable[n2 - 1][0][j] + 1) + th1) for i in range(n1)]) for j in range(n2)])



def main(inputs):
    n1=inputs['n_1']
    n2=inputs['n_2']
    r1=inputs['r_1']
    r2=inputs['r_2']
    th1=inputs['th_1']*np.pi
    th2=inputs['th_2']*np.pi
    n1=2
    n2=2
    s=IGAL(f,n1,r1,r2,n2,th1,th2)
    #print(s)
    return {"Sol":s,"n1":n1,"n2":n2,"r1":r1,"r2":r2,"th1":inputs['th_1'],"th2":inputs['th_2']}