/* Monte Carlo integration for functions in R^3.
 * 
 * Note: This code only samples above z = 0. If you wish to do otherwise, 
 * please modify the function mc_int().
 *
 * Implemented in ANSI C.
 *
 * Author: Michael Mazack, <mazack @ yahoo . com>
 *
 * Date: February 3rd, 2010
 *
 * License: Public Domain. Redistribution and modification without 
 * restriction is granted. 
 */
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

/* Prototype for the function to integrate. */
double func(double, double);

/* Prototype for the domain function. */
int    domain(double, double);

/* Prototype for Monte Carlo Integration. */
double mc_int(double (*)(double, double), int (*)(double, double),
	      unsigned long, double, double, double);

int main()
{
  double vol;
  double pi;

  /* Pi. */
  pi = acos(-1.0);

  /* Seed the random number generator. */
  srand(time(NULL));

  /* Get the Monte Carlo volume. */
  vol = mc_int(&func, &domain, 10000, -1, -1, 1);

  /* Display results. */
  printf("\nMonte Carlo integration using 10000 samples.\n");
  printf("Approximate volume of a hemisphere    : %lf\n", vol);
  printf("Actual volume of a hemisphere         : %lf\n", 0.5*(4.0/3.0)*pi);
  printf("\n");

  return 0;
}

/* Note: Any function in R^3 can be used here. Which coordinate system
 * you use is up to you. However, it's better to sample in cartesian
 * coordinates if you are using a uniform "random" distribution.
 */
double func(double x, double y)
{
  /* Top half of unit sphere. */
  return sqrt(1 - x*x - y*y);
}

/* The domain function returns 1 (true) if the point is in the domain,
 * otherwise it returns 0 (false). This is expected by mc_int().
 */
int domain(double x, double y)
{
  /* Unit disc. */
  if(x*x + y*y <= 1)
    return 1;
  else
    return 0;
}

/* Implementation of Monte Carlo integration for functions in R^3.
 *
 * First argument: Function to integrate.
 * Second argument: Domain to integrate over.
 * Third argument:  Number of samples.
 * Fourth argument: Magnitude of x boundary to sample. (1/2 length of cube).
 * Fifth argument: Magnitude of y boundary to sample. (1/2 width of cube).
 * Sixth argument: Magnitude of z boundary to sample. (height of cube from 0).
 */
double mc_int(double (*f)(double, double), int (*d)(double, double),
	      unsigned long samps, double xbound, double ybound, double zbound)
{
  unsigned long i, inside;
  double x, y, z;

  /* Initialize. */
  inside = 0;

  /* Sample as many times as requested. */
  for(i = 0; i < samps; i++)
    {
      /* Generate random x and y inside the cube. */
      x = -xbound + 2.0*xbound*(rand()/(RAND_MAX + 1.0));
      y = -ybound + 2.0*ybound*(rand()/(RAND_MAX + 1.0));

      /* See if the sample is inside the domain. */
      if(domain(x, y))
	{
	  /* Generate random z inside the cube. */
	  z = zbound*(rand()/(RAND_MAX + 1.0));

	  /* Check to see if z lies below the function and increment if so. */
	  if(z < f(x, y) && z > 0)
	    inside++;
	}
    }

  /* Return the probability of sampling inside the function times the volume
   * of the box (i.e. the volume of what we are integrating).
   */
  return ((double)inside/(double)samps)*4.0*xbound*ybound*zbound;
}