// =====================================================================
// Phong tessellation. Boubekeur & Alexa, "Phong Tessellation" (2008).
//
// Cheaper than PN triangles -- needs only the 3 corner positions and
// 3 corner normals (no Bezier control point precompute), and the
// per-vertex evaluation is one linear interpolate plus three plane
// projections. 18 floats per source tri vs PN's 48; same crack-free
// C0 surface across edges (since the projections agree along the edge
// where two adjacent corners' tangent planes meet).
//
// Bary convention here matches CalcPlainBaryPtsForTess: u = 1-x-y for
// corner P0 = (0,0), v = y for P1 = (0,1), w = x for P2 = (1,0).
//
// At barycentric (u, v, w) (u+v+w=1):
//   P_lin   = u*P0 + v*P1 + w*P2
//   For each corner i: pi_i = P_lin - dot(P_lin - Pi, Ni) * Ni
//                              (= projection of P_lin onto the tangent
//                               plane through Pi)
//   P_phong = u*pi_0 + v*pi_1 + w*pi_2
//   N_phong = normalize(u*N0 + v*N1 + w*N2)
//
// Setup is just the six vec3; no precomputed control points.
// =====================================================================
function phongSetup(p0, p1, p2, n0, n1, n2) {
  return { p0, p1, p2, n0, n1, n2 };
}

function phongEvalPos(s, u, v, w) {
  // Linear interpolate.
  const lx = u*s.p0[0] + v*s.p1[0] + w*s.p2[0];
  const ly = u*s.p0[1] + v*s.p1[1] + w*s.p2[1];
  const lz = u*s.p0[2] + v*s.p1[2] + w*s.p2[2];
  // Project onto tangent plane at P0.
  const d0 = (lx - s.p0[0])*s.n0[0] + (ly - s.p0[1])*s.n0[1] + (lz - s.p0[2])*s.n0[2];
  const x0 = lx - d0*s.n0[0]; const y0 = ly - d0*s.n0[1]; const z0 = lz - d0*s.n0[2];
  // Project onto tangent plane at P1.
  const d1 = (lx - s.p1[0])*s.n1[0] + (ly - s.p1[1])*s.n1[1] + (lz - s.p1[2])*s.n1[2];
  const x1 = lx - d1*s.n1[0]; const y1 = ly - d1*s.n1[1]; const z1 = lz - d1*s.n1[2];
  // Project onto tangent plane at P2.
  const d2 = (lx - s.p2[0])*s.n2[0] + (ly - s.p2[1])*s.n2[1] + (lz - s.p2[2])*s.n2[2];
  const x2 = lx - d2*s.n2[0]; const y2 = ly - d2*s.n2[1]; const z2 = lz - d2*s.n2[2];
  // Barycentric blend of the three projections.
  return [
    u*x0 + v*x1 + w*x2,
    u*y0 + v*y1 + w*y2,
    u*z0 + v*z1 + w*z2,
  ];
}

function phongEvalNrm(s, u, v, w) {
  const x = u*s.n0[0] + v*s.n1[0] + w*s.n2[0];
  const y = u*s.n0[1] + v*s.n1[1] + w*s.n2[1];
  const z = u*s.n0[2] + v*s.n1[2] + w*s.n2[2];
  const l = Math.hypot(x, y, z) || 1;
  return [x/l, y/l, z/l];
}
