// Draws the paper's figures as SVG into figures/, from data/*.json and numerics/*.json.
// Colours and type come from classes defined in paper/page.css, so the drawings follow
// the page theme. Every mark carries a <title> with its value and source.

import fs from 'node:fs'
import path from 'node:path'
import { fileURLToPath } from 'node:url'

const root = path.dirname(fileURLToPath(import.meta.url))
const json = p => JSON.parse(fs.readFileSync(path.join(root, p), 'utf8'))
const out = (name, svg) => fs.writeFileSync(path.join(root, 'figures', `${name}.svg`), svg)

const esc = s => String(s).replace(/&/g, '&amp;').replace(/</g, '&lt;')
const r1 = x => Math.round(x * 10) / 10
const lin = (d0, d1, r0, r1_) => v => r0 + (v - d0) / (d1 - d0) * (r1_ - r0)
const log = (d0, d1, r0, r1_) => v => r0 + (Math.log(v) - Math.log(d0)) / (Math.log(d1) - Math.log(d0)) * (r1_ - r0)
const text = (x, y, s, o = {}) =>
  `<text x="${r1(x)}" y="${r1(y)}"${o.anchor ? ` text-anchor="${o.anchor}"` : ''}${o.cls ? ` class="${o.cls}"` : ''}>${esc(s)}</text>`
const line = (x1, y1, x2, y2, cls) => `<line x1="${r1(x1)}" y1="${r1(y1)}" x2="${r1(x2)}" y2="${r1(y2)}" class="${cls}"/>`
const dot = (x, y, cls, title, r = 4.5) =>
  `<circle cx="${r1(x)}" cy="${r1(y)}" r="${r}" class="${cls} ring"><title>${esc(title)}</title></circle>`
const hollow = (x, y, strokeCls, title, r = 4.5) =>
  `<circle cx="${r1(x)}" cy="${r1(y)}" r="${r}" class="paper-fill ${strokeCls}" stroke-width="2"><title>${esc(title)}</title></circle>`
const frame = (w, h, label, body) =>
  `<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 ${w} ${h}" role="img" aria-label="${esc(label)}">${body}</svg>\n`

