{"spec_id":"stereonet-equal-area","library":"highcharts","language":"javascript","code":"// anyplot.ai\n// stereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)\n// Library: highcharts 12.6.0 | JavaScript 22.22.3\n// Quality: 88/100 | Created: 2026-06-16\n//# anyplot-orientation: square\n// anyplot.ai\n// stereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)\n// Library: Highcharts 12.6.0 | Node 22\n// License: Highcharts — commercial license, free for non-commercial use (highcharts.com/license)\n// Quality: pending | Created: 2026-06-16\n\nconst t = window.ANYPLOT_TOKENS;\nconst DEG = Math.PI / 180;\nconst R = 1.0; // primitive (perimeter) radius in axis units\n\n// --- Deterministic RNG (fixed-seed LCG + Box–Muller) -----------------------\nlet _seed = 20260616;\nfunction rand() {\n  _seed = (_seed * 1103515245 + 12345) & 0x7fffffff;\n  return _seed / 0x7fffffff;\n}\nfunction gauss(mean, sd) {\n  let u = 0,\n    v = 0;\n  while (u === 0) u = rand();\n  while (v === 0) v = rand();\n  return mean + sd * Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v);\n}\n\n// --- Equal-area (Schmidt, lower-hemisphere) projection ----------------------\n// A line of given trend/plunge maps to a planar point. North = +y (top),\n// East = +x (right); azimuth measured clockwise from North.\nfunction projLine(trendDeg, plungeDeg) {\n  const tr = trendDeg * DEG;\n  const rho = (90 - plungeDeg) * DEG; // angular distance from net centre\n  const r = R * Math.SQRT2 * Math.sin(rho / 2);\n  return [r * Math.sin(tr), r * Math.cos(tr)];\n}\n// Project a downward direction vector (E, N, U) with U <= 0.\nfunction projVec(e, n, u) {\n  const h = Math.hypot(e, n);\n  const rho = Math.acos(Math.max(-1, Math.min(1, -u))); // angle from nadir\n  const r = R * Math.SQRT2 * Math.sin(rho / 2);\n  if (h < 1e-9) return [0, 0];\n  return [(r * e) / h, (r * n) / h];\n}\n// Pole to a plane: normal plunges (90 - dip), trend = dipDir + 180.\nfunction poleXY(dipDir, dip) {\n  return projLine(dipDir + 180, 90 - dip);\n}\n// Great circle of a plane: sweep in-plane lines from one strike end to the other.\nfunction greatCircle(dipDir, dip, steps) {\n  const a = dipDir * DEG,\n    d = dip * DEG;\n  const sx = Math.sin(a - Math.PI / 2),\n    sy = Math.cos(a - Math.PI / 2); // strike line (horizontal)\n  const tx = Math.cos(d) * Math.sin(a),\n    ty = Math.cos(d) * Math.cos(a),\n    tz = -Math.sin(d); // dip vector\n  const pts = [];\n  for (let i = 0; i <= steps; i++) {\n    const b = (Math.PI * i) / steps;\n    const cb = Math.cos(b),\n      sb = Math.sin(b);\n    pts.push(projVec(cb * sx + sb * tx, cb * sy + sb * ty, sb * tz));\n  }\n  return pts;\n}\n\n// --- Field data: bedding, joint set, fault (deterministic) ------------------\nconst features = [\n  { name: \"Bedding\", color: t.palette[0], symbol: \"circle\", count: 46, ddMean: 118, ddSd: 11, dipMean: 24, dipSd: 6, nGreat: 7 },\n  { name: \"Joint set\", color: t.palette[1], symbol: \"triangle\", count: 34, ddMean: 212, ddSd: 14, dipMean: 74, dipSd: 8, nGreat: 8 },\n  { name: \"Fault\", color: t.palette[2], symbol: \"diamond\", count: 17, ddMean: 47, ddSd: 12, dipMean: 60, dipSd: 11, nGreat: 8 },\n];\n\nconst clamp = (x, lo, hi) => Math.max(lo, Math.min(hi, x));\n\n// Generate every measurement first; collect all poles for the density pass.\nconst allPoles = [];\nfeatures.forEach((f) => {\n  f.planes = [];\n  for (let i = 0; i < f.count; i++) {\n    const dipDir = (gauss(f.ddMean, f.ddSd) + 360) % 360;\n    const dip = clamp(gauss(f.dipMean, f.dipSd), 2, 88);\n    f.planes.push({ dipDir, dip });\n    allPoles.push(poleXY(dipDir, dip));\n  }\n});\n\n// --- Kamb density contours over the pole population -------------------------\n// Kamb (1959): count poles falling inside a counting circle whose area is sized\n// so the expected uniform count is 3σ. Contours are drawn at successive σ levels\n// above that uniform background via marching squares over an equal-area grid.\nfunction kambContours(poles) {\n  const n = poles.length;\n  const rc = Math.sqrt(9 / (n + 9)); // counting-circle radius (3σ), primitive R=1\n  const rc2 = rc * rc;\n  const E = n * rc2; // expected count under a uniform distribution\n  const sigma = Math.sqrt(E * (1 - rc2));\n  const NX = 120,\n    NY = 120,\n    span = 2.1,\n    x0 = -1.05,\n    y0 = -1.05;\n  const dx = span / (NX - 1),\n    dy = span / (NY - 1);\n  const grid = new Float64Array(NX * NY);\n  let maxC = 0;\n  for (let j = 0; j < NY; j++) {\n    const gy = y0 + j * dy;\n    for (let i = 0; i < NX; i++) {\n      const gx = x0 + i * dx;\n      if (gx * gx + gy * gy > 1.0) continue; // only count inside the primitive\n      let c = 0;\n      for (const p of poles) {\n        const ex = p[0] - gx,\n          ey = p[1] - gy;\n        if (ex * ex + ey * ey <= rc2) c++;\n      }\n      grid[j * NX + i] = c;\n      if (c > maxC) maxC = c;\n    }\n  }\n  // Contour levels in σ above the uniform background; keep those the data reaches.\n  let levels = [2, 4, 6, 8].map((k) => E + k * sigma).filter((l) => l < maxC);\n  if (levels.length === 0) levels = [0.5 * maxC];\n\n  const frac = (a, b, level) => (level - a) / (b - a);\n  function edgePt(e, v0, v1, v2, v3, level, xL, yB) {\n    if (e === 0) return [xL + dx * frac(v0, v1, level), yB]; // bottom: BL→BR\n    if (e === 1) return [xL + dx, yB + dy * frac(v1, v2, level)]; // right: BR→TR\n    if (e === 2) return [xL + dx * frac(v3, v2, level), yB + dy]; // top: TL→TR\n    return [xL, yB + dy * frac(v0, v3, level)]; // left: BL→TL\n  }\n  // Marching-squares segment table (edges 0=bottom 1=right 2=top 3=left).\n  const tbl = {\n    1: [[3, 0]], 2: [[0, 1]], 3: [[3, 1]], 4: [[1, 2]],\n    5: [[3, 0], [1, 2]], 6: [[0, 2]], 7: [[3, 2]], 8: [[2, 3]],\n    9: [[0, 2]], 10: [[0, 1], [2, 3]], 11: [[1, 2]], 12: [[3, 1]],\n    13: [[0, 1]], 14: [[3, 0]],\n  };\n  const data = [];\n  for (const level of levels) {\n    for (let j = 0; j < NY - 1; j++) {\n      for (let i = 0; i < NX - 1; i++) {\n        const v0 = grid[j * NX + i],\n          v1 = grid[j * NX + i + 1],\n          v2 = grid[(j + 1) * NX + i + 1],\n          v3 = grid[(j + 1) * NX + i];\n        let idx = 0;\n        if (v0 > level) idx |= 1;\n        if (v1 > level) idx |= 2;\n        if (v2 > level) idx |= 4;\n        if (v3 > level) idx |= 8;\n        const segs = tbl[idx];\n        if (!segs) continue;\n        const xL = x0 + i * dx,\n          yB = y0 + j * dy;\n        for (const [a, b] of segs) {\n          data.push(edgePt(a, v0, v1, v2, v3, level, xL, yB));\n          data.push(edgePt(b, v0, v1, v2, v3, level, xL, yB));\n          data.push([null, null]);\n        }\n      }\n    }\n  }\n  return data;\n}\nconst contourData = kambContours(allPoles);\n\n// --- Subtle equal-area net graticule (10° spacing about the N–S axis) -------\nconst netData = [];\nfunction pushPath(pts) {\n  for (const p of pts) netData.push(p);\n  netData.push([null, null]);\n}\n// Meridians: planes striking N–S, dipping E/W — they pass through N and S.\nfor (let dip = 10; dip <= 80; dip += 10) {\n  pushPath(greatCircle(90, dip, 90));\n  pushPath(greatCircle(270, dip, 90));\n}\npushPath(greatCircle(90, 90, 90)); // straight N–S diameter\npushPath(greatCircle(0, 90, 90)); // straight E–W diameter\n// Small circles: cones about the horizontal N–S axis (latitude lines).\nfor (let colat = 10; colat <= 170; colat += 10) {\n  if (colat === 90) continue; // coincides with the E–W diameter\n  const k = colat * DEG;\n  const arc = [];\n  for (let j = 0; j <= 90; j++) {\n    const phi = Math.PI + (Math.PI * j) / 90; // lower hemisphere (U <= 0)\n    arc.push(projVec(Math.sin(k) * Math.cos(phi), Math.cos(k), Math.sin(k) * Math.sin(phi)));\n  }\n  pushPath(arc);\n}\n\n// --- Primitive circle + 10° perimeter tick marks ---------------------------\nconst frameData = [];\nfor (let deg = 0; deg <= 360; deg += 2) {\n  frameData.push([R * Math.sin(deg * DEG), R * Math.cos(deg * DEG)]);\n}\nframeData.push([null, null]);\nfor (let deg = 0; deg < 360; deg += 10) {\n  const outer = deg % 90 === 0 ? 1.05 : 1.03;\n  frameData.push([R * Math.sin(deg * DEG), R * Math.cos(deg * DEG)]);\n  frameData.push([outer * Math.sin(deg * DEG), outer * Math.cos(deg * DEG)]);\n  frameData.push([null, null]);\n}\n\n// --- Build per-feature pole (scatter) + great-circle (line) series ----------\nconst series = [];\n\n// Net graticule and frame sit beneath the data.\nseries.push({\n  type: \"line\",\n  name: \"Net grid\",\n  data: netData,\n  color: t.grid,\n  lineWidth: 1.1,\n  enableMouseTracking: false,\n  showInLegend: false,\n  marker: { enabled: false },\n  zIndex: 0,\n});\nseries.push({\n  type: \"line\",\n  name: \"Primitive\",\n  data: frameData,\n  color: t.inkSoft,\n  lineWidth: 1.6,\n  enableMouseTracking: false,\n  showInLegend: false,\n  marker: { enabled: false },\n  zIndex: 1,\n});\n\nfeatures.forEach((f, idx) => {\n  // Poles (scatter) — every measurement; clustering reveals preferred orientation.\n  const poleId = \"poles-\" + idx;\n  series.push({\n    type: \"scatter\",\n    id: poleId,\n    name: f.name,\n    data: f.planes.map((p) => poleXY(p.dipDir, p.dip)),\n    color: f.color,\n    marker: {\n      symbol: f.symbol,\n      radius: 5,\n      lineWidth: 1,\n      lineColor: t.pageBg,\n    },\n    zIndex: 5,\n  });\n\n  // Great circles (line) — an evenly sampled subset keeps the net legible.\n  const gcData = [];\n  const stride = Math.max(1, Math.floor(f.count / f.nGreat));\n  for (let i = 0; i < f.count; i += stride) {\n    for (const pt of greatCircle(f.planes[i].dipDir, f.planes[i].dip, 80)) gcData.push(pt);\n    gcData.push([null, null]);\n  }\n  series.push({\n    type: \"line\",\n    linkedTo: poleId,\n    name: f.name + \" planes\",\n    data: gcData,\n    color: f.color,\n    lineWidth: 1.9,\n    opacity: 0.7,\n    enableMouseTracking: false,\n    marker: { enabled: false },\n    zIndex: 3,\n  });\n});\n\n// Kamb density contours — nested isolines over the pole clusters (zIndex below\n// the poles, above the graticule) directly mark preferred orientations.\nseries.push({\n  type: \"line\",\n  name: \"Kamb density\",\n  data: contourData,\n  color: t.inkSoft,\n  opacity: 0.6,\n  lineWidth: 1.3,\n  cropThreshold: Infinity, // unsorted marching-squares segments must not be cropped\n  enableMouseTracking: false,\n  marker: { enabled: false },\n  zIndex: 4, // above the graticule & great circles, beneath the pole markers\n});\n\n// --- Cardinal labels (N/E/S/W) just outside the primitive ------------------\nseries.push({\n  type: \"scatter\",\n  name: \"Cardinals\",\n  showInLegend: false,\n  enableMouseTracking: false,\n  color: t.ink,\n  marker: { enabled: false },\n  dataLabels: {\n    enabled: true,\n    allowOverlap: true,\n    style: { color: t.ink, fontSize: \"16px\", fontWeight: \"600\", textOutline: \"none\" },\n  },\n  data: [\n    { x: 0, y: 1.12, dataLabels: { format: \"N\" } },\n    { x: 1.13, y: 0, dataLabels: { format: \"E\" } },\n    { x: 0, y: -1.13, dataLabels: { format: \"S\" } },\n    { x: -1.13, y: 0, dataLabels: { format: \"W\" } },\n  ],\n  zIndex: 6,\n});\n\n// --- Chart ------------------------------------------------------------------\nHighcharts.chart(\"container\", {\n  chart: {\n    type: \"scatter\",\n    backgroundColor: \"transparent\",\n    animation: false,\n    style: { fontFamily: \"inherit\" },\n    // Symmetric margins make the plot box square (1000×1000) so the net is round.\n    marginTop: 70,\n    marginBottom: 130,\n    marginLeft: 100,\n    marginRight: 100,\n  },\n  credits: { enabled: false },\n  colors: t.palette,\n  title: {\n    text: \"stereonet-equal-area · javascript · highcharts · anyplot.ai\",\n    style: { color: t.ink, fontSize: \"22px\", fontWeight: \"600\" },\n  },\n  subtitle: {\n    text: \"Lower-hemisphere Schmidt net — poles, great circles & Kamb density contours\",\n    style: { color: t.inkSoft, fontSize: \"14px\" },\n  },\n  xAxis: {\n    min: -1.18,\n    max: 1.18,\n    visible: false,\n    startOnTick: false,\n    endOnTick: false,\n  },\n  yAxis: {\n    min: -1.18,\n    max: 1.18,\n    visible: false,\n    startOnTick: false,\n    endOnTick: false,\n    gridLineWidth: 0,\n  },\n  legend: {\n    align: \"center\",\n    verticalAlign: \"bottom\",\n    itemStyle: { color: t.inkSoft, fontSize: \"14px\" },\n    itemHoverStyle: { color: t.ink },\n    symbolRadius: 6,\n  },\n  plotOptions: {\n    series: { animation: false, states: { inactive: { opacity: 1 } } },\n    scatter: { stickyTracking: false },\n  },\n  series,\n});\n"}