Revision d606122dc24b56ecf537d55eda38f4bf5ac4de1f authored by Adrian Baddeley on 25 October 2010, 10:40:51 UTC, committed by cran-robot on 25 October 2010, 10:40:51 UTC
1 parent 66bc933
sphevol.c
#include <math.h>
#include <R.h>
#include "geom3.h"
/*
$Revision: 1.1 $ $Date: 2009/11/04 23:54:15 $
Routine for calculating ABSOLUTE volume of intersection
between sphere and box
Arbitrary positions: point is allowed to be inside or outside 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
#define TRUE (0 == 0)
#define FALSE (!TRUE)
#define ABS(X) ((X >= 0.0) ? (X) : -(X))
static double rcubed, spherevol;
double
sphevol(point, box, r)
Point *point;
Box *box;
double r;
{
double sum, p[4], q[4];
double v1(), v2(), v3();
int i, j;
rcubed = r * r * r;
spherevol = (4.0/3.0) * PI * rcubed;
p[1] = box->x0 - point->x;
p[2] = box->y0 - point->y;
p[3] = box->z0 - point->z;
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 += v1(p[i], -1, r) + v1(q[i], 1, r);
#ifdef DEBUG
Rprintf("i = %d, v1 = %f, v1 = %f\n", i, v1(p[i], -1, r), v1(q[i], 1, r));
#endif
}
DBG("Past v1", sum)
for(i = 1; i < 3; i++)
for(j = i+1; j <= 3; j++)
{
sum -= v2(p[i], -1, p[j], -1, r) + v2(p[i], -1, q[j], 1, r)
+ v2(q[i], 1, p[j], -1, r) + v2(q[i], 1, q[j], 1, r);
#ifdef DEBUG
Rprintf("i = %d, j = %d, sum = %f\n", i, j, sum);
#endif
}
DBG("Past v2", sum)
sum += v3(p[1], -1, p[2], -1, p[3], -1, r)
+ v3(p[1], -1, p[2], -1, q[3], 1, r);
DBG("sum", sum)
sum += v3(p[1], -1, q[2], 1, p[3], -1, r)
+ v3(p[1], -1, q[2], 1, q[3], 1, r);
DBG("sum", sum)
sum += v3(q[1], 1, p[2], -1, p[3], -1, r)
+ v3(q[1], 1, p[2], -1, q[3], 1, r);
DBG("sum", sum)
sum += v3(q[1], 1, q[2], 1, p[3], -1, r)
+ v3(q[1], 1, q[2], 1, q[3], 1, r);
DBG("Past v3", sum)
DBG("sphere volume", spherevol)
return(spherevol - sum);
}
double
v1(a,s,r)
double a, r;
int s;
{
double value;
double u();
short sign;
value = 4.0 * rcubed * u(ABS(a)/r, 0.0, 0.0);
sign = (a >= 0.0) ? 1 : -1;
if(sign == s)
return(value);
else
return(spherevol - value);
}
double
v2(a, sa, b, sb, r)
double a, b, r;
int sa, sb;
{
short sign;
double u();
sign = (b >= 0.0) ? 1 : -1;
if(sign != sb )
return(v1(a, sa, r) - v2(a, sa, ABS(b), 1, r));
b = ABS(b);
sb = 1;
sign = (a >= 0.0) ? 1 : -1;
if(sign != sa)
return(v1(b, sb, r) - v2(ABS(a), 1, b, sb, r));
a = ABS(a);
return(2.0 * rcubed * u(a/r, b/r, 0.0));
}
double
v3(a, sa, b, sb, c, sc, r)
double a, b, c, r;
int sa, sb, sc;
{
short sign;
double u();
sign = (c >= 0.0) ? 1 : -1;
if(sign != sc)
return(v2(a,sa,b,sb,r) - v3(a,sa,b,sb, ABS(c), 1, r));
c = ABS(c);
sc = 1;
sign = (b >= 0.0) ? 1 : -1;
if(sign != sb)
return(v2(a,sa,c,sc,r) - v3(a,sa,ABS(b),1,c,sc,r));
b = ABS(b);
sb = 1;
sign = (a >= 0.0) ? 1 : -1;
if(sign != sa)
return(v2(b,sb, c, sc, r) - v3(ABS(a),1, b, sb, c, sc, r));
a = ABS(a);
return(rcubed * u(a/r, b/r, c/r));
}
double
u(a, b, c)
double a, b, c;
{
double w();
if(a * a + b * b + c * c >= 1.0)
return(0.0);
return(
(PI/12.0) * (2.0 - 3.0 * (a + b + c)
+ (a * a * a + b * b * b + c * c * c))
+ w(a,b) + w(b,c) + w(a,c)
- a * b * c
);
}
double
w(x,y)
double x,y; /* Arguments assumed >= 0 */
{
double z;
z = sqrt(1 - x * x - y * y);
return(
(x / 2.0 - x * x * x / 6.0) * atan2(y, z)
+ (y / 2.0 - y * y * y / 6.0) * atan2(x, z)
- ( atan2(x * y , z) - x * y * z )/3.0
);
}
Computing file changes ...