```/* nag_quad_md_sphere (d01fdc) Example Program.
*
* NAGPRODCODE Version.
*
* Copyright 2016 Numerical Algorithms Group.
*
* Mark 26, 2016.
*/
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd01.h>

#ifdef __cplusplus
extern "C"
{
#endif
static double NAG_CALL f(Integer ndim, const double x[], Nag_Comm *comm);
static void NAG_CALL region(Integer ndim, const double x[], Integer j,
double *c, double *d, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{
static double ruser[2] = { -1.0, -1.0 };
Integer exit_status = 0;
double r0, result, sigma, u;
Integer i, limit, ncalls, ndim;
Nag_Comm comm;
NagError fail;

INIT_FAIL(fail);

/* For communication with user-supplied functions: */
comm.user = ruser;

/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%" NAG_IFMT "", &ndim);
scanf("%" NAG_IFMT "", &limit);
scanf("%lf %*[^\n] ", &u);

for (i = 1; i <= 2; i++) {
* general product region or n-sphere.
*/
switch (i) {
case 1:
printf("\nSphere-to-sphere transformation\n");
sigma = 1.5;
r0 = 0.9;
nag_quad_md_sphere(ndim, f, sigma, NULLFN, limit, r0, u, &result,
&ncalls, &comm, &fail);
break;
case 2:
printf("\nProduct region transformation\n");
sigma = -1.0;
r0 = 0.8;
nag_quad_md_sphere(ndim, f, sigma, region, limit, r0, u, &result,
&ncalls, &comm, &fail);
break;
}
if (fail.code != NE_NOERROR) {
exit_status = 1;
goto END;
}

printf("\nEstimated value of the integral = %9.3f"
"\nNumber of integrand evaluations = %4" NAG_IFMT "\n",
result, ncalls);
}

END:
return exit_status;
}

static double NAG_CALL f(Integer ndim, const double x[], Nag_Comm *comm)
{
Integer i;
double x_sq = 0.0;

if (comm->user[0] == -1.0) {
printf("(User-supplied callback f, first invocation.)\n");
comm->user[0] = 0.0;
}

for (i = 0; i < ndim; i++)
x_sq += pow(x[i], 2.0);

return 1.0 / sqrt(fabs(pow(1.5, 2) - x_sq));
}

static void NAG_CALL region(Integer ndim, const double x[], Integer j,
double *c, double *d, Nag_Comm *comm)
{
Integer i;
double x_sq = 0.0;

if (comm->user[1] == -1.0) {
printf("(User-supplied callback region, first invocation.)\n");
comm->user[1] = 0.0;
}
if (j > 1) {
for (i = 0; i < (j - 1); i++)
x_sq += pow(x[i], 2.0);
*d = sqrt(fabs(pow(1.5, 2) - x_sq));
*c = -*d;
}
else {
*c = -1.5;
*d = 1.5;
}
}
```