#include <math.h>
#include <R.h>
#include "geom3.h"
/*

	$Revision: 1.1 $ 	$Date: 2009/11/04 23:54:15 $

	Routine for calculating surface area of sphere
	intersected with box

# /////////////////////////////////////////////
# AUTHOR: Adrian Baddeley, CWI, Amsterdam, 1991.
# 
# MODIFIED BY: Adrian Baddeley, Perth 2009
#
# This software is distributed free
# under the conditions that
# 	(1) it shall not be incorporated
# 	in software that is subsequently sold
# 	(2) the authorship of the software shall
# 	be acknowledged in any publication that 
# 	uses results generated by the software
# 	(3) this notice shall remain in place
# 	in each file.
# //////////////////////////////////////////////

*/

#ifdef DEBUG
#define DBG(X,Y) Rprintf("%s: %f\n", (X), (Y)) 
#else
#define DBG(X,Y) 
#endif

static double pi = 3.141592653589793;

/* Factor of 4 * pi * r * r IS ALREADY TAKEN OUT */

double
sphesfrac(point, box, r)
Point *point;
Box *box;
double r;
{
	double sum, p[4], q[4];
	double a1(), a2(), a3();

	int i, j;

	p[1] = point->x - box->x0;
	p[2] = point->y - box->y0;
	p[3] = point->z - box->z0;

	q[1] = box->x1 - point->x;
	q[2] = box->y1 - point->y;
	q[3] = box->z1 - point->z;

	sum = 0;
	for(i = 1; i <= 3; i++)
	{
		sum += a1(p[i],r) + a1(q[i],r);
#ifdef DEBUG
		Rprintf("i = %d, a1 = %f, a1 = %f\n", i, a1(p[i],r), a1(q[i],r));
#endif
	}
	DBG("Past a1", sum)
	
	for(i = 1; i < 3; i++)
		for(j = i+1; j <= 3; j++)
		{
			sum -= a2(p[i], p[j], r) + a2(p[i], q[j], r)
			 + a2(q[i], p[j], r) + a2(q[i], q[j], r);
#ifdef DEBUG
			Rprintf("i = %d, j = %d, sum = %f\n", i, j, sum);
#endif
		}
	DBG("Past a2", sum)

	sum += a3(p[1], p[2], p[3], r) + a3(p[1], p[2], q[3], r);

	DBG("sum", sum)

	sum += a3(p[1], q[2], p[3], r) + a3(p[1], q[2], q[3], r);

	DBG("sum", sum)

	sum += a3(q[1], p[2], p[3], r) + a3(q[1], p[2], q[3], r);

	DBG("sum", sum)

	sum += a3(q[1], q[2], p[3], r) + a3(q[1], q[2], q[3], r);

	DBG("Past a3", sum)

	return(1 - sum);
}

double 
a1(t, r)
double t, r;
{
	/* This is the function A1 divided by 4 pi r^2 */

	if(t >= r) 
		return(0.0);

	return((1 - t/r) * 0.5);
}

double
a2(t1, t2, r)
double t1, t2, r;
{
	double c2();
	/* This is A2 divided by 4 pi r^2 because c2 is C divided by pi */
	
	return(c2( t1 / r, t2 / r) / 2.0);
}

double
a3(t1, t2, t3, r)
double t1, t2, t3, r;
{
	double c3();
	/* This is A3 divided by 4 pi r^2 because c3 is C divided by pi */

	return(c3(t1 / r, t2 / r, t3 / r) / 4.0);
}

double 
c2(a, b)
double a, b;
{
	double z, z2;
	double c2();

	/* This is the function C(a, b, 0) divided by pi 
		- assumes a, b > 0  */

	if( ( z2 = 1.0 - a * a - b * b) < 0.0 )
		return(0.0);
	z = sqrt(z2);
	return((atan2(z, a * b) - a * atan2(z, b) - b * atan2(z, a)) / pi);
}

double
c3(a, b, c)
double a, b, c;
{
	double za, zb, zc, sum;
	/* This is C(a,b,c) divided by pi. Arguments assumed > 0 */

	if(a * a + b * b + c * c >= 1.0)
		return(0.0);

	za = sqrt(1 - b * b - c * c);
	zb = sqrt(1 - a * a - c * c);
	zc = sqrt(1 - a * a - b * b);

	sum =  atan2(zb, a * c) + atan2(za, b * c) 
		+ atan2(zc, a * b)
		- a * atan2(zb, c)
		+ a * atan2(b, zc)
		- b * atan2(za, c)
		+ b * atan2(a, zc)
		- c * atan2(zb, a)
		+ c * atan2(b, za);

	return(sum / pi  - 1);
}