/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:10 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_dintmf s=dbov str=l x=f - prototypes */
#include
#include "fcrt.h"
#include
#include
#include "p_dintmf.h"
/* DRDINTM
*>> 2003-03-30 DRDINTMF Krogh Changed way of dealing with 0 denominator.
*>> 2001-05-22 DRDINTMF Krogh Minor change for making .f90 version.
*>> 1994-11-02 DRDINTMF Krogh Changes to use M77CON
*>> 1994-08-08 DRDINTMF Snyder Took '0' out of formats for C conversion
*>> 1993-05-05 DRDINTMF Krogh Adjusted to simplify conversion to C.
*>> 1987-12-09 DRDINTMF Snyder Initial Code.
*--D replaces "?": DR?INTM, DR?INTMF, ?INTM, ?INTF
*
* DEMO DRIVER for multi-dimensional quadrature subprogram DINTM.
*
* Compute the integral for Y = 0.0 to Y = PI of
* the integral for X = 0.0 to X = Y of
* X * COS(Y) /(X*X+Y*Y). The ANSWER should be zero.
*
* The integrand and all limits are provided by the subprogram DINTF,
* although the limits of the outer dimension, and the lower limit of
* the inner dimension, being constants, could have been provided
* here.
* */
/* PARAMETER translations */
#define NDIMI 2
#define NWORK (3*NDIMI + NWPD*(NDIMI - 1))
#define NWPD 217
/* end of PARAMETER translations */
int main( )
{
long int iopt[10];
double answer, work[NWORK];
/* OFFSET Vectors w/subscript range: 1 to dimension */
long *const Iopt = &iopt[0] - 1;
double *const Work = &work[0] - 1;
/* end of OFFSET VECTORS */
printf("\n DRDINTMF:\n Compute the integral for Y = 0.0 to Y = PI of\n the integral for X = 0.0 to X = Y of\n X * COS(Y) /(X*X+Y*Y). The ANSWER should be zero.\n");
Iopt[2] = 10;
Iopt[3] = 0;
Iopt[4] = 0;
dintm( NDIMI, &answer, work, NWORK, iopt );
printf("\n ANSWER = %15.8g\n ERROR ESTIMATE =%15.8g\n STATUS FLAG = %3ld\n FUNCTION VALUES (INNERMOST INTEGRALS) =%6ld\n",
answer, Work[1], Iopt[1], Iopt[3]);
exit(0);
} /* end of function */
/* End of DRDINTMF */
void /*FUNCTION*/ dintf(
double *answer,
double work[],
long *iflag)
{
double denom, xdy;
static double cosy, ysq;
static long ndimi = 2;
/* OFFSET Vectors w/subscript range: 1 to dimension */
double *const Work = &work[0] - 1;
/* end of OFFSET VECTORS */
/* Subroutine to provide integrand and limits for DINTM.
*
* This sample subroutine demonstrates the use of functional
* transformation of the inner integral to reduce the cost of the
* overall computation. The factor "cos y" does not depend on the
* variable of integration for the inner integral, and therefore is
* factored out. Similarly, the term y*y in the denominator of the
* inner integrand is pre-computed when the limits of the inner
* integral are requested.
*
* This example also demonstrates the necessity to cope with overflow
* of the integrand, even though the integral is well behaved.
* */
if (*iflag == 0)
{
/* IFLAG = 0, compute innermost integrand.
* */
denom = Work[1]*Work[1] + ysq;
if (denom != 0.0e0)
{
/* This test will not detect overflow of the integrand if the
* arithmetic underflows gradually. */
*answer = Work[1]/denom;
}
else
{
/* Special care to avoid a zero denominator. */
xdy = Work[1]/Work[2];
*answer = xdy/(Work[2]*(1.e0 + SQ(xdy)));
}
}
else if (*iflag == 1)
{
/* Compute limits of inner dimension.
*
* Set WORK(1) = COS(WORK(2)) = Partial derivative, with respect to
* the integral over the inner dimension, of the transformation
* applied when integration over the inner dimension is complete.
* This is so dintm knows how much accuracy it needs for this integral.
* */
Work[ndimi + *iflag] = 0.0e0;
Work[2*ndimi + *iflag] = Work[2];
ysq = Work[2]*Work[2];
cosy = cos( Work[2] );
Work[1] = cosy;
}
else if (*iflag == 2)
{
/* Compute limits of outer dimension.
* */
Work[ndimi + *iflag] = 0.0e0;
Work[2*ndimi + *iflag] = 4.0e0*atan( 1.0e0 );
}
else
{
/* IFLAG < 0, transform inner integrand.
* */
*answer *= cosy;
Work[1] *= cosy;
}
return;
} /* end of function */