вопрос рекурсии (2-я квадратура)

Итак, у меня есть программа, которая реализует адаптивное 2D трапециевидное правило для функции x ^ 2 + y ^ 2 < 1, но кажется, что рекурсия не работает — программа здесь представляет собой модифицированную форму (рабочего) трапециевидного 1D-метода, поэтому я не уверен, где код ломается, он должен возвращать PI:

double trapezoidal(d_fp_d f,
double a, double b,
double c, double d) { //helper function
return 0.25*(b-a)*(d-c)*
(f(a, c)+f(a, d) +
f(b, c)+f(b, d));
}

double atrap( double a, double b, double c, double d, d_fp_d f, double tol )
{// helper function

return atrap1(a, b, c, d, f, tol );
}
double atrap1( double a, double b, double c, double d, d_fp_d f, double tol)
{
//implements 2D trap rule
static int level = 0;
const static int minLevel = 4;
const static int maxLevel = 30;
++level;
double m1 = (a + b)/2.0;
double m2 = (c + d)/2.0;
double coarse = trapezoidal(f,a,b,c,d);
double fine =
trapezoidal(f, a, m1, c, m2)
+ trapezoidal(f, a, m1, m2, d)
+ trapezoidal(f, m1, b, c, m2)
+ trapezoidal(f, m1, b, m2, d);
++fnEvals;
if( level< minLevel
|| ( abs( fine - coarse ) > 3.0*tol && level < maxLevel ) ){

fine =  atrap1( a, m1, c, m2, f,tol/4.0)
+ atrap1( a, m1, m2, d, f, tol/4.0)
+ atrap1(m1, b, c, m2, f, tol/4.0)
+ atrap1(m1, b, m2, d, f,tol/4.0);

}

--level;
return fine;
}

где функция определяется как

double ucircle( double x, double y)
{
return x*x + y*y < 1 ? 1.0 : 0.0;
}

и моя основная функция

int main()
{
double a, b, c, d;
cout << "Enter a: " <<endl;
cin >> a;
cout << "Enter b: " <<endl;
cin >> b;
cout << "Enter c: " <<endl;
cin >> c;
cout << "Enter d: " <<endl;
cin >> d;

cout << "The approximate integral is: " << atrap( a, b, c, d, ucircle, 1.0e-5) << endl;

return 0;
}