// Record Tc against year.
{
  const d = json('data/records.json')
  const W = 720, H = 380, L = 44, R = 150, T = 34, B = 34
  const x = lin(1905, d.end_year + 1, L, W - R), y = lin(0, 320, H - B, T)
  let s = ''
  for (const t of [0, 50, 100, 150, 200, 250, 300]) {
    s += line(L, y(t), W - R, y(t), t === 0 ? 'axis' : 'grid') + text(L - 8, y(t) + 4, t, { anchor: 'end' })
  }
  for (const yr of [1920, 1940, 1960, 1980, 2000, 2020]) s += text(x(yr), H - B + 18, yr, { anchor: 'middle' })
  s += text(2, T - 14, 'Tc (K)')

  const step = (pts, end, startAt) => {
    let p = startAt ? `M${r1(x(startAt.year))},${r1(y(startAt.tc))}` : `M${r1(x(pts[0].year))},${r1(y(0))}`
    let last = startAt ? startAt.tc : 0
    for (const q of pts) { p += `L${r1(x(q.year))},${r1(y(last))}L${r1(x(q.year))},${r1(y(q.tc))}`; last = q.tc }
    return p + `L${r1(x(end))},${r1(y(last))}`
  }
  // the any-pressure series starts from the ambient record standing when it first rose above it
  const ambTop = d.ambient.filter(q => q.year <= d.any[0].year).pop()
  s += `<path d="${step(d.any, d.end_year, ambTop)}" class="line s2"/>`
  s += `<path d="${step(d.ambient, d.end_year)}" class="line s1"/>`
  for (const q of d.ambient.filter(q => q.year >= 1986)) s += dot(x(q.year), y(q.tc), 'f1', `${q.material}: ${q.tc} K, ${q.year}`, 3.5)
  for (const q of d.any) s += dot(x(q.year), y(q.tc), 'f2', `${q.material}: ${q.tc} K, ${q.year}`, 3.5)
  for (const q of d.unreproduced) s += hollow(x(q.year), y(q.tc), q.series === 'ambient' ? 's1' : 's2', `${q.material}: ${q.tc} K, ${q.year}, not independently reproduced`)

  const xr = x(d.end_year) + 10
  s += text(xr, y(250) + 4, 'Any pressure: 250 K', { cls: 't-ink' })
  s += text(xr, y(250) + 19, 'LaH₁₀ at 170 GPa')
  s += text(xr, y(138) + 4, 'Ambient pressure: 138 K', { cls: 't-ink' })
  s += text(xr, y(138) + 19, 'Hg cuprates, since 1993')
  s += text(xr, y(298) + 4, 'LaSc₂H₂₄, not reproduced')
  s += text(xr, y(151) - 15, 'quenched Hg-1223,')
  s += text(xr, y(151) - 1, 'one group')
  s += text(x(1987) - 8, y(93) + 4, 'YBa₂Cu₃O₇', { anchor: 'end' })
  s += text(x(2015) - 8, y(203) - 6, 'H₃S', { anchor: 'end' })
  s += text(x(1994) - 8, y(164) - 6, 'Hg-1223 at 31 GPa', { anchor: 'end' })
  s += text(x(1973), y(22.3) - 10, 'Nb₃Ge', { anchor: 'middle' })

  // legend
  const ly = 14, LX = 90
  s += `<path d="M${LX},${ly}h22" class="line s1"/>` + text(LX + 28, ly + 4, 'Ambient pressure')
  s += `<path d="M${LX + 150},${ly}h22" class="line s2"/>` + text(LX + 178, ly + 4, 'Any pressure')
  s += hollow(LX + 290, ly, 'muted-line', 'single group, not reproduced', 4) + text(LX + 300, ly + 4, 'Single group, not reproduced')
  out('records', frame(W, H, 'Record superconducting transition temperature by year at ambient pressure and at any pressure', s))
}

