in packages/miew/src/gfx/geometries/ContactSurface.js [454:548]
function projectTorus(a, b) {
const aIdx = itemSize * a;
const bIdx = itemSize * b;
const xa = posRad[aIdx];
const ya = posRad[aIdx + 1];
const za = posRad[aIdx + 2];
const r1 = posRad[aIdx + 3];
let dx = mid.x = posRad[bIdx] - xa;
let dy = mid.y = posRad[bIdx + 1] - ya;
let dz = mid.z = posRad[bIdx + 2] - za;
const innR2 = posRad[bIdx + 3];
let d2 = dx * dx + dy * dy + dz * dz;
// This check now redundant as already done in AVHash.withinRadii
// if( d2 > (( r1 + r2 ) * ( r1 + r2 )) ){ return; }
const d = Math.sqrt(d2);
// Find angle between a->b vector and the circle
// of their intersection by cosine rule
const cosA = (r1 * r1 + d * d - innR2 * innR2) / (2.0 * r1 * d);
// distance along a->b at intersection
const dmp = r1 * cosA;
mid.normalize();
// Create normal to line
normalToLine(n1, mid);
n1.normalize();
// Cross together for second normal vector
n2.crossVectors(mid, n1);
n2.normalize();
// r is radius of circle of intersection
const rInt = Math.sqrt(r1 * r1 - dmp * dmp);
n1.multiplyScalar(rInt);
n2.multiplyScalar(rInt);
mid.multiplyScalar(dmp);
mid.x += xa;
mid.y += ya;
mid.z += za;
lastClip = -1;
const ng = ngTorus;
for (let innI = 0; innI < probePositions; innI++) {
const cost = cosTable[innI];
const sint = sinTable[innI];
const px = mid.x + cost * n1.x + sint * n2.x;
const py = mid.y + cost * n1.y + sint * n2.y;
const pz = mid.z + cost * n1.z + sint * n2.z;
if (obscured(px, py, pz, a, b) === -1) {
// As above, iterate over our grid...
// px, py, pz in grid coords
const iax = Math.floor(scaleFactor * (px - min[0]));
const iay = Math.floor(scaleFactor * (py - min[1]));
const iaz = Math.floor(scaleFactor * (pz - min[2]));
const minx = Math.max(0, iax - ng);
const miny = Math.max(0, iay - ng);
const minz = Math.max(0, iaz - ng);
const maxx = Math.min(dim[0], iax + ng + 2);
const maxy = Math.min(dim[1], iay + ng + 2);
const maxz = Math.min(dim[2], iaz + ng + 2);
for (let iz = minz; iz < maxz; iz++) {
dz = pz - gridz[iz];
const zOffset = dim[1] * dim[0] * iz;
for (let iy = miny; iy < maxy; iy++) {
dy = py - gridy[iy];
const dzy2 = dz * dz + dy * dy;
const zyOffset = zOffset + dim[0] * iy;
for (let ix = minx; ix < maxx; ix++) {
dx = px - gridx[ix];
d2 = dzy2 + dx * dx;
const idx = ix + zyOffset;
const current = grid[idx];
if (current > 0.0 && d2 < (current * current)) {
grid[idx] = Math.sqrt(d2);
}
}
}
}
}
}
}