# include <math.h>
# include <stdlib.h>
# include <stdio.h>
# include <time.h>
# include "disk01_integrals.h"
/******************************************************************************/
double disk01_area ( )
/******************************************************************************/
/*
Purpose:
DISK01_AREA returns the area of the unit disk in 2D.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
03 January 2014
Author:
John Burkardt
Parameters:
Output, double DISK01_AREA, the area of the unit disk.
*/
{
double area;
const double r = 1.0;
const double r8_pi = 3.141592653589793;
area = r8_pi * r * r;
return area;
}
/******************************************************************************/
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;
}
/******************************************************************************/
double *disk01_sample ( int n, int *seed )
/******************************************************************************/
/*
Purpose:
DISK01_SAMPLE uniformly samples the unit disk in 2D.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
03 January 2014
Author:
John Burkardt
Parameters:
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 norm;
double r;
double *v;
double *x;
x = ( double * ) malloc ( 2 * n * sizeof ( double ) );
for ( j = 0; j < n; j++ )
{
v = r8vec_normal_01_new ( 2, seed );
/*
Compute the length of the vector.
*/
norm = sqrt ( pow ( v[0], 2 ) + pow ( v[1], 2 ) );
/*
Normalize the vector.
*/
for ( i = 0; i < 2; i++ )
{
v[i] = v[i] / norm;
}
/*
Now compute a value to map the point ON the circle INTO the circle.
*/
r = r8_uniform_01 ( seed );
for ( i = 0; i < 2; i++ )
{
x[i+j*2] = sqrt ( r ) * v[i];
}
free ( v );
}
return x;
}
/******************************************************************************/
int *i4vec_uniform_ab_new ( int n, int a, int b, int *seed )
/******************************************************************************/
/*
Purpose:
I4VEC_UNIFORM_AB_NEW returns a scaled pseudorandom I4VEC.
Discussion:
The pseudorandom numbers should be uniformly distributed
between A and B.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
24 May 2012
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, integer N, the dimension of the vector.
Input, int A, B, the limits of the interval.
Input/output, int *SEED, the "seed" value, which should NOT be 0.
On output, SEED has been updated.
Output, int I4VEC_UNIFORM_AB_NEW[N], a vector of random values
between A and B.
*/
{
int c;
int i;
int i4_huge = 2147483647;
int k;
float r;
int value;
int *x;
if ( *seed == 0 )
{
fprintf ( stderr, "\n" );
fprintf ( stderr, "I4VEC_UNIFORM_AB_NEW - Fatal error!\n" );
fprintf ( stderr, " Input value of SEED = 0.\n" );
exit ( 1 );
}
/*
Guaranteee A <= B.
*/
if ( b < a )
{
c = a;
a = b;
b = c;
}
x = ( int * ) malloc ( n * sizeof ( int ) );
for ( i = 0; i < n; i++ )
{
k = *seed / 127773;
*seed = 16807 * ( *seed - k * 127773 ) - k * 2836;
if ( *seed < 0 )
{
*seed = *seed + i4_huge;
}
r = ( float ) ( *seed ) * 4.656612875E-10;
/*
Scale R to lie between A-0.5 and B+0.5.
*/
r = ( 1.0 - r ) * ( ( float ) a - 0.5 )
+ r * ( ( float ) b + 0.5 );
/*
Use rounding to convert R to an integer between A and B.
*/
value = round ( r );
/*
Guarantee A <= VALUE <= B.
*/
if ( value < a )
{
value = a;
}
if ( b < value )
{
value = b;
}
x[i] = value;
}
return x;
}
/******************************************************************************/
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 GNU LGPL 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_un