// Coupling strength against phonon frequency, with Eliashberg isotherms.
{
  const d = json('data/plane.json'), req = json('numerics/requirements.json')
  const W = 720, H = 430, L = 56, R = 64, T = 34, B = 44
  const XMAX = 4.1
  const x = lin(0.5, XMAX, L, W - R), y = log(40, 4000, H - B, T)
  let s = `<defs><clipPath id="plane-clip"><rect x="${L}" y="${T}" width="${W - R - L}" height="${H - B - T}"/></clipPath></defs>`
  for (const t of [50, 100, 200, 500, 1000, 2000]) s += line(L, y(t), W - R, y(t), 'grid') + text(L - 8, y(t) + 4, t, { anchor: 'end' })
  s += line(L, H - B, W - R, H - B, 'axis')
  for (const t of [0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4]) s += line(x(t), H - B, x(t), H - B + 4, 'axis') + text(x(t), H - B + 18, t.toFixed(1), { anchor: 'middle' })
  s += text((L + W - R) / 2, H - 6, 'Electron-phonon coupling λ', { anchor: 'middle' })
  s += text(2, T - 14, 'ω_log (K)')

  const iso = k => req.isotherms[k].filter(p => p.lambda >= 0.5 && p.lambda <= XMAX + 0.01)
  const pathOf = pts => pts.map((p, i) => `${i ? 'L' : 'M'}${r1(x(p.lambda))},${r1(y(p.omega_K))}`).join('')
  const top = iso('300')
  s += `<path d="${pathOf(top)}L${r1(x(top[top.length - 1].lambda))},${T}L${r1(x(top[0].lambda))},${T}Z" class="wash2" clip-path="url(#plane-clip)"/>`
  for (const k of ['39', '100', '200', '300']) {
    const pts = iso(k), end = pts[pts.length - 1]
    s += `<path d="${pathOf(pts)}" class="muted-line" stroke-width="${k === '300' ? 2 : 1.25}" clip-path="url(#plane-clip)"/>`
    s += text(W - R + 6, y(end.omega_K) + 4, `${k} K`, { cls: k === '300' ? 't-ink' : '' })
  }
  s += text(x(3.3), y(2900), 'Tc above 300 K', { anchor: 'middle' })

  const cls = { pressure: 'f2', predicted: 'f1', ambient: 'f3' }
  for (const p of d.points) s += dot(x(p.lambda), y(p.omega_log), cls[p.group], `${p.label}: λ = ${p.lambda}, ω_log = ${p.omega_log} K. ${p.detail}`)
  const at = (lam, om, dx, dy, str, anchor) => text(x(lam) + dx, y(om) + dy, str, { anchor })
  s += at(2.67, 1119, 9, 4, 'LaH₁₀') + at(1.84, 1078, -9, 4, 'H₃S', 'end') + at(2.0, 1253, 0, -10, 'CaH₆, YH₆', 'middle')
  s += `<line x1="${r1(x(2.1))}" y1="${r1(y(580))}" x2="${r1(x(2.56))}" y2="${r1(y(391))}" class="s1" stroke-width="1"/>`
  s += at(2.56, 391, 4, 18, 'Mg₂IrH₆, two calculations', 'middle') + at(1.96, 491, -9, 12, 'Li₂CuH₆', 'end')
  s += at(2.21, 609, 9, -2, 'PdH₄')
  s += at(1.35, 766, 0, -11, 'Mg₂RhH₆, Mg₂PtH₆', 'middle')
  s += at(3.82, 328, 0, 22, 'Li₂AgH₆, Li₂AuH₆', 'middle') + at(3.28, 294, -9, 14, 'B₂C₈Cl', 'end')
  s += at(0.87, 725, -9, 4, 'MgB₂', 'end') + at(1.24, 154, 9, 4, 'Nb') + at(1.45, 58, 9, 4, 'Pb')

  const ly = 14, LX = 110
  s += dot(LX + 5, ly, 'f2', 'measured; exists only under pressure') + text(LX + 15, ly + 4, 'Measured, 150 GPa or more')
  s += dot(LX + 205, ly, 'f1', 'predicted at 1 atm') + text(LX + 215, ly + 4, 'Predicted at 1 atm')
  s += dot(LX + 355, ly, 'f3', 'measured at 1 atm') + text(LX + 365, ly + 4, 'Measured at 1 atm')
  out('plane', frame(W, H, 'Published coupling strength and logarithmic phonon frequency of superconductors against isotherms of transition temperature from the Eliashberg equations', s))
}

