/*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 <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#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 */