{"spec_id":"stereonet-equal-area","library":"muix","language":"javascript","code":"// anyplot.ai\n// stereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)\n// Library: muix 7.29.1 | JavaScript 22.22.3\n// Quality: 94/100 | Created: 2026-06-16\n//# anyplot-orientation: square\n// anyplot.ai\n// stereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)\n// Library: MUI X Charts | React | Node 22\n// License: @mui/x-charts — MIT (community). Pro/Premium are out of scope.\n// Quality: pending | Created: 2026-06-16\nimport { ChartContainer } from \"@mui/x-charts/ChartContainer\";\nimport { ScatterPlot } from \"@mui/x-charts/ScatterChart\";\nimport { ChartsTooltip } from \"@mui/x-charts/ChartsTooltip\";\nimport { useXScale, useYScale } from \"@mui/x-charts/hooks\";\n\nconst t = window.ANYPLOT_TOKENS;\nconst SIZE = window.ANYPLOT_SIZE;\nconst TITLE = \"stereonet-equal-area · javascript · muix · anyplot.ai\";\n\n// Axis half-extent: the unit primitive circle (radius 1) sits centred, leaving\n// top/bottom bands for the title and legend. Square drawing area + equal domains\n// keep the projection perfectly circular.\nconst R = 1.45;\n\n// --- Stereonet projection helpers ------------------------------------------\n// Lower-hemisphere Lambert equal-area (Schmidt) net. A line of plunge p (from\n// horizontal) and trend a (clockwise from North) maps to polar radius\n// r = √2·sin((90°−p)/2), normalised so the horizon (p=0) lands on the circle.\nfunction projectVec(n, e, d) {\n  // Force the direction into the lower hemisphere (poles/lines are undirected).\n  if (d < 0) {\n    n = -n;\n    e = -e;\n    d = -d;\n  }\n  const plunge = Math.asin(Math.max(-1, Math.min(1, d)));\n  const trend = Math.atan2(e, n);\n  const theta = Math.PI / 2 - plunge; // colatitude from the downward vertical\n  const r = Math.SQRT2 * Math.sin(theta / 2);\n  return { x: r * Math.sin(trend), y: r * Math.cos(trend) };\n}\n\n// Pole (normal) to a plane given its dip and dip direction.\nfunction poleVec(dipDeg, dipDirDeg) {\n  const pl = ((90 - dipDeg) * Math.PI) / 180;\n  const tr = ((dipDirDeg + 180) * Math.PI) / 180;\n  return {\n    n: Math.cos(pl) * Math.cos(tr),\n    e: Math.cos(pl) * Math.sin(tr),\n    d: Math.sin(pl),\n  };\n}\n\nfunction norm3(v) {\n  const m = Math.hypot(v[0], v[1], v[2]) || 1;\n  return [v[0] / m, v[1] / m, v[2] / m];\n}\nfunction cross3(a, b) {\n  return [\n    a[1] * b[2] - a[2] * b[1],\n    a[2] * b[0] - a[0] * b[2],\n    a[0] * b[1] - a[1] * b[0],\n  ];\n}\n\n// Great circle (the planar feature itself) as polyline segments, split where the\n// arc wraps across the primitive circle so the path never draws a stray chord.\nfunction greatCirclePath(pole) {\n  const P = norm3([pole.n, pole.e, pole.d]);\n  const seed = Math.abs(P[2]) < 0.9 ? [0, 0, 1] : [0, 1, 0];\n  const u = norm3(cross3(P, seed));\n  const v = cross3(P, u);\n  const segs = [];\n  let cur = [];\n  let prev = null;\n  const STEPS = 200;\n  for (let i = 0; i <= STEPS; i++) {\n    const phi = (i / STEPS) * 2 * Math.PI;\n    const cphi = Math.cos(phi);\n    const sphi = Math.sin(phi);\n    const c = [\n      cphi * u[0] + sphi * v[0],\n      cphi * u[1] + sphi * v[1],\n      cphi * u[2] + sphi * v[2],\n    ];\n    const pt = projectVec(c[0], c[1], c[2]);\n    if (prev && Math.hypot(pt.x - prev.x, pt.y - prev.y) > 0.25) {\n      if (cur.length > 1) segs.push(cur);\n      cur = [];\n    }\n    cur.push(pt);\n    prev = pt;\n  }\n  if (cur.length > 1) segs.push(cur);\n  return segs;\n}\n\n// Small circle (a cone of fixed half-angle about an axis) — used for the net's\n// reference graticule. Only the lower-hemisphere portion is kept.\nfunction smallCirclePath(axis, beta) {\n  const A = norm3(axis);\n  const seed = Math.abs(A[2]) < 0.9 ? [0, 0, 1] : [0, 1, 0];\n  const e1 = norm3(cross3(A, seed));\n  const e2 = cross3(A, e1);\n  const cb = Math.cos(beta);\n  const sb = Math.sin(beta);\n  const segs = [];\n  let cur = [];\n  let prev = null;\n  const STEPS = 200;\n  for (let i = 0; i <= STEPS; i++) {\n    const phi = (i / STEPS) * 2 * Math.PI;\n    const cphi = Math.cos(phi);\n    const sphi = Math.sin(phi);\n    const c = [\n      cb * A[0] + sb * (cphi * e1[0] + sphi * e2[0]),\n      cb * A[1] + sb * (cphi * e1[1] + sphi * e2[1]),\n      cb * A[2] + sb * (cphi * e1[2] + sphi * e2[2]),\n    ];\n    if (c[2] < 0) {\n      if (cur.length > 1) segs.push(cur);\n      cur = [];\n      prev = null;\n      continue;\n    }\n    const pt = projectVec(c[0], c[1], c[2]);\n    if (prev && Math.hypot(pt.x - prev.x, pt.y - prev.y) > 0.25) {\n      if (cur.length > 1) segs.push(cur);\n      cur = [];\n    }\n    cur.push(pt);\n    prev = pt;\n  }\n  if (cur.length > 1) segs.push(cur);\n  return segs;\n}\n\n// --- Data: a geological mapping campaign (deterministic, in-memory) ---------\n// Four fabric elements, each a clustered orientation population. Colours follow\n// the Imprint categorical order (Bedding is brand green, position 1).\nlet seed = 20260616 >>> 0;\nfunction rand() {\n  seed = (1664525 * seed + 1013904223) >>> 0;\n  return seed / 4294967296;\n}\nfunction randn() {\n  let u = 0;\n  let v = 0;\n  while (u === 0) u = rand();\n  while (v === 0) v = rand();\n  return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);\n}\n\nconst FEATURES = [\n  { name: \"Bedding\", color: t.palette[0], n: 46, dip: 24, dipDir: 118, sDip: 6, sDir: 14 },\n  { name: \"Joint set\", color: t.palette[1], n: 38, dip: 78, dipDir: 48, sDip: 5, sDir: 10 },\n  { name: \"Foliation\", color: t.palette[2], n: 34, dip: 62, dipDir: 305, sDip: 7, sDir: 12 },\n  { name: \"Fault\", color: t.palette[3], n: 26, dip: 54, dipDir: 212, sDip: 6, sDir: 11 },\n];\n\nconst SERIES = [];\nconst allPoles = [];\nfor (const f of FEATURES) {\n  const data = [];\n  for (let i = 0; i < f.n; i++) {\n    const dip = Math.max(2, Math.min(89, f.dip + randn() * f.sDip));\n    const dipDir = (((f.dipDir + randn() * f.sDir) % 360) + 360) % 360;\n    allPoles.push({ dip, dipDir });\n    const v = poleVec(dip, dipDir);\n    const pt = projectVec(v.n, v.e, v.d);\n    data.push({ x: pt.x, y: pt.y, id: `${f.name}-${i}` });\n  }\n  SERIES.push({ type: \"scatter\", label: f.name, color: f.color, markerSize: 5, data });\n}\n\n// Representative mean plane per feature, drawn as a great circle.\nconst MEAN_PLANES = FEATURES.map((f) => ({\n  color: f.color,\n  segs: greatCirclePath(poleVec(f.dip, f.dipDir)),\n}));\n\n// Equatorial Schmidt net graticule (meridians + small circles), kept subtle.\nconst NET = [];\nfor (let dip = 10; dip <= 80; dip += 10) {\n  NET.push(greatCirclePath(poleVec(dip, 90)));\n  NET.push(greatCirclePath(poleVec(dip, 270)));\n}\nfor (let beta = 10; beta <= 80; beta += 10) {\n  NET.push(smallCirclePath([1, 0, 0], (beta * Math.PI) / 180));\n  NET.push(smallCirclePath([1, 0, 0], ((180 - beta) * Math.PI) / 180));\n}\n\n// --- Kamb-style density field + contour lines -------------------------------\n// Exponential smoothing of pole axes over an equal-area grid, then marching\n// squares for iso-density lines. Coloured with the Imprint sequential ramp.\nconst NX = 80;\nconst NY = 80;\nconst gx = (i) => -1.0 + (i * 2) / (NX - 1);\nconst gy = (j) => -1.0 + (j * 2) / (NY - 1);\nconst poleUnits = allPoles.map((p) => {\n  const v = poleVec(p.dip, p.dipDir);\n  return norm3([v.n, v.e, v.d]);\n});\nconst K = 28;\nconst density = new Float64Array(NX * NY);\nfor (let j = 0; j < NY; j++) {\n  for (let i = 0; i < NX; i++) {\n    const X = gx(i);\n    const Y = gy(j);\n    const r = Math.hypot(X, Y);\n    if (r > 1.0) continue;\n    const theta = 2 * Math.asin(Math.min(1, r / Math.SQRT2));\n    const az = Math.atan2(X, Y);\n    const V = [\n      Math.sin(theta) * Math.cos(az),\n      Math.sin(theta) * Math.sin(az),\n      Math.cos(theta),\n    ];\n    let s = 0;\n    for (const p of poleUnits) {\n      const dot = V[0] * p[0] + V[1] * p[1] + V[2] * p[2];\n      s += Math.exp(K * (dot * dot - 1));\n    }\n    density[j * NX + i] = s;\n  }\n}\n\nfunction marchingSquares(values, nx, ny, level) {\n  const segs = [];\n  const at = (i, j) => values[j * nx + i];\n  const lerp = (x1, y1, v1, x2, y2, v2) => {\n    const tt = (level - v1) / (v2 - v1);\n    return [x1 + tt * (x2 - x1), y1 + tt * (y2 - y1)];\n  };\n  for (let j = 0; j < ny - 1; j++) {\n    for (let i = 0; i < nx - 1; i++) {\n      const v0 = at(i, j);\n      const v1 = at(i + 1, j);\n      const v2 = at(i + 1, j + 1);\n      const v3 = at(i, j + 1);\n      const x0 = gx(i);\n      const x1 = gx(i + 1);\n      const y0 = gy(j);\n      const y1 = gy(j + 1);\n      let code = 0;\n      if (v0 >= level) code |= 1;\n      if (v1 >= level) code |= 2;\n      if (v2 >= level) code |= 4;\n      if (v3 >= level) code |= 8;\n      if (code === 0 || code === 15) continue;\n      const eB = lerp(x0, y0, v0, x1, y0, v1);\n      const eR = lerp(x1, y0, v1, x1, y1, v2);\n      const eT = lerp(x1, y1, v2, x0, y1, v3);\n      const eL = lerp(x0, y1, v3, x0, y0, v0);\n      const push = (a, b) => segs.push([a[0], a[1], b[0], b[1]]);\n      switch (code) {\n        case 1:\n        case 14:\n          push(eL, eB);\n          break;\n        case 2:\n        case 13:\n          push(eB, eR);\n          break;\n        case 3:\n        case 12:\n          push(eL, eR);\n          break;\n        case 4:\n        case 11:\n          push(eR, eT);\n          break;\n        case 5:\n          push(eL, eT);\n          push(eB, eR);\n          break;\n        case 6:\n        case 9:\n          push(eB, eT);\n          break;\n        case 7:\n        case 8:\n          push(eL, eT);\n          break;\n        case 10:\n          push(eL, eB);\n          push(eR, eT);\n          break;\n        default:\n          break;\n      }\n    }\n  }\n  return segs;\n}\n\nfunction hexLerp(a, b, tt) {\n  const pa = [parseInt(a.slice(1, 3), 16), parseInt(a.slice(3, 5), 16), parseInt(a.slice(5, 7), 16)];\n  const pb = [parseInt(b.slice(1, 3), 16), parseInt(b.slice(3, 5), 16), parseInt(b.slice(5, 7), 16)];\n  const mix = pa.map((v, i) => Math.round(v + (pb[i] - v) * tt));\n  return \"#\" + mix.map((v) => v.toString(16).padStart(2, \"0\")).join(\"\");\n}\n\nconst maxD = density.reduce((m, v) => (v > m ? v : m), 0);\nconst LEVEL_FRACS = [0.16, 0.3, 0.48, 0.7];\nconst CONTOURS = LEVEL_FRACS.map((lf, idx) => ({\n  segs: marchingSquares(density, NX, NY, lf * maxD),\n  color: hexLerp(t.seq[0], t.seq[1], idx / (LEVEL_FRACS.length - 1)),\n  w: 1.4 + idx * 0.6,\n  op: 0.45 + idx * 0.16,\n}));\n\n// --- SVG overlays (data → pixels via the chart's scales) --------------------\nconst FONT = '-apple-system, BlinkMacSystemFont, \"Segoe UI\", Roboto, Helvetica, Arial, sans-serif';\n\nfunction polyToPath(xs, ys, segs) {\n  return segs\n    .map(\n      (poly) =>\n        \"M\" + poly.map((pt) => `${xs(pt.x).toFixed(1)} ${ys(pt.y).toFixed(1)}`).join(\" L\"),\n    )\n    .join(\" \");\n}\n\nfunction NetGrid() {\n  const xs = useXScale();\n  const ys = useYScale();\n  const cx = xs(0);\n  const cy = ys(0);\n  const rPx = Math.abs(xs(1) - xs(0));\n  return (\n    <g>\n      {NET.map((segs, k) => (\n        <path key={k} d={polyToPath(xs, ys, segs)} fill=\"none\" stroke={t.grid} strokeWidth={1} />\n      ))}\n      <line x1={xs(-1)} y1={cy} x2={xs(1)} y2={cy} stroke={t.grid} strokeWidth={1} />\n      <line x1={cx} y1={ys(-1)} x2={cx} y2={ys(1)} stroke={t.grid} strokeWidth={1} />\n      <circle cx={cx} cy={cy} r={rPx} fill=\"none\" stroke={t.ink} strokeWidth={2.5} />\n    </g>\n  );\n}\n\nfunction DensityContours() {\n  const xs = useXScale();\n  const ys = useYScale();\n  return (\n    <g>\n      {CONTOURS.map((c, k) => (\n        <path\n          key={k}\n          d={c.segs\n            .map((s) => `M${xs(s[0]).toFixed(1)} ${ys(s[1]).toFixed(1)} L${xs(s[2]).toFixed(1)} ${ys(s[3]).toFixed(1)}`)\n            .join(\" \")}\n          fill=\"none\"\n          stroke={c.color}\n          strokeWidth={c.w}\n          strokeOpacity={c.op}\n          strokeLinecap=\"round\"\n        />\n      ))}\n    </g>\n  );\n}\n\nfunction GreatCircles() {\n  const xs = useXScale();\n  const ys = useYScale();\n  return (\n    <g>\n      {MEAN_PLANES.map((mp, k) => (\n        <path\n          key={k}\n          d={polyToPath(xs, ys, mp.segs)}\n          fill=\"none\"\n          stroke={mp.color}\n          strokeWidth={3.5}\n          strokeOpacity={0.95}\n          strokeLinejoin=\"round\"\n          strokeLinecap=\"round\"\n        />\n      ))}\n    </g>\n  );\n}\n\nconst CARDINALS = [\n  [\"N\", 0],\n  [\"E\", 90],\n  [\"S\", 180],\n  [\"W\", 270],\n];\nconst DEG_LABELS = [30, 60, 120, 150, 210, 240, 300, 330];\n\nfunction Chrome() {\n  const xs = useXScale();\n  const ys = useYScale();\n  const radial = (deg, rr) => {\n    const a = (deg * Math.PI) / 180;\n    return [xs(rr * Math.sin(a)), ys(rr * Math.cos(a))];\n  };\n  const ticks = [];\n  for (let deg = 0; deg < 360; deg += 10) {\n    const major = deg % 30 === 0;\n    const [x1, y1] = radial(deg, 1.0);\n    const [x2, y2] = radial(deg, major ? 1.04 : 1.022);\n    ticks.push(\n      <line key={`t${deg}`} x1={x1} y1={y1} x2={x2} y2={y2} stroke={t.inkSoft} strokeWidth={major ? 1.6 : 1} />,\n    );\n  }\n  return (\n    <g fontFamily={FONT}>\n      {ticks}\n      {DEG_LABELS.map((deg) => {\n        const [x, y] = radial(deg, 1.11);\n        return (\n          <text key={`d${deg}`} x={x} y={y} textAnchor=\"middle\" dominantBaseline=\"middle\" fontSize={13} fill={t.inkSoft}>\n            {deg}°\n          </text>\n        );\n      })}\n      {CARDINALS.map(([label, deg]) => {\n        const [x, y] = radial(deg, 1.085);\n        return (\n          <text key={label} x={x} y={y} textAnchor=\"middle\" dominantBaseline=\"middle\" fontSize={24} fontWeight={700} fill={t.ink}>\n            {label}\n          </text>\n        );\n      })}\n    </g>\n  );\n}\n\nfunction Annotations() {\n  const xs = useXScale();\n  const ys = useYScale();\n  const legendStart = -1.18;\n  const legendStep = 0.62;\n  return (\n    <g fontFamily={FONT}>\n      <text x={xs(0)} y={ys(1.34)} textAnchor=\"middle\" dominantBaseline=\"middle\" fontSize={25} fontWeight={600} fill={t.ink}>\n        {TITLE}\n      </text>\n      <text x={xs(0)} y={ys(1.2)} textAnchor=\"middle\" dominantBaseline=\"middle\" fontSize={15} fill={t.inkSoft}>\n        Lower-hemisphere equal-area (Schmidt) net · poles to planes with mean great circles\n      </text>\n      {FEATURES.map((f, i) => {\n        const x = legendStart + i * legendStep;\n        return (\n          <g key={f.name}>\n            <circle cx={xs(x)} cy={ys(-1.18)} r={7} fill={f.color} />\n            <text x={xs(x) + 14} y={ys(-1.18)} dominantBaseline=\"middle\" fontSize={15} fill={t.inkSoft}>\n              {f.name}\n            </text>\n          </g>\n        );\n      })}\n      <text x={xs(0)} y={ys(-1.34)} textAnchor=\"middle\" dominantBaseline=\"middle\" fontSize={13} fill={t.inkSoft}>\n        Density contours: Imprint sequential ramp (exponential Kamb smoothing of pole axes)\n      </text>\n    </g>\n  );\n}\n\n// --- Chart (default-exported component — the harness mounts it) -------------\nexport default function Chart() {\n  return (\n    <ChartContainer\n      width={SIZE.width}\n      height={SIZE.height}\n      margin={{ top: 12, right: 12, bottom: 12, left: 12 }}\n      series={SERIES}\n      xAxis={[{ scaleType: \"linear\", min: -R, max: R }]}\n      yAxis={[{ scaleType: \"linear\", min: -R, max: R }]}\n      skipAnimation\n    >\n      <NetGrid />\n      <DensityContours />\n      <GreatCircles />\n      <ScatterPlot />\n      <Chrome />\n      <Annotations />\n      <ChartsTooltip trigger=\"item\" />\n    </ChartContainer>\n  );\n}\n"}