// Predicted Tc against energy above the hull.
{
  const d = json('data/frontier.json')
  const W = 720, H = 420, L = 44, R = 24, T = 34, B = 44
  const x = lin(-62, 340, L, W - R), y = lin(0, 190, H - B, T)
  let s = ''
  for (const t of [0, 50, 100, 150]) s += line(L, y(t), W - R, y(t), t === 0 ? 'axis' : 'grid') + text(L - 8, y(t) + 4, t, { anchor: 'end' })
  for (const t of [0, 100, 200, 300]) s += line(x(t), H - B, x(t), H - B + 4, 'axis') + text(x(t), H - B + 18, t, { anchor: 'middle' })
  s += text((x(0) + W - R) / 2, H - 6, 'Energy above the convex hull (meV per atom)', { anchor: 'middle' })
  s += text(2, T - 14, 'Tc (K)')

  s += line(x(d.metastability_median), T + 6, x(d.metastability_median), H - B, 'muted-line')
  s += line(x(d.metastability_p90), T + 6, x(d.metastability_p90), H - B, 'muted-line')
  s += text(x(d.metastability_median) + 5, y(4), 'median')
  s += text(x(d.metastability_p90) + 5, y(4), '90th percentile of known metastable phases')

  // Each compound is drawn at its current hull distance. A thin whisker runs back to the
  // value first published. Points that share a hull distance are nudged apart.
  const nudge = { 'MgB₂': -5, 'LiMoN₂': 5, 'Mg₂PtH₆': -6, 'Mg₂PdH₆': 6, 'Li₂AgH₆': -6, 'PdH₄': 6 }
  const whiskerAt = { 'Mg₂PdH₆': 0, 'KInH₃': 0 }  // 0 = lower end of the Tc range, default upper end
  for (const p of d.points) {
    const cx = x(p.hull[1]) + (nudge[p.label] || 0)
    const c = p.kind === 'measured' ? 3 : 1
    const tcText = p.tc[0] === p.tc[1] ? p.tc[0] : `${p.tc[0]} to ${p.tc[1]}`
    const hullText = p.hull[0] === p.hull[1] ? p.hull[1] : `${p.hull[1]} (first published as ${p.hull[0]})`
    const title = `${p.label}: Tc ${tcText} K, ${hullText} meV per atom above the hull. ${p.status}`
    if (p.hull[0] !== p.hull[1]) {
      const yw = y(p.tc[p.label in whiskerAt ? whiskerAt[p.label] : 1])
      s += `<line x1="${r1(x(p.hull[0]))}" y1="${r1(yw)}" x2="${r1(cx)}" y2="${r1(yw)}" class="muted-line" stroke-width="1"/>`
      s += `<line x1="${r1(x(p.hull[0]))}" y1="${r1(yw - 4)}" x2="${r1(x(p.hull[0]))}" y2="${r1(yw + 4)}" class="muted-line" stroke-width="1"/>`
    }
    if (p.tc[0] !== p.tc[1]) s += `<line x1="${r1(cx)}" y1="${r1(y(p.tc[0]))}" x2="${r1(cx)}" y2="${r1(y(p.tc[1]))}" class="s${c}" stroke-width="2" stroke-linecap="round"><title>${esc(title)}</title></line>`
    s += dot(cx, y(p.tc[1]), `f${c}`, title, 4)
    if (p.tc[0] !== p.tc[1]) s += dot(cx, y(p.tc[0]), `f${c}`, title, 4)
  }
  const at = (h, t, dx, dy, str, anchor, cls) => text(x(h) + dx, y(t) + dy, str, { anchor, cls })
  s += at(0, 39, -13, 6, 'MgB₂, measured', 'end', 't-ink') + at(6, 46, -6, 0, 'LiMoN₂', 'end') + at(0, 22.9, -10, 4, 'Mo₂Pr₃N₆', 'end')
  s += at(67, 45, 6, 12, 'Mg₂RhH₆', null, 't-ink') + at(67, 45, 6, 26, 'lost on pressure release')
  s += at(85.5, 175, 0, -22, 'Mg₂IrH₆', 'end', 't-ink') + at(85.5, 175, 0, -8, 'not formed', 'end')
  s += at(114.6, 73, 9, -8, 'KInH₃')
  s += at(195.7, 80, 16, 0, 'Mg₂PtH₆', null, 't-ink') + at(195.7, 80, 16, 14, 'Mg₄Pt₃H₆ formed')
  s += at(195.7, 51, 16, 8, 'Mg₂PdH₆')
  s += at(171.5, 140, -9, 0, 'Li₂AuH₆', 'end', 't-ink') + at(171.5, 140, -9, 14, 'unstable in dynamics', 'end')
  s += at(319.1, 108, -16, 0, 'Li₂AgH₆', 'end', 't-ink') + at(319.1, 108, -16, 14, 'collapses in dynamics', 'end')
  s += at(318.9, 146, 10, -22, 'PdH₄', 'end', 't-ink') + at(318.9, 146, 10, -8, 'harmonic only', 'end')

  const ly = 14, LX = 90
  s += dot(LX + 5, ly, 'f1', 'calculated') + text(LX + 15, ly + 4, 'Calculated, published range')
  s += dot(LX + 205, ly, 'f3', 'measured') + text(LX + 215, ly + 4, 'Measured')
  s += `<line x1="${LX + 290}" y1="${ly}" x2="${LX + 312}" y2="${ly}" class="muted-line" stroke-width="1"/>` + text(LX + 318, ly + 4, 'Hull distance when first published')
  out('frontier', frame(W, H, 'Calculated transition temperature against energy above the convex hull for ambient-pressure candidates', s))
}

