/* nag_quad_md_sphere (d01fdc) Example Program. * * Copyright 2011, Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #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); printf("nag_quad_md_sphere (d01fdc) Example Program Results\n"); /* For communication with user-supplied functions: */ comm.user = ruser; /* Skip heading in data file */ scanf("%*[^\n] "); scanf("%ld", &ndim); scanf("%ld", &limit); scanf("%lf %*[^\n] ", &u); for (i = 1; i <= 2; i++) { /* nag_quad_md_sphere (d01fdc). * Multidimensional quadrature, Sag-Szekeres method, * 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) { printf("Error from nag_quad_md_sphere (d01fdc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\nEstimated value of the integral = %9.3f" "\nNumber of integrand evaluations = %4ld\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; } }