Resources Numerical accuracy Calculation policy
Works out how much of the sky each window sees and the daylight factor on a grid of points in a room.
Triangles are 64-bit, in meters. What hides the sky is found by cutting the triangles in vertical sections; the sky is then summed over azimuth bands 2° wide with closed-form integrals in altitude. The daylight factor at a point sums patches of the window opening by the midpoint rule, so its accuracy is set by the band width and the patch count.
19 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Tests tools/bimcraft-test/ops-vsc.mjs, tools/bimcraft-test/ops-grid.mjs
The horizontal illuminance of an unobstructed overcast sky, which every component is divided by.
E_h = 7π ÷ 9, in units of the zenith luminance
Evaluated once in 64-bit from Math.PI.
export const EH = 7 * Math.PI / 9; // horizontal illuminance of the whole CIE overcast sky, in LzThe light a vertical plane gets from the sky between the horizon and an altitude.
F(θ) = (θ ÷ 2 + sin 2θ ÷ 4) ÷ 3 + (2 ÷ 9)(1 − cos³θ)
The closed form of the integral, one 64-bit expression. For a small altitude 1 − cos³θ subtracts nearly equal numbers; the term is itself small there. Callers take differences F(b) − F(a), which also lose digits for a very narrow interval.
/** ∫₀^θ L(t)/Lz cos²t dt for the CIE overcast sky (θ radians). */
export const cieF = (th) => (th / 2 + Math.sin(2 * th) / 4) / 3 + (2 / 9) * (1 - Math.cos(th) ** 3);The sky's brightness at an altitude, relative to the zenith.
L(θ) = (1 + 2 sin θ) ÷ 3
One 64-bit expression with exact whole-number factors.
export const L = (th) => (1 + 2 * Math.sin(th)) / 3;Joins the altitude ranges hidden by separate objects into one list.
Intervals in degrees are clipped to 0…90, sorted and joined when they overlap or touch within 1e-9°, an absolute tolerance.
/** Merge intervals (degrees), clipped to [lo, hi]. */
export function merge(list, lo = 0, hi = 90) {
const a = list.map(([x, y]) => [Math.max(lo, Math.min(x, y)), Math.min(hi, Math.max(x, y))]).filter(([x, y]) => y > x).sort((p, q) => p[0] - q[0]);
const out = [];
for (const iv of a) { const last = out[out.length - 1]; if (last && iv[0] <= last[1] + 1e-9) last[1] = Math.max(last[1], iv[1]); else out.push([...iv]); }
return out;
}The altitudes hidden along one level direction from a window.
A product is cut only if its plan box straddles the section line and reaches ahead of the window, tested with EPS = 1e-6 m. The hidden altitudes themselves come from daylight-bim.js sectionTheta, which cuts each triangle.
/**
* The altitudes hidden along the level direction d from origin (merged): only products whose plan box the vertical
* half-plane along d crosses are cut (a big model has thousands; one section crosses a few). A bare triangle array is cut whole.
*/
export function cut(origin, d, list) {
const s = [d[1], -d[0]], all = [];
for (const o of list) {
if (o.box) {
const b = o.box;
let neg = false, pos = false, ahead = false;
for (const x of [b.lo[0], b.hi[0]]) for (const y of [b.lo[1], b.hi[1]]) {
const side = (x - origin[0]) * s[0] + (y - origin[1]) * s[1];
if (side < EPS) neg = true;
if (side > -EPS) pos = true;
if ((x - origin[0]) * d[0] + (y - origin[1]) * d[1] > EPS) ahead = true;
}
if (!(neg && pos) || !ahead) continue;
}
all.push(...sectionTheta(origin, d, o.tr || o).blocked);
}
return merge(all);
}For each direction across the window's front, the altitudes that are hidden.
φᵢ = −90° + step ÷ 2 + i × step, i = 0 … 180 ÷ step − 1
Band centers are computed from the index, not accumulated. One vertical section at the center stands for each band; one more, square to the window, gives the angle of obstruction.
export function profile(win, tris, { step = 2 } = {}) {
const k = Math.round(180 / step), phi = [], blocked = [];
for (let i = 0; i < k; i++) {
const f = -90 + step / 2 + i * step;
phi.push(f); blocked.push(cut(win.origin, turn(win.n, win.s, f), tris));
}
// and the section square to the window (φ = 0), exact, for the angle of obstruction
return { step, phi, blocked, square: cut(win.origin, win.n, tris) };
}How much a band of directions counts, by how squarely it faces the window.
w = sin(φ + step ÷ 2) − sin(φ − step ÷ 2)
The exact integral of cos φ over the band, as a difference of two sines. Near ±90° the two sines agree to three digits, so about three digits go in the subtraction; the weight is smallest there.
const bandWeight = (f, step) => Math.sin((f + step / 2) * D2R) - Math.sin((f - step / 2) * D2R); // ∫cos φ dφ over the bandThe share of the unobstructed sky's light that reaches the window, in percent.
VSC = 100 × Σ over bands of w × Σ over visible intervals of (F(b) − F(a)) ÷ (7π ÷ 9)
A double sum in 64-bit, added in band order. The altitude integral is exact for each visible interval; the azimuth is sampled once per band. With nothing hidden the sum is 2(π ÷ 12 + 2 ÷ 9) ÷ (7π ÷ 9) = 39.62 %.
export function vsc(prof) {
let e = 0;
prof.phi.forEach((f, i) => { let s = 0; for (const [a, b] of visible(prof.blocked[i])) s += cieF(b * D2R) - cieF(a * D2R); e += bandWeight(f, prof.step) * s; });
return (e / EH) * 100;
}The share of a uniform sky's light a window still gets, used for diffuse sun gain.
fraction = Σ w × Σ (G(b) − G(a)) ÷ Σ w × G(π ÷ 2), G(θ) = θ ÷ 2 + sin 2θ ÷ 4
The same double sum for a uniform sky, divided by the unobstructed sum built from the same band weights, so rounding in the weights largely cancels in the ratio.
export function skyFraction(prof) {
let e = 0, full = 0;
prof.phi.forEach((f, i) => { const w = bandWeight(f, prof.step); let s = 0; for (const [a, b] of visible(prof.blocked[i])) s += isoF(b * D2R) - isoF(a * D2R); e += w * s; full += w * isoF(Math.PI / 2); });
return full > 0 ? e / full : 0;
}The altitude of the top of whatever hides the horizon square to the window.
The top of the first hidden interval that starts at or below 1e-9°, in the section square to the window; 0 when the horizon is clear there.
/** The top of whatever hides the horizon square to the window (0 when the horizon is clear there). */
export function obstructionAngle(prof) {
const first = (prof.square || []).find(([a]) => a <= 1e-9);
return first ? first[1] : 0;
}Decides whether a ray from a point through the window meets sky or an obstruction.
The band is found by Math.floor((φ + 90) ÷ step), clamped to the first and last band, and the altitude is tested against that band's hidden intervals with inclusive ends. A ray takes its band's section; nothing is interpolated between bands.
const blockedAt = (prof, phiDeg, altDeg) => {
const i = Math.max(0, Math.min(prof.phi.length - 1, Math.floor((phiDeg + 90) / prof.step)));
return prof.blocked[i].some(([a, b]) => altDeg >= a && altDeg <= b);
};Picks the floor triangles out of a room's shell.
floor when |n_z| > 0.7 |n| and the triangle's mean height is below the shell's middle
The normal is a 64-bit cross product of two edges. A triangle is floor when its normal is within about 46° of vertical and it lies in the lower half of the shell.
// ---------------------------------------------------------------- the grid
function floorOf(tr) {
let zlo = Infinity, zhi = -Infinity;
for (let i = 2; i < tr.length; i += 3) { zlo = Math.min(zlo, tr[i]); zhi = Math.max(zhi, tr[i]); }
const zMid = (zlo + zhi) / 2, tris = [];
for (let i = 0; i + 8 < tr.length; i += 9) {
const ax = tr[i + 3] - tr[i], ay = tr[i + 4] - tr[i + 1], az = tr[i + 5] - tr[i + 2], bx = tr[i + 6] - tr[i], by = tr[i + 7] - tr[i + 1], bz = tr[i + 8] - tr[i + 2];
const cz = ax * by - ay * bx, l = Math.hypot(ay * bz - az * by, az * bx - ax * bz, cz);
if (l > 0 && Math.abs(cz / l) > 0.7 && (tr[i + 2] + tr[i + 5] + tr[i + 8]) / 3 < zMid) tris.push([[tr[i], tr[i + 1]], [tr[i + 3], tr[i + 4]], [tr[i + 6], tr[i + 7]], (tr[i + 2] + tr[i + 5] + tr[i + 8]) / 3]);
}
return { tris, zlo };
}Tests a plan point against one floor triangle.
Barycentric coordinates in 64-bit. A triangle whose determinant is under 1e-12 in absolute value is skipped; a point counts as inside up to 1e-9 outside an edge in barycentric terms.
const inTri = (p, [a, b, c]) => { const d = (b[1] - c[1]) * (a[0] - c[0]) + (c[0] - b[0]) * (a[1] - c[1]); if (Math.abs(d) < 1e-12) return false; const l1 = ((b[1] - c[1]) * (p[0] - c[0]) + (c[0] - b[0]) * (p[1] - c[1])) / d, l2 = ((c[1] - a[1]) * (p[0] - c[0]) + (a[0] - c[0]) * (p[1] - c[1])) / d; return l1 >= -1e-9 && l2 >= -1e-9 && 1 - l1 - l2 >= -1e-9; };The points at working-plane height where the daylight factor is worked out.
spacing = max(0.25, ceil(√(area ÷ max) × 20) ÷ 20) m; points at lo + s ÷ 2, then every s
Unless a spacing is given, it is rounded up to a multiple of 0.05 m; it then grows by 0.05 m until at most max points fit. Grid coordinates are built by repeated addition of the spacing. Floor edges are matched on coordinates rounded to 1 mm; points closer to an edge than the border, less 1e-9 m, are dropped.
export function roomGrid(room, prods, { spacing = null, height = 0.85, border = 0.5, max = 400 } = {}) {
const sp = prods.get(room.id);
if (!sp) return { points: [], spacing: 0, floor: [], edges: [] };
const { tris } = floorOf(sp.tr);
if (!tris.length) return { points: [], spacing: 0, floor: [], edges: [] };
// the floor's outline: edges used by one floor triangle only
const key = (p) => Math.round(p[0] * 1000) + ',' + Math.round(p[1] * 1000);
const count = new Map();
for (const t of tris) for (let k = 0; k < 3; k++) { const a = t[k], b = t[(k + 1) % 3]; const e = [key(a), key(b)].sort().join('|'); count.set(e, (count.get(e) || { n: 0, a, b })); count.get(e).n++; }
const edges = [...count.values()].filter((e) => e.n === 1).map((e) => [e.a, e.b]);
let lo = [Infinity, Infinity], hi = [-Infinity, -Infinity], area = 0;
for (const t of tris) { for (let k = 0; k < 3; k++) { lo = [Math.min(lo[0], t[k][0]), Math.min(lo[1], t[k][1])]; hi = [Math.max(hi[0], t[k][0]), Math.max(hi[1], t[k][1])]; } area += Math.abs((t[1][0] - t[0][0]) * (t[2][1] - t[0][1]) - (t[2][0] - t[0][0]) * (t[1][1] - t[0][1])) / 2; }
const place = (s, b) => {
const pts = [];
for (let x = lo[0] + s / 2; x < hi[0]; x += s) for (let y = lo[1] + s / 2; y < hi[1]; y += s) {
const p = [x, y], t = tris.find((tt) => inTri(p, tt));
if (!t) continue;
if (b > 0 && edges.some(([a, c]) => segDist(p, a, c) < b - 1e-9)) continue;
pts.push([x, y, t[3] + height]);
}
return pts;
};
let s = spacing && spacing > 0 ? spacing : Math.max(0.25, Math.ceil(Math.sqrt(area / max) * 20) / 20);
let b = border, pts = place(s, b);
if (!pts.length && b > 0.1) { b = 0.1; pts = place(s, b); }
while (pts.length > max) { s = Math.round((s + 0.05) * 100) / 100; pts = place(s, b); }
return { points: pts, spacing: s, border: b, floor: tris.map((t) => [t[0], t[1], t[2]]), edges };
}The rectangle that stands for a window's opening.
k = √(area ÷ (width × height)); half sizes = k × width ÷ 2 and k × height ÷ 2
The window's span and height are scaled by one 64-bit factor so the rectangle's area equals the window's stated area.
/** The window's opening as a rectangle on its outer face: centre, across (s), up, half sizes (keeping the window's area). */
function aperture(w) {
const wid = w.span[1] - w.span[0], hgt = w.hz;
const k = wid > 0 && hgt > 0 && w.area > 0 ? Math.sqrt(w.area / (wid * hgt)) : 1;
const mid = (w.span[0] + w.span[1]) / 2;
return { c: [w.origin[0] + w.s[0] * mid, w.origin[1] + w.s[1] * mid, w.origin[2]], hw: (wid * k) / 2, hh: (hgt * k) / 2 };
}The light a point gets through one window straight from the sky and from obstructions, in percent.
SC = 100 ÷ E_h × Σ over patches of L(θ) sin θ × dA cos ω ÷ r²
A midpoint-rule sum over N × N patches of the opening, N from 8 to 48 by the opening's size over its distance, in 64-bit. The altitude is asin(d_z ÷ r) with r the hypot of the same vector; it is not clamped. Patches hidden in the window's profile count at 0.1 of the sky's value.
export function pointComponents(P, w, prof, { sky = 'cie', obstruction = 0.1, patches = 0 } = {}) {
const ap = aperture(w), n = w.n, s = w.s;
// only the part of the opening above the point's own horizon lights the working plane from the sky: clip it first, so no
// patch straddles the horizon
const v0 = Math.max(-ap.hh, P[2] - ap.c[2]), v1 = ap.hh;
if (v1 <= v0) return { sc: 0, erc: 0 };
const dist = Math.hypot(ap.c[0] - P[0], ap.c[1] - P[1], ap.c[2] - P[2]);
const N = patches || Math.max(8, Math.min(48, Math.ceil((24 * 2 * Math.max(ap.hw, ap.hh)) / Math.max(dist, 1e-3)))); // within 0.1 % of a converged sum (tests)
const du = (2 * ap.hw) / N, dv = (v1 - v0) / N, dA = du * dv;
let sc = 0, erc = 0;
for (let i = 0; i < N; i++) for (let j = 0; j < N; j++) {
const u = -ap.hw + (i + 0.5) * du, v = v0 + (j + 0.5) * dv;
const q = [ap.c[0] + s[0] * u, ap.c[1] + s[1] * u, ap.c[2] + v];
const d = [q[0] - P[0], q[1] - P[1], q[2] - P[2]], r = Math.hypot(d[0], d[1], d[2]);
const cosW = (d[0] * n[0] + d[1] * n[1]) / r;
if (cosW <= 0 || d[2] <= 0) continue; // behind the glass, or looking down at the ground (the inside reflected part has that)
const alt = Math.asin(d[2] / r), dw = (dA * cosW) / (r * r);
const lum = sky === 'uniform' ? 1 : L(alt);
const e = lum * Math.sin(alt) * dw;
const phi = Math.atan2(d[0] * s[0] + d[1] * s[1], d[0] * n[0] + d[1] * n[1]) * R2D;
if (prof && blockedAt(prof, phi, alt * R2D)) erc += obstruction * e; else sc += e;
}
const eh = sky === 'uniform' ? Math.PI : EH;
return { sc: (sc / eh) * 100, erc: (erc / eh) * 100 };
}The coefficient C of the internally reflected component.
straight lines between table rows at 0°, 10° … 80°
Linear interpolation in a 9-row table of whole numbers, with the angle clamped to 0…80°.
export function cOf(alpha) {
const a = Math.max(0, Math.min(80, alpha));
for (let i = 1; i < C_TABLE.length; i++) if (a <= C_TABLE[i][0]) { const [x0, y0] = C_TABLE[i - 1], [x1, y1] = C_TABLE[i]; return y0 + ((a - x0) / (x1 - x0)) * (y1 - y0); }
return 5;
}The light a room gets by reflection from its own surfaces, in percent.
IRC = T W ÷ (A (1 − R)) × (C R_fw + 5 R_cw) × M
One 64-bit expression. It returns 0 when there is no glazing, no room area or a mean reflectance of 1 or more, which guards the division by 1 − R. C is the mean of the table value weighted by window area.
export function irc(room, { T = DEFAULTS.T, glazedShare = DEFAULTS.glazedShare, M = DEFAULTS.M, refl = DEFAULTS.refl } = {}) {
const A = room.areas, R = meanReflectance(A, refl);
const Wg = room.windows.reduce((t, w) => t + w.area * glazedShare, 0);
if (!(Wg > 0) || !(room.A > 0) || R >= 1) return 0;
const C = room.windows.reduce((t, w) => t + w.area * cOf(w.alpha ?? 0), 0) / room.windows.reduce((t, w) => t + w.area, 0);
const rfw = (A.floor * refl.floor + 0.5 * A.wall * refl.wall) / (A.floor + 0.5 * A.wall || 1), rcw = (A.ceiling * refl.ceiling + 0.5 * A.wall * refl.wall) / (A.ceiling + 0.5 * A.wall || 1);
return ((T * Wg) / (room.A * (1 - R))) * (C * rfw + 5 * rcw) * M;
}The daylight factor at every grid point of a room and its mean, median, minimum, maximum and uniformity.
DF = (SC + ERC) × T × glazed share × M + IRC; uniformity = min ÷ mean
Each point's components are summed over the room's windows in 64-bit. The mean, median, minimum, maximum and the share at or above 2 % are taken from the sorted values without rounding.
export function gridDaylight(g, rooms, { T = DEFAULTS.T, glazedShare = DEFAULTS.glazedShare, M = DEFAULTS.M, refl = DEFAULTS.refl, height = 0.85, border = 0.5, spacing = null, max = 400, obstruction = 0.1 } = {}) {
return rooms.map((r) => {
const grid = roomGrid(r, g.prods, { spacing, height, border, max });
const ir = irc(r, { T, glazedShare, M, refl });
const pts = grid.points.map(([x, y, z]) => {
let sc = 0, erc = 0;
for (const w of r.windows) { const c = pointComponents([x, y, z], w, w.prof, { obstruction }); sc += c.sc; erc += c.erc; }
const f = T * glazedShare * M;
return { x, y, z, sc: sc * f, erc: erc * f, df: (sc + erc) * f + ir };
});
const dfs = pts.map((p) => p.df).sort((a, b) => a - b);
const mean = dfs.length ? dfs.reduce((t, x) => t + x, 0) / dfs.length : 0;
const median = dfs.length ? (dfs.length % 2 ? dfs[(dfs.length - 1) / 2] : (dfs[dfs.length / 2 - 1] + dfs[dfs.length / 2]) / 2) : 0;
return { ...r, irc: ir, grid: { points: pts, spacing: grid.spacing, border: grid.border, floor: grid.floor, edges: grid.edges, mean, median, min: dfs[0] ?? 0, max: dfs[dfs.length - 1] ?? 0,
uniformity: mean > 0 ? (dfs[0] ?? 0) / mean : 0, over2: dfs.length ? dfs.filter((x) => x >= 2).length / dfs.length : 0 } };
});
}The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.