// Survival function of calculated Tc by distance from the hull.
{
  const d = json('numerics/tails.json')
  const W = 720, H = 400, L = 56, R = 24, T = 34, B = 44
  const x = lin(0, 118, L, W - R), y = log(1e-4, 1, H - B, T)
  let s = `<defs><clipPath id="tail-clip"><rect x="${L}" y="${T}" width="${W - R - L}" height="${H - B - T}"/></clipPath></defs>`
  for (const [v, lab] of [[1, '1'], [0.1, '0.1'], [0.01, '0.01'], [0.001, '0.001'], [0.0001, '0.0001']]) s += line(L, y(v), W - R, y(v), v === 1e-4 ? 'axis' : 'grid') + text(L - 8, y(v) + 4, lab, { anchor: 'end' })
  for (const t of [0, 20, 40, 60, 80, 100]) s += line(x(t), H - B, x(t), H - B + 4, 'axis') + text(x(t), H - B + 18, t, { anchor: 'middle' })
  s += text((L + W - R) / 2, H - 6, 'Calculated Tc (K)', { anchor: 'middle' })
  s += text(2, T - 14, 'Fraction at or above')

  const grid = d.survival.grid
  const merge = (a, na, b, nb) => a.map((v, i) => (v * na + b[i] * nb) / (na + nb))
  const st = Object.fromEntries(d.strata.map(q => [`${q.lo}-${q.hi}`, q]))
  const near = merge(d.survival['0-25'], st['0-25'].n, d.survival['25-50'], st['25-50'].n)
  const series = [
    { key: 'ord1', vals: near, label: '0 to 50' },
    { key: 'ord2', vals: d.survival['50-100'], label: '50 to 100' },
    { key: 'ord3', vals: d.survival['100-200'], label: '100 to 200' },
    { key: 'ord4', vals: d.survival['200-376'], label: '200 to 376' },
  ]
  // Exponential fit to the first curve: solid over the range it describes (2 to 30 K),
  // dashed beyond, where the two highest compounds lie above it.
  const f = d.near_hull_fit
  const fit = t => f.p_u * Math.exp(-(t - f.u) / f.tau)
  s += `<path d="M${r1(x(2))},${r1(y(fit(2)))}L${r1(x(30))},${r1(y(fit(30)))}" class="ink-line" stroke-width="1" clip-path="url(#tail-clip)"/>`
  s += `<path d="M${r1(x(30))},${r1(y(fit(30)))}L${r1(x(60))},${r1(y(fit(60)))}" class="ink-line" stroke-width="1" stroke-dasharray="4 3" clip-path="url(#tail-clip)"/>`
  for (const q of series.slice().reverse()) {
    let p = ''
    grid.forEach((g, i) => { if (q.vals[i] > 0) p += `${p ? 'L' : 'M'}${r1(x(g))},${r1(y(Math.max(q.vals[i], 1e-4)))}` })
    s += `<path d="${p}" class="line ${q.key}" clip-path="url(#tail-clip)"><title>${esc(`${q.label} meV per atom above the hull`)}</title></path>`
  }
  // The two compounds above the fit.
  const top = d.pareto_front.filter(q => q.e_hull <= 50).sort((a, b) => b.tc - a.tc)
  s += text(x(top[1].tc) + 5, y(2 / f.n) - 5, `LiMoN₂, ${top[1].tc.toFixed(1)} K`)
  s += text(x(top[0].tc) + 6, y(1 / f.n) + 1, `Mg₂RhH₆, ${top[0].tc.toFixed(1)} K`)

  // Key, in the empty corner above the curves.
  const kx = x(66), ky = T + 12
  s += text(kx, ky + 4, 'Energy above the convex hull (meV per atom)', { cls: 't-ink' })
  series.forEach((q, i) => { s += `<path d="M${r1(kx)},${ky + 18 * (i + 1)}h22" class="line ${q.key}"/>` + text(kx + 28, ky + 18 * (i + 1) + 4, q.label) })
  s += `<path d="M${r1(kx)},${ky + 90}h22" class="ink-line" stroke-width="1"/>` + text(kx + 28, ky + 94, `Exponential fit to 0 to 50, scale ${f.tau.toFixed(1)} K;`)
  s += text(kx + 28, ky + 108, 'dashed beyond 30 K')
  out('tail', frame(W, H, 'Fraction of compounds with calculated transition temperature at or above a given value, in four ranges of energy above the convex hull', s))
}

