# include # include # include # include # include "annulus_monte_carlo.h" /******************************************************************************/ double disk01_monomial_integral ( int e[2] ) /******************************************************************************/ /* Purpose: disk01_monomial_integral() returns monomial integrals in the unit disk in 2D. Discussion: The integration region is X^2 + Y^2 <= 1. The monomial is F(X,Y) = X^E(1) * Y^E(2). Licensing: This code is distributed under the GNU LGPL license. Modified: 03 January 2014 Author: John Burkardt Reference: Philip Davis, Philip Rabinowitz, Methods of Numerical Integration, Second Edition, Academic Press, 1984, page 263. Parameters: Input, int E[2], the exponents of X and Y in the monomial. Each exponent must be nonnegative. Output, double DISK01_MONOMIAL_INTEGRAL, the integral. */ { double arg; int i; double integral; const double r = 1.0; double s; if ( e[0] < 0 || e[1] < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "disk01_monomial_integral(): Fatal error!\n" ); fprintf ( stderr, " All exponents must be nonnegative.\n" ); fprintf ( stderr, " E[0] = %d\n", e[0] ); fprintf ( stderr, " E[1] = %d\n", e[1] ); exit ( 1 ); } if ( ( e[0] % 2 ) == 1 || ( e[1] % 2 ) == 1 ) { integral = 0.0; } else { integral = 2.0; for ( i = 0; i < 2; i++ ) { arg = 0.5 * ( double ) ( e[i] + 1 ); integral = integral * tgamma ( arg ); } arg = 0.5 * ( double ) ( e[0] + e[1] + 2 ); integral = integral / tgamma ( arg ); } /* Adjust the surface integral to get the volume integral. */ s = e[0] + e[1] + 2; integral = integral * pow ( r, s ) / ( double ) ( s ); return integral; }