//=========================================================================
#define PROGRAMMER "BAHN, William"
#define PROG_CODE  "bahw"
#define COURSE     "ECE-1021"
#define YEAR       (2003)
#define TERM       "Fall"
#define SECTION    (0)
#define ASSIGNMENT "RWA 5.7"
#define REVISION   (0)
#define TITLE      "Real World Application from Section 5.7"
#define SUBTITLE   "Monte Carlo Integration"
#define EMAIL      "wbahn@eas.uccs.edu"
#define FILENAME   "rwa0507.c"
//=========================================================================

// PROBLEM:
//
// MODIFIED FROM TEXT
//
// Using Monte Carlo integration, evaluate the integral
//
// Integral{x = 0 to 1} ( e^[-(x^2)/2] )
//
// BASIC APPROACH
//
// 1) REM: Get limits and number of trials from user.
// 2) LOOP: Until the number of trails of been run
// 2.1) REM: Generate random (x,y) test pair within the set of limits.
// 2.2) REM: Determine if the test pair is bounded by the function.
// 2.3) REM: Keep count of total number of pairs within function's area.
// 3) Area = (Area of test box) * (fraction of pairs within function)
//
// NOTE: The assumption is that the function be integrated is never
//       negative within the test region.

//== INCLUDE FILES ========================================================
#include <stdio.h>  // printf(), scanf()
#include <math.h>   // exp()
#include <stdlib.h> // rand(), srand(), RAND_MAX
#include <limits.h> // LONG_MAX

//== FUNCTION PROTOTYPES ==================================================
void PrintHeader(void);
double GetMaxY(void);
double GetMinX(void);
double GetMaxX(void);
double GetTrials(void);
int    UnderCurve(double x, double y);
double RandFloat(double min, double max);
void   PrintResult(double min, double max, double area);

//== MAIN FUNCTION ========================================================
int main(void)
{
	long i, trials, undercurve;
	double ymax, xmin, xmax;
	double x, y;
	double area;

	PrintHeader();
	printf("\n"); // Print a blank line

	ymax   = GetMaxY();
	xmin   = GetMinX();
	xmax   = GetMaxX();
	trials = GetTrials();

	for(undercurve = 0, i = 0; i < trials; i++)
	{
		x = RandFloat(xmin, xmax);
		y = RandFloat(0.0, ymax);
		undercurve += UnderCurve(x,y);
	}

	area = ymax * (xmax - xmin) * ((double) undercurve / (double) trials );
	PrintResult(xmin, xmax, area);

	return(0);
}

//== SUPPORT FUNCTIONS ====================================================
void PrintHeader(void)
{
	printf("========================================"
			 "=======================================\n");
	printf("Course....... %s-%i (%s %i)\n", COURSE, SECTION, TERM, YEAR);
	printf("Programmer... %s (%s)\n", PROGRAMMER, PROG_CODE);
	printf("Assignment... %s (Rev %i) (Source Code in %s)\n",
		ASSIGNMENT, REVISION, FILENAME);
	printf("Description.. %s\n", TITLE);
	printf("              %s\n", SUBTITLE);
	printf("========================================"
			 "=======================================\n");

	return;
}

double GetMaxY(void)
{
	double v;

	printf("\n");
	printf("Enter a value that is GREATER than any value of the function\n");
	printf("within the limits of integration. The closer the value is to\n");
	printf("actual maximum the better the results will be.\n");
	printf("\n");
	printf("Limiting y-value: ");
	scanf("%lf", &v);

	return(v);
}

double GetMinX(void)
{
	double v;

	printf("\n");
	printf("Enter the LOWER integration limit: ");
	scanf("%lf", &v);

	return(v);
}

double GetMaxX(void)
{
	double v;

	printf("\n");
	printf("Enter the UPPER integration limit: ");
	scanf("%lf", &v);

	return(v);
}

double GetTrials(void)
{
	long int i;

	printf("\n");
	printf("Number of trials (1 to %li): ", LONG_MAX);
	scanf("%li", &i);

	return(i);
}

double function(double x)
{
	return( exp(-x*x/2.0));
}

int UnderCurve(double x, double y)
{
	return(y <= function(x));
}

double RandFloat(double min, double max)
{
	return( min + (max - min)*((double) rand() / (double) RAND_MAX) );
}

void PrintResult(double min, double max, double area)
{
	printf("\n");
	printf("The integral evaluated from %f to %f is %f.\n",
		min, max, area);
	return;
}