// Hopfield parameter of hydrogen against hydrogen number density: the release and the megabar hydrides.
{
  const pts = json('numerics/h_alexandria_points.json')
  const stat = json('numerics/h_alexandria.json')
  const mb = json('data/h_census.json').rows.filter(r => r.include_in_stats && r.group === 'megabar')
  const W = 720, H = 460, L = 52, R = 96, T = 34, B = 44
  const x = log(0.002, 0.6, L, W - R), y = log(0.01, 20, H - B, T)
  const K = 136.5
  let s = `<defs><clipPath id="pack-clip"><rect x="${L}" y="${T}" width="${W - R - L}" height="${H - B - T}"/></clipPath></defs>`
  for (const t of [0.01, 0.1, 1, 10]) s += line(L, y(t), W - R, y(t), t === 0.01 ? 'axis' : 'grid') + text(L - 8, y(t) + 4, t, { anchor: 'end' })
  for (const t of [0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5]) s += line(x(t), H - B, x(t), H - B + 4, 'axis') + text(x(t), H - B + 18, t, { anchor: 'middle' })
  s += text((L + W - R) / 2, H - 6, 'Hydrogen number density (atoms per ų)', { anchor: 'middle' })
  s += text(2, T - 14, 'Hopfield parameter of hydrogen (eV/Ų)')

  s += `<rect x="${L}" y="${r1(y(12))}" width="${W - R - L}" height="${r1(y(8.7) - y(12))}" class="wash2"/>`
  s += text(L + 8, y(10.2) + 4, 'needed for 300 K, single mode')
  for (const t of [50, 100, 200, 300, 400, 500]) { const sv = (t / K) ** 2; s += line(W - R, y(sv), W - R + 4, y(sv), 'axis') + text(W - R + 8, y(sv) + 4, `${t} K`) }
  s += text(W - R + 8, T - 14, 'asymptote')

  const hl = (h, cls, w) => `<path d="M${r1(x(0.002))},${r1(y(0.002 * h))}L${r1(x(0.6))},${r1(y(0.6 * h))}" class="${cls}" stroke-width="${w}" clip-path="url(#pack-clip)"/>`
  const gm = stat.h.geometric_mean
  s += hl(36.6, 'ink-line', 1.25) + hl(gm, 'muted-line', 1)
  const kx = x(0.105), ky = y(0.075)
  s += `<line x1="${r1(kx)}" y1="${r1(ky)}" x2="${r1(kx + 22)}" y2="${r1(ky)}" class="ink-line" stroke-width="1.25"/>` + text(kx + 28, ky + 4, '36.6 eV Å per proton:', { cls: 't-ink' }) + text(kx + 28, ky + 18, 'one proton in an electron gas', { cls: 't-ink' })
  s += `<line x1="${r1(kx)}" y1="${r1(ky + 36)}" x2="${r1(kx + 22)}" y2="${r1(ky + 36)}" class="muted-line" stroke-width="1"/>` + text(kx + 28, ky + 40, `${gm.toFixed(1)} eV Å per proton:`) + text(kx + 28, ky + 54, 'geometric mean at 1 atm')

  s += line(x(0.090), y(7.5), x(0.090), H - B, 'muted-line')
  s += text(x(0.090) + 6, H - B - 22, 'densest hydrogen stable')
  s += text(x(0.090) + 6, H - B - 8, 'at 1 atm (TiH₂, Mg₂FeH₆)')

  // the release: one small mark per compound, no tooltip (4,619 of them)
  let dots = ''
  for (let i = 0; i < pts.rho.length; i++) dots += `<circle cx="${r1(x(pts.rho[i]))}" cy="${r1(y(Math.max(pts.eta[i], 0.0101)))}" r="1.5"/>`
  s += `<g class="f1" fill-opacity="0.45" clip-path="url(#pack-clip)">${dots}</g>`

  const val = r => r.eta_H_eV_A2 ?? r.S_eV_A2
  for (const name of [...new Set(mb.map(r => r.compound))]) {
    const q = mb.filter(r => r.compound === name).sort((a, b) => a.pressure_GPa - b.pressure_GPa)
    if (q.length > 1) s += `<path d="${q.map((r, i) => `${i ? 'L' : 'M'}${r1(x(r.rho_H))},${r1(y(val(r)))}`).join('')}" class="s2" fill="none" stroke-width="1.25"/>`
  }
  for (const r of mb) s += dot(x(r.rho_H), y(val(r)), 'f2', `${r.compound} at ${r.pressure_GPa} GPa (calculated): ${r.rho_H.toFixed(3)} H per ų, ${val(r).toFixed(1)} eV/Ų, h = ${r.h_eV_A.toFixed(0)} eV Å`, 4)

  // the five largest values at one atmosphere, drawn over the cloud and named
  const top = stat.largest_eta_H.slice(0, 5)
  for (const r of top) s += dot(x(r.rho_H), y(r.eta_H), 'f1', `${r.formula}, 1 atm (calculated): ${r.rho_H} H per ų, ${r.eta_H} eV/Ų, h = ${r.h} eV Å`, 3.5)
  s += text(x(0.064) - 10, y(5.5) + 4, 'KPtH₆, RbPtH₆, CsPtH₆,', { anchor: 'end' }) + text(x(0.064) - 10, y(5.5) + 18, 'RbNiH₆, CsNiH₆', { anchor: 'end' })

  const find = (name, p) => mb.find(r => r.compound === name && (p === undefined || r.pressure_GPa === p))
  const lab = (r, dx, dy, str, anchor) => r ? text(x(r.rho_H) + dx, y(val(r)) + dy, str, { anchor }) : ''
  s += lab(find('SH3', 220), -8, -8, 'H₃S', 'end') + lab(find('CaH6', 150), -8, 14, 'CaH₆', 'end') + lab(find('LaH10', 300), 7, -7, 'LaH₁₀')
  s += lab(find('MgH6', 300), -8, -5, 'MgH₆', 'end') + lab(find('YH10', 400), -4, -10, 'YH₁₀', 'middle')

  const ly = 14, LX = 250
  s += `<circle cx="${LX + 5}" cy="${ly}" r="2.5" class="f1"/>` + text(LX + 15, ly + 4, `Calculated at 1 atm (${pts.rho.length.toLocaleString('en-US')})`)
  s += dot(LX + 215, ly, 'f2', 'calculated at megabar pressure', 4) + text(LX + 225, ly + 4, 'Calculated, megabar')
  out('packing', frame(W, H, `Hopfield parameter of hydrogen against hydrogen number density for ${pts.rho.length} hydrides calculated at one atmosphere and five superhydrides at megabar pressure`, s))
}

