I have found two bugs in the code in SurfaceAreaAndVolume.cpp. These are simplified from a real case.
Bug 1: zero area&volume when no intersections at all
For a single sphere at (0, 0, 0) radius 2
it computes area=0 (expected 4 pi r*r)
and vol=0 (expected 4/3 pi r*r*r)
Bug 2: incorrect areas when little spheres intersect inside bigger sphere
For three spheres:
center ( 0, 0, 1.0) radius 2
center ( delta, 2*delta, 0) radius 1
center (-delta, 2*delta, 0) radius 1
Let delta take on the increasing values in the second column below. The SASAs for the second and third spheres are pretty wacko, until the spheres move far enough apart that are no longer intersecting. The SASA values appear in the third and fourth columns below.
n delta sasa1 sasa2
0 0.0250 0.01713 -6.11051
1 0.1250 0.40944 -5.12501
2 0.2250 1.20920 -3.76828
3 0.3250 2.23546 -2.19990
4 0.4250 3.35080 -0.53164
5 0.5250 4.48483 1.19139
6 0.6250 5.61398 2.96453
7 0.7250 6.74001 4.79053
8 0.8250 7.88015 6.65789
9 0.9250 9.07469 9.07469
10 1.0250 10.37016 10.37016
11 1.1250 11.30607 11.30607
12 1.2250 12.21203 12.21203
13 1.3250 0.00000 0.00000
14 1.4250 0.00000 0.00000
Problems:
- The areas shouldn't be negative.
- They should be equal to each other.
- The zeros in rows 13-14 are presumably the same problem as the single-sphere Bug #1 above.
The above SASA values were printed by the code below. It was pasted into testSurfaceAreaAndVolume.cpp at the end of main().
fprintf(stdout, "\n\nTEST6\n\n");
av.clear();
radii.clear();
double delta = 0.025;
a.setCoor( 0, 0, 1.0); av.push_back(&a); radii.push_back(2);
b.setCoor( delta, 2*delta, 0); av.push_back(&b); radii.push_back(1);
c.setCoor(-delta, 2*delta, 0); av.push_back(&c); radii.push_back(1);
for (int step = 0; step < 15; step++) {
b.setCoor( delta, 2*delta, 0);
c.setCoor(-delta, 2*delta, 0);
sav.computeSurfaceAreaAndVolume(av, radii);
double sa1 = sav.getRadiiSurfaceAreaAndVolume(&b)[1];
double sa2 = sav.getRadiiSurfaceAreaAndVolume(&c)[1];
fprintf(stdout, "%3d %.4f %.5f %.5f\n", step, delta, sa1, sa2);
delta += 0.1;
}