/*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_dintmr s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_dintmr.h"
/*     program DRDINTMR
 *>> 2003-03-30 DRDINTMR Krogh Changed way of dealing with 0 denominator.
 *>> 2001-05-22 DRDINTMR Krogh Minor change for making .f90 version.
 *>> 1994-10-19 DRDINTMR Krogh  Changes to use M77CON
 *>> 1994-08-31 DRDINTMR Snyder  Moved formats for C conversion
 *>> 1994-08-08 DRDINTMR Snyder  Took '0' out of formats for C conversion
 *>> 1991-11-20 CLL Edited for Fortran 90.
 *>> 1987-12-09 DRDINTMR Snyder  Initial Code.
 *--D replaces "?": DR?INTMR, ?INTM, ?INTMA, ?INTA
 *
 *     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 reverse
 *     communication.
 *
 *     This sample 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.
 * */
		/* PARAMETER translations */
#define	NDIMI	2
#define	NWORK	(3*NDIMI + NWPD*(NDIMI - 1))
#define	NWPD	217
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int iflag, iopt[10];
	double answer, cosy, denom, work[NWORK], xdy, ysq;
		/* 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 DRDINTMR:\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] = 6;
	Iopt[5] = 0;
	dintm( NDIMI, &answer, work, NWORK, iopt );
L_30:
	;
	dintma( &answer, work, iopt );
	iflag = Iopt[1];
	if (iflag == 0)
	{
 
		/*     IFLAG = 0, compute innermost integrand.
		 *
		 *     Iterate on the innermost integrand by calling DINTA here.
		 *     The code would be slightly simpler, but slightly slower,
		 *     if control of the iteration were relegated to DINTMA.
		 * */
L_40:
		;
		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)));
		}
		dinta( &answer, work, iopt );
		if (Iopt[1] == 0)
			goto L_40;
 
		if (Iopt[1] > 0)
		{
			/*                                 Done with the integration */
			Iopt[1] = -(Iopt[1] + NDIMI);
			goto L_60;
		}
 
	}
	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
	{
		if (iflag + NDIMI <= 0)
			goto L_60;
 
		/*     IFLAG < 0, transform inner integrand.
		 * */
		answer *= cosy;
		Work[1] *= cosy;
 
	}
	goto L_30;
 
L_60:
	;
 
	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 */