Revision c9b2c621c3bff55aaa77646dc1ba7316765cd7e4 authored by Adrian Baddeley on 25 April 2013, 00:00:00 UTC, committed by Gabor Csardi on 25 April 2013, 00:00:00 UTC
1 parent f86606a
sphefrac.c
#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);
}
Computing file changes ...