# include # include # include # include # include "annulus_rule.h" /******************************************************************************/ void annulus_rule_compute ( double center[2], double r1, double r2, int nr, int nt, double w[], double x[], double y[] ) /******************************************************************************/ /* Purpose: ANNULUS_RULE_COMPUTE computes a quadrature rule for an annulus. Discussion: The integration region is points (X,Y) such that R1^2 <= ( X - CENTER(1) )^2 + ( Y - CENTER(2) )^2 <= R2^2 Licensing: This code is distributed under the GNU LGPL license. Modified: 06 July 2018 Author: John Burkardt Reference: William Peirce, Numerical Integration Over the Planar Annulus, Journal of the Society for Industrial and Applied Mathematics, Volume 5, Issue 2, June 1957, pages 66-73. Parameters: Input, double CENTER(2), the coordinates of the center. Input, double R1, R2, the inner and outer radii of the annulus. 0.0 <= R1 <= R2. Input, int NR, the number of points in the radial rule. Input, int NT, the number of angles to use. The value NT=4*NR is recommended. Output, double W[NR*NT], the weights for the rule. Output, double X[NR*NT], Y[NR*NT], the points for the rule. */ { double a; double area; double b; double c; double d; int i; int j; int k; const double r8_pi = 3.141592653589793; double *ra; double *rw; double t; double tw; /* Request a Legendre rule for [-1,+1]. */ ra = ( double * ) malloc ( nr * sizeof ( double ) ); rw = ( double * ) malloc ( nr * sizeof ( double ) ); legendre_ek_compute ( nr, ra, rw ); /* Adjust the rule from [-1,+1] to [r1^2,r2^2]. */ a = -1.0; b = +1.0; c = r1 * r1; d = r2 * r2; rule_adjust ( a, b, c, d, nr, ra, rw ); /* Convert from R^2 to R. */ for ( i = 0; i < nr; i++ ) { ra[i] = sqrt ( ra[i] ); } for ( i = 0; i < nr; i++ ) { rw[i] = rw[i] / ( r2 + r1 ) / ( r2 - r1 ); } /* Set the angular weight. */ tw = 1.0 / ( double ) ( nt ); /* Get area of annulus. */ area = annulus_area ( center, r1, r2 ); /* Form the abscissa and weight vectors. */ k = 0; for ( i = 0; i < nt; i++ ) { t = 2.0 * r8_pi * ( double ) ( i ) / ( double ) ( nt ); for ( j = 0; j < nr; j++ ) { x[k] = center[0] + ra[j] * cos ( t ); y[k] = center[1] + ra[j] * sin ( t ); w[k] = area * tw * rw[j]; k = k + 1; } } free ( ra ); free ( rw ); return; }