Stock simulation and machining time

constructor(model, { cell = 0, tolV = 0.02, tolLat = null } = {}) {
  this.model = model;
  const s = model.stock;
  const w = s.max[0] - s.min[0], h = s.max[1] - s.min[1];
  this.cell = cell > 0 ? cell : Math.min(2, Math.max(0.1, Math.max(w, h) / 500));
  this.nx = Math.max(1, Math.ceil(w / this.cell));
  this.ny = Math.max(1, Math.ceil(h / this.cell));
  if (this.nx * this.ny > 4e6) throw new Error(t('cam.sim.grid_too_fine'));
  this.x0 = s.min[0]; this.y0 = s.min[1];
  this.top = s.max[2]; this.bottom = s.min[2];
  const n = this.nx * this.ny;
  this.z = new Float64Array(n);
  this.inStock = new Uint8Array(n);
  for (let j = 0; j < this.ny; j++) for (let i = 0; i < this.nx; i++) {
    const [x, y] = this.centre(i, j);
    const ins = s.kind === 'cylinder' ? Math.hypot(x - s.c[0], y - s.c[1]) <= s.r : x <= s.max[0] && y <= s.max[1];
    this.inStock[j * this.nx + i] = ins ? 1 : 0;
    this.z[j * this.nx + i] = ins ? this.top : -Infinity;
  }
  const diag = Math.hypot(...model.box.max.map((v, k) => v - model.box.min[k]));
  this.tolV = tolV;
  this.tolLat = tolLat ?? 2e-4 * diag + 0.01;   // the exact kernel's display mesh: linear deflection 2e-4 × diagonal
  this.hp = null;
  this.gouge = new Uint8Array(n);
  this.log = [];
}

/** Lower the cells under one flat XY track (a prim) at height z. → cells lowered by more than tolV */
stamp(prim, z, prof) {
  const R = prof.R, b = O.boxOf(prim);
  const i0 = Math.max(0, Math.floor((b[0] - R - this.x0) / this.cell)), i1 = Math.min(this.nx - 1, Math.floor((b[2] + R - this.x0) / this.cell));
  const j0 = Math.max(0, Math.floor((b[1] - R - this.y0) / this.cell)), j1 = Math.min(this.ny - 1, Math.floor((b[3] + R - this.y0) / this.cell));
  let hit = 0;
  for (let j = j0; j <= j1; j++) for (let i = i0; i <= i1; i++) {
    const k = j * this.nx + i;
    if (!this.inStock[k]) continue;
    const zc = z + prof.h(O.distToPrim(this.centre(i, j), prim));
    if (zc < this.z[k]) { if (this.z[k] - zc > this.tolV) hit++; this.z[k] = zc; }
  }
  return hit;
}

/** Sweep one segment ({k, a, b, c?, sw?}). */
sweep(s, prof) {
  const dz = s.b[2] - s.a[2];
  const flat = Math.abs(dz) < 1e-9;
  const xyLen = s.k === 'arc' ? Math.abs(s.sw) * Math.hypot(s.a[0] - s.c[0], s.a[1] - s.c[1]) : Math.hypot(s.b[0] - s.a[0], s.b[1] - s.a[1]);
  const primOf = (t0, t1) => {
    if (s.k === 'arc') {
      const r = Math.hypot(s.a[0] - s.c[0], s.a[1] - s.c[1]), a0 = Math.atan2(s.a[1] - s.c[1], s.a[0] - s.c[0]);
      return O.arc(s.c, r, a0 + s.sw * t0, s.sw * (t1 - t0));
    }
    const p = (t) => [s.a[0] + (s.b[0] - s.a[0]) * t, s.a[1] + (s.b[1] - s.a[1]) * t];
    return O.line(p(t0), p(t1));
  };
  if (flat) return this.stamp(primOf(0, 1), s.a[2], prof);
  if (xyLen < 1e-9) return this.stamp(O.line([s.a[0], s.a[1]], [s.a[0], s.a[1]]), Math.min(s.a[2], s.b[2]), prof);
  const n = Math.max(1, Math.ceil(Math.abs(dz) / DZ));
  let hit = 0;
  for (let i = 0; i < n; i++) hit += this.stamp(primOf(i / n, (i + 1) / n), Math.min(s.a[2] + (dz * i) / n, s.a[2] + (dz * (i + 1)) / n), prof);
  return hit;
}

