# include # include # include # include "annulus_monte_carlo.h" /******************************************************************************/ double annulus_area ( double center[], double r1, double r2 ) /******************************************************************************/ /* Purpose: annulus_area() returns the area of an annulus in 2D. Discussion: A circular annulus with center (XC,YC), inner radius R1 and outer radius R2, is the set of points (X,Y) so that R1^2 <= (X-XC)^2 + (Y-YC)^2 <= R2^2 Licensing: This code is distributed under the MIT license. Modified: 06 July 2018 Author: John Burkardt Parameters: Input, double CENTER[2], coordinates of the center. This information is not actually needed to compute the area. Input, double R1, R2, the inner and outer radius of the disk. 0.0 <= R1 <= R2. Output, double ANNULUS_AREA, the area of the annulus. */ { double area; const double r8_pi = 3.141592653589793; if ( r1 < 0.0 ) { printf ( "\n" ); printf ( "ANNULUS_AREA - Fatal error!\n" ); printf ( " Inner radius R1 < 0.0.\n" ); exit ( 1 ); } if ( r2 < r1 ) { printf ( "\n" ); printf ( "ANNULUS_AREA - Fatal error!\n" ); printf ( " Outer radius R2 < R1 = inner radius.\n" ); exit ( 1 ); } area = r8_pi * ( r2 + r1 ) * ( r2 - r1 ); return area; } /******************************************************************************/ double *annulus_sample ( double center[], double r1, double r2, int n, int *seed ) /******************************************************************************/ /* Purpose: annulus_sample() uniformly samples a general annulus in 2D. Licensing: This code is distributed under the MIT license. Modified: 06 July 2018 Author: John Burkardt Parameters: Input, double CENTER[2], coordinates of the center. Input, double R1, R2, the inner and outer radius. Input, int N, the number of points. Input/output, int *SEED, a seed for the random number generator. Output, double X[2*N], the points. */ { int i; int j; double *p; double *r; const double r8_pi = 3.141592653589793; double *theta; if ( r1 < 0.0 ) { printf ( "\n" ); printf ( "ANNULUS_SAMPLE - Fatal error!\n" ); printf ( " Inner radius R1 < 0.0.\n" ); exit ( 1 ); } if ( r2 < r1 ) { printf ( "\n" ); printf ( "ANNULUS_SAMPLE - Fatal error!\n" ); printf ( " Outer radius R2 < R1 = inner radius.\n" ); exit ( 1 ); } theta = r8vec_uniform_01_new ( n, seed ); for ( i = 0; i < n; i++ ) { theta[i] = theta[i] * 2.0 * r8_pi; } r = r8vec_uniform_01_new ( n, seed ); for ( i = 0; i < n; i++ ) { r[i] = sqrt ( ( 1.0 - r[i] ) * r1 * r1 + r[i] * r2 * r2 ); } p = ( double * ) malloc ( 2 * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { p[0+j*2] = center[0] + r[j] * cos ( theta[j] ); p[1+j*2] = center[1] + r[j] * sin ( theta[j] ); } free ( r ); free ( theta ); return p; } /******************************************************************************/ 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 MIT 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; } /******************************************************************************/ double *monomial_value ( int m, int n, int e[], double x[] ) /******************************************************************************/ /* Purpose: monomial_value() evaluates a monomial. Discussion: This routine evaluates a monomial of the form product ( 1 <= i <= m ) x(i)^e(i) where the exponents are nonnegative integers. Note that if the combination 0^0 is encountered, it should be treated as 1. Licensing: This code is distributed under the MIT license. Modified: 08 May 2014 Author: John Burkardt Parameters: Input, int M, the spatial dimension. Input, int N, the number of points at which the monomial is to be evaluated. Input, int E[M], the exponents. Input, double X[M*N], the point coordinates. Output, double MONOMIAL_VALUE[N], the value of the monomial. */ { int i; int j; double *v; v = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { v[j] = 1.0; } for ( i = 0; i < m; i++ ) { if ( 0 != e[i] ) { for ( j = 0; j < n; j++ ) { v[j] = v[j] * pow ( x[i+j*m], e[i] ); } } } return v; } /******************************************************************************/ double r8_uniform_01 ( int *seed ) /******************************************************************************/ /* Purpose: r8_uniform_01() returns a pseudorandom R8 scaled to [0,1]. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) r8_uniform_01 = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. If the initial seed is 12345, then the first three computations are Input Output R8_UNIFORM_01 SEED SEED 12345 207482415 0.096616 207482415 1790989824 0.833995 1790989824 2035175616 0.947702 Licensing: This code is distributed under the MIT license. Modified: 11 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation edited by Jerry Banks, Wiley Interscience, page 95, 1998. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. P A Lewis, A S Goodman, J M Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input/output, int *SEED, the "seed" value. Normally, this value should not be 0. On output, SEED has been updated. Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { int i4_huge = 2147483647; int k; double r; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r = ( ( double ) ( *seed ) ) * 4.656612875E-10; return r; } /******************************************************************************/ double *r8vec_normal_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: r8vec_normal_01_new() returns a unit pseudonormal R8VEC. Discussion: The standard normal probability distribution function (PDF) has mean 0 and standard deviation 1. Licensing: This code is distributed under the MIT license. Modified: 06 August 2013 Author: John Burkardt Parameters: Input, int N, the number of values desired. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_NORMAL_01_NEW[N], a sample of the standard normal PDF. Local parameters: Local, double R[N+1], is used to store some uniform random values. Its dimension is N+1, but really it is only needed to be the smallest even number greater than or equal to N. Local, int X_LO, X_HI, records the range of entries of X that we need to compute. */ { int i; int m; const double pi = 3.141592653589793; double *r; double *x; int x_hi; int x_lo; x = ( double * ) malloc ( n * sizeof ( double ) ); /* Record the range of X we need to fill in. */ x_lo = 1; x_hi = n; /* If we need just one new value, do that here to avoid null arrays. */ if ( x_hi - x_lo + 1 == 1 ) { r = r8vec_uniform_01_new ( 2, seed ); x[x_hi-1] = sqrt ( - 2.0 * log ( r[0] ) ) * cos ( 2.0 * pi * r[1] ); free ( r ); } /* If we require an even number of values, that's easy. */ else if ( ( x_hi - x_lo + 1 ) % 2 == 0 ) { m = ( x_hi - x_lo + 1 ) / 2; r = r8vec_uniform_01_new ( 2*m, seed ); for ( i = 0; i <= 2*m-2; i = i + 2 ) { x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * pi * r[i+1] ); x[x_lo+i ] = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * pi * r[i+1] ); } free ( r ); } /* If we require an odd number of values, we generate an even number, and handle the last pair specially, storing one in X(N). */ else { x_hi = x_hi - 1; m = ( x_hi - x_lo + 1 ) / 2 + 1; r = r8vec_uniform_01_new ( 2*m, seed ); for ( i = 0; i <= 2*m-4; i = i + 2 ) { x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * pi * r[i+1] ); x[x_lo+i ] = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * pi * r[i+1] ); } i = 2*m - 2; x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * pi * r[i+1] ); free ( r ); } return x; } /******************************************************************************/ double r8vec_sum ( int n, double a[] ) /******************************************************************************/ /* Purpose: r8vec_sum() returns the sum of an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double A[N], the vector. Output, double R8VEC_SUM, the sum of the vector. */ { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + a[i]; } return value; } /******************************************************************************/ double *r8vec_uniform_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: r8vec_uniform_01_new() returns a unit pseudorandom R8VEC. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the MIT license. Modified: 19 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Second Edition, Springer, 1987, ISBN: 0387964673, LC: QA76.9.C65.B73. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, December 1986, pages 362-376. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation, edited by Jerry Banks, Wiley, 1998, ISBN: 0471134031, LC: T57.62.H37. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, Number 2, 1969, pages 136-143. Parameters: Input, int N, the number of entries in the vector. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values. */ { int i; int i4_huge = 2147483647; int k; double *r; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8VEC_UNIFORM_01_NEW - Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r[i] = ( double ) ( *seed ) * 4.656612875E-10; } return r; }