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
dist2.c
# include <math.h>
#define YES (0 == 0)
#define NO (!YES)
/*
dist2: squared distance in torus
dist2thresh: faster code for testing whether dist2 < r2
*/
double dist2(u,v,x,y,period)
double u, v, x, y, *period;
{
double wide, high, a1, a2, a3, b1, b2, b3, amin, bmin, d2;
wide = period[0];
high = period[1];
/*
The following line is commented out,
as we now avoid calling 'dist2' in this case:
if(wide < 0.) return pow(u-x,2) + pow(v-y,2);
*/
a1 = pow(u-x,2);
a2 = pow(u-x+wide,2);
a3 = pow(u-x-wide,2);
b1 = pow(v-y,2);
b2 = pow(v-y+high,2);
b3 = pow(v-y-high,2);
amin = a1;
if(a2 < amin) amin=a2;
if(a3 < amin) amin=a3;
bmin = b1;
if(b2 < bmin) bmin=b2;
if(b3 < bmin) bmin=b3;
d2 = amin + bmin;
return d2;
}
double dist2either(u,v,x,y,period)
double u, v, x, y, *period;
{
if(period[0] < 0.) return pow(u-x,2) + pow(v-y,2);
return(dist2(u,v,x,y,period));
}
int dist2thresh(u,v,x,y,period,r2)
double u, v, x, y, *period, r2;
{
double wide, high, amin, bmin, aa, bb, residue, dx, dy;
wide = period[0];
high = period[1];
dx = u - x;
amin = dx * dx;
aa = (dx - wide) * (dx - wide);
if(aa < amin) amin = aa;
aa = (dx + wide) * (dx + wide);
if(aa < amin) amin = aa;
residue = r2 - amin;
if(residue < 0)
return NO;
dy = v - y;
bmin = dy * dy;
bb = (dy - high) * (dy - high);
if(bb < bmin) bmin= bb;
bb = (dy + high) * (dy + high);
if(bb < bmin) bmin= bb;
residue = residue - bmin;
return (residue >= 0);
}
Computing file changes ...