// Distribution of the scattering strength per proton at one atmosphere.
{
  const pts = json('numerics/h_alexandria_points.json')
  const mb = json('data/h_census.json').rows.filter(r => r.include_in_stats && r.group === 'megabar')
  const W = 720, H = 360, L = 52, R = 24, T = 44, B = 44
  const per = 8, lo = 0.4, n = Math.ceil(Math.log10(120 / lo) * per)
  const bin = h => Math.floor(Math.log10(h / lo) * per)
  const all = new Array(n).fill(0), hot = new Array(n).fill(0)
  for (let i = 0; i < pts.rho.length; i++) {
    const b = bin(pts.eta[i] / pts.rho[i])
    if (b >= 0 && b < n) { all[b]++; if (pts.tc20[i]) hot[b]++ }
  }
  const ymax = Math.ceil(Math.max(...all) / 100) * 100
  const x = log(lo, lo * 10 ** (n / per), L, W - R), y = lin(0, ymax, H - B, T)
  let s = ''
  for (let t = 0; t <= ymax; t += 200) s += line(L, y(t), W - R, y(t), t === 0 ? 'axis' : 'grid') + text(L - 8, y(t) + 4, t, { anchor: 'end' })
  for (const t of [0.5, 1, 2, 5, 10, 20, 50, 100]) s += line(x(t), H - B, x(t), H - B + 4, 'axis') + text(x(t), H - B + 18, t, { anchor: 'middle' })
  s += text((L + W - R) / 2, H - 6, 'Scattering strength per proton, h (eV Å)', { anchor: 'middle' })
  s += text(2, T - 24, 'Compounds')
  const edge = b => lo * 10 ** (b / per)
  for (let b = 0; b < n; b++) {
    const x0 = x(edge(b)) + 1, w = x(edge(b + 1)) - x(edge(b)) - 2
    if (all[b]) s += `<rect x="${r1(x0)}" y="${r1(y(all[b]))}" width="${r1(w)}" height="${r1(y(0) - y(all[b]))}" class="fo4"><title>${edge(b).toFixed(1)} to ${edge(b + 1).toFixed(1)} eV Å: ${all[b]} compounds, ${hot[b]} of them calculated above 20 K</title></rect>`
    if (hot[b]) s += `<rect x="${r1(x0)}" y="${r1(y(hot[b]))}" width="${r1(w)}" height="${r1(y(0) - y(hot[b]))}" class="fo1"/>`
  }
  s += line(x(36.6), T - 6, x(36.6), H - B, 'ink-line')
  s += text(x(36.6) + 6, T + 34, 'one proton in an', { cls: 't-ink' }) + text(x(36.6) + 6, T + 48, 'electron gas, 36.6', { cls: 't-ink' })
  // megabar compounds, one tick per pressure point
  for (const r of mb) s += `<line x1="${r1(x(r.h_eV_A))}" y1="${T + 2}" x2="${r1(x(r.h_eV_A))}" y2="${T + 14}" class="s2" stroke-width="2"><title>${esc(r.compound)} at ${r.pressure_GPa} GPa: h = ${r.h_eV_A.toFixed(1)} eV Å</title></line>`
  s += text(x(23) - 6, T + 12, 'megabar hydrides', { anchor: 'end' })
  const ly = T - 24, LX = 150
  s += `<rect x="${LX}" y="${ly - 6}" width="12" height="12" class="fo4"/>` + text(LX + 18, ly + 4, `All ${pts.rho.length.toLocaleString('en-US')} labelled hydrides at 1 atm`)
  s += `<rect x="${LX + 270}" y="${ly - 6}" width="12" height="12" class="fo1"/>` + text(LX + 288, ly + 4, `Calculated Tc above 20 K (${pts.tc20.reduce((a, b) => a + b, 0)})`)
  out('hdist', frame(W, H, 'Distribution of the scattering strength per proton among hydrides calculated at one atmosphere, with the subset calculated above 20 K', s))
}

console.log('figures written')