/** Stock volume (above the stock bottom). */
volume() {
  let v = 0;
  for (let k = 0; k < this.z.length; k++) if (this.inStock[k] && this.z[k] > this.bottom) v += this.z[k] - this.bottom;
  return v * this.cell * this.cell;
}

/** The part's height at a cell with the lateral tolerance (min over the centre and 8 points tolLat away). */
erodedHeight(i, j) {
  const [x, y] = this.centre(i, j);
  let m = this.partHeights()[j * this.nx + i];
  for (let k = 0; k < 8 && m > -Infinity; k++) { const a = (k * Math.PI) / 4; m = Math.min(m, this.model.heightAt(x + this.tolLat * Math.cos(a), y + this.tolLat * Math.sin(a))); }
  return m;
}

/** Mark new gouge cells → {count (new), at: [[x, y, depth]] up to 50 spots} */
findGouges() {
  const hp = this.partHeights();
  let count = 0;
  const at = [];
  for (let j = 0; j < this.ny; j++) for (let i = 0; i < this.nx; i++) {
    const k = j * this.nx + i;
    if (this.gouge[k] || !(hp[k] > -Infinity) || this.z[k] >= hp[k] - this.tolV) continue;
    const he = this.erodedHeight(i, j);
    if (!(this.z[k] < he - this.tolV)) continue;
    this.gouge[k] = 1; count++;
    if (at.length < 50) { const c = this.centre(i, j); at.push([c[0], c[1], he - this.z[k]]); }
  }
  return { count, at };
}

/** Summary: remaining stock above the part, deepest gouge, cells, method. */
report() {
  const hp = this.partHeights();
  let remaining = 0, maxLeft = 0, gouges = 0, deepest = 0;
  for (let j = 0; j < this.ny; j++) for (let i = 0; i < this.nx; i++) {
    const k = j * this.nx + i;
    if (!this.inStock[k]) continue;
    const floor = hp[k] > -Infinity ? Math.max(hp[k], this.bottom) : this.bottom;
    const left = this.z[k] - floor;
    if (left > this.tolV) { remaining += left; maxLeft = Math.max(maxLeft, left); }
    if (this.gouge[k]) { gouges++; deepest = Math.max(deepest, this.erodedHeight(i, j) - this.z[k]); }
  }
  const a = this.cell * this.cell;
  return { remaining: remaining * a, maxLeft, gouges, deepest, cell: this.cell, nx: this.nx, ny: this.ny, tolV: this.tolV, tolLat: this.tolLat, volume: this.volume() };
}

/** Machining time of a program: feed moves at their feed, rapids at the setup's rapid rate, tool changes. */
export function cycleTime(program, setup) {
  const rapid = setup.machine.rapid, tc = setup.machine.toolChange;
  let feed = 0, rap = 0, changes = 0, prevTool = null;
  const perOp = {};
  for (const ph of program.paths) {
    if (!ph.moves || !ph.moves.length) continue;
    if (ph.tool && ph.tool.id !== prevTool) { if (prevTool !== null) changes++; prevTool = ph.tool.id; }
    let f = 0, r = 0;
    for (const s of segments([ph])) {
      const L = segLength(s);
      if (s.k === 'rapid') r += L / rapid; else f += L / s.f;
    }
    perOp[ph.op.id] = (f + r) * 60;
    feed += f; rap += r;
  }
  return { seconds: (feed + rap) * 60 + changes * tc, feed: feed * 60, rapid: rap * 60, toolChanges: changes, perOp, rapidRate: rapid, toolChangeTime: tc };
}