#!/usr/bin/env python3 """Vasilyev et al. (2024)'s Table 1 against their own Figs. 3, 4 and S6: the readings on Fact-Checked Science's page https://factchecked.science/tess-superflare-rate/ (section 6). Our analysis. Python 3 with pypdf (pip install pypdf); nothing else. Run it with `python3 -I`: python3 -I vasilyev_figures.py 2412.12265v1.pdf [output.json] The PDF is the accepted version of Vasilyev et al., "Sun-like stars produce superflares roughly once per century", Science 386, 1301-1305 (2024), doi:10.1126/science.adl5441, at https://arxiv.org/pdf/2412.12265v1 (40 pages: the paper, then its supplementary materials). Our copy's SHA-256 was 48462741cd7159f7636e53042663d6e87291d59dedc3016b13c6ffcc8872358a; a different copy is read anyway, with a warning. Its Figs. 3, 4 and S6 (PDF pages 11, 12 and 35) are vector drawings: the script reads their paths and markers with pypdf, as data, and calibrates each axis from the figure's own tick marks. It prints, and writes as JSON, the numbers the page quotes: the integrals and their bounds, the applied completeness bin by bin from 10^34.4 erg up, and Fig. S6's detection fractions from 10^34.5 erg up (42 marker values), but not the figures' curves. These are our readings of figures for which Crossref lists no licence, offered as facts. What it reads Fig. 3 (p. 11): the flare frequency per star per year per erg of the full sample (thick blue histogram) and the same corrected for the detection threshold (thin blue), and the same for the tighter cuts (orange). Integrated over 1e34-1e35 erg and above 1e35 erg, taking each bin as flat, for comparison with Table 1 (p. 10): 8.63e-3 and 3.90e-4 per star per year for the full sample, 9.31e-3 and 3.59e-4 for the tighter cuts. Fig. 4 (p. 12): the cumulative frequencies; each step drawn over a bin is the rate above that bin's upper edge (its top step equals Fig. 3's top bin alone), so the rate above 1e34 erg is read between the two steps that straddle it. Fig. S6 (supplement p. 17, PDF p. 35): the detected fraction r of injected flares at 17 energies (5e32 to 5e35 erg, log-spaced, in the colour order of the figure) and six temperatures; compared with the supplement's "Flares with energies E > 5 * 10^34 erg are detected on all stars (r = 1)" (p. 11) and with the completeness the correction in Fig. 3 applied (observed over corrected, bin by bin). """ import hashlib import json import sys from pypdf import PdfReader from pypdf.generic import ContentStream SHA = "48462741cd7159f7636e53042663d6e87291d59dedc3016b13c6ffcc8872358a" STAR_YEARS = 220000.0 # "our full sample corresponds to ≈ 220,000 years of stellar activity" (p. 4) def ops(form, reader): for operands, op in ContentStream(form, reader).operations: yield operands, (op.decode("latin-1") if isinstance(op, bytes) else op) def strokes(reader, page, name): """Stroked subpaths of a Form XObject: colour, width and vertices (translations applied).""" form = reader.pages[page - 1]["/Resources"]["/XObject"][name].get_object() ctm, stack, rg, w, sub, cur, out = [1, 0, 0, 1, 0, 0], [], None, None, [], [], [] def ap(x, y): a, b, c, d, e, f = ctm return (a * x + c * y + e, b * x + d * y + f) for operands, op in ops(form, reader): if op == "q": stack.append((list(ctm), rg, w)) elif op == "Q": ctm, rg, w = stack.pop() elif op == "cm": a2, b2, c2, d2, e2, f2 = (float(v) for v in operands) a, b, c, d, e, f = ctm ctm = [a * a2 + b * c2, a * b2 + b * d2, c * a2 + d * c2, c * b2 + d * d2, e * a2 + f * c2 + e2, e * b2 + f * d2 + f2] elif op == "RG": rg = tuple(round(float(v), 4) for v in operands) elif op == "G": rg = (round(float(operands[0]), 4),) * 3 elif op == "w": w = float(operands[0]) elif op == "m": if sub: cur.append(sub) sub = [ap(float(operands[0]), float(operands[1]))] elif op == "l": sub.append(ap(float(operands[0]), float(operands[1]))) elif op in ("S", "s", "B", "b", "f", "n"): if sub: cur.append(sub) if op in ("S", "s", "B", "b"): out += [(rg, w, p) for p in cur] cur, sub = [], [] return out def ticks(paths, length): xs, ys = set(), set() for rg, _w, p in paths: if rg == (0.0, 0.0, 0.0) and len(p) == 2: (xa, ya), (xb, yb) = p if abs(xa - xb) < 1e-6 and abs(abs(ya - yb) - length) < 1e-6: xs.add(round(xa, 4)) if abs(ya - yb) < 1e-6 and abs(abs(xa - xb) - length) < 1e-6: ys.add(round(ya, 4)) return sorted(xs), sorted(ys) def steps(p, X, Y): return [(X(p[i][0]), X(p[i + 1][0]), Y(p[i][1])) for i in range(len(p) - 1) if abs(p[i][1] - p[i + 1][1]) < 1e-6 and p[i + 1][0] > p[i][0]] def integral(bins, lo, hi): return sum(10 ** lf * (10 ** min(b, hi) - 10 ** max(a, lo)) for a, b, lf in bins if min(b, hi) > max(a, lo)) def main(): pdf = sys.argv[1] out_path = sys.argv[2] if len(sys.argv) > 2 else "vasilyev-figures.json" sha = hashlib.sha256(open(pdf, "rb").read()).hexdigest() if sha != SHA: print(f"warning: {pdf} has SHA-256 {sha}, not our copy's {SHA}; reading it anyway") reader = PdfReader(pdf) res = {"pdf_sha256": sha} # Fig. 3: four major ticks on each axis; x 1e33..1e36 erg, y 1e-41..1e-35 per star per year per erg f3 = strokes(reader, 11, "/Im3") xs, ys = ticks(f3, 5.0) X = lambda x: 33 + 3 * (x - xs[0]) / (xs[-1] - xs[0]) yt = ys # four major ticks, 1e-41 to 1e-35, two decades apart if len(xs) != 4 or len(yt) != 4: sys.exit(f"Fig. 3: expected 4 major ticks on each axis, found {xs} and {yt}") Y = lambda y: -41 + 2 * (y - yt[0]) / ((yt[-1] - yt[0]) / 3) curves = {} names = {(0.0, 0.4471, 0.698): "full sample", (0.3373, 0.7059, 0.9137): "full sample, corrected", (0.8353, 0.3686, 0.0): "tighter cuts", (0.949, 0.5569, 0.1686): "tighter cuts, corrected"} for rg, _w, p in f3: if rg in names and len(p) > 20: curves[names[rg]] = steps(p, X, Y) res["fig3_ticks"] = {"x": xs, "y": yt} res["fig3"] = {} for k, bins in curves.items(): inner = sum(10 ** lf * (10 ** b - 10 ** a) for a, b, lf in bins if a >= 34 and b <= 35) edges = sum(10 ** lf * (10 ** b - 10 ** a) for a, b, lf in bins if a < 34 < b or a < 35 < b) above = sum(10 ** lf * (10 ** b - 10 ** a) for a, b, lf in bins if a >= 35) edge35 = sum(10 ** lf * (10 ** b - 10 ** a) for a, b, lf in bins if a < 35 < b) res["fig3"][k] = {"rate_1e34_to_1e35": float(f"{integral(bins, 34, 35):.4g}"), "rate_above_1e35": float(f"{integral(bins, 35, 40):.4g}"), "rate_above_1e34": float(f"{integral(bins, 34, 40):.4g}"), # bounds: the bins that straddle 1e34 and 1e35 erg left out, or counted in full "rate_1e34_to_1e35_bounds": [float(f"{inner:.3g}"), float(f"{inner + edges:.3g}")], "rate_above_1e35_bounds": [float(f"{above:.3g}"), float(f"{above + edge35:.3g}")]} full, corr = curves["full sample"], curves["full sample, corrected"] res["fig3"]["full sample"]["total_rate"] = float(f"{integral(full, 0, 99):.4g}") res["fig3"]["full sample"]["flares_at_220000_star_years"] = round(integral(full, 0, 99) * STAR_YEARS) res["fig3"]["full sample"]["share_of_2889_detected"] = round(integral(full, 0, 99) * STAR_YEARS / 2889, 3) obs = {round(a, 2): lf for a, b, lf in full} res["fig3_applied_completeness"] = [{"bin_log_erg": [round(a, 3), round(b, 3)], "observed_over_corrected": round(10 ** (obs[round(a, 2)] - lf), 3)} for a, b, lf in corr if round(a, 2) in obs and a >= 34.4] res["table1"] = {"full sample": {"rate_1e34_to_1e35": 8.63e-3, "rate_above_1e35": 3.90e-4}, "tighter cuts": {"rate_1e34_to_1e35": 9.31e-3, "rate_above_1e35": 3.59e-4}} # Fig. 4: x ticks 1e29..1e36 erg, y ticks 1e-6..1e4 per star per year (two decades per tick) f4 = strokes(reader, 12, "/Im4") xs4, ys4 = ticks(f4, 5.0) X4 = lambda x: 29 + 7 * (x - xs4[0]) / (xs4[-1] - xs4[0]) Y4 = lambda y: -6 + 10 * (y - ys4[0]) / (ys4[-1] - ys4[0]) res["fig4"] = {} for rg, _w, p in f4: label = {(0.0, 0.4471, 0.698): "full sample", (0.3373, 0.7059, 0.9137): "full sample, corrected"}.get(rg) if label and len(p) > 20: st = steps(p, X4, Y4) above = [(b, lf) for a, b, lf in st] # the rate above each step's upper edge lo = max((b, lf) for b, lf in above if b <= 34) hi = min((b, lf) for b, lf in above if b > 34) lf34 = lo[1] + (34 - lo[0]) / (hi[0] - lo[0]) * (hi[1] - lo[1]) # each step of Fig. 4 against the sum of Fig. 3's bins above the step's upper edge f3c = curves[label] diffs = [abs(10 ** lf - sum(10 ** l3 * (10 ** b3 - 10 ** a3) for a3, b3, l3 in f3c if a3 >= b - 1e-3)) for a, b, lf in st if any(a3 >= b - 1e-3 for a3, _b3, _l3 in f3c)] res["fig4"][label] = {"rate_above_1e34": float(f"{10 ** lf34:.4g}"), "max_difference_from_fig3_cumulated": float(f"{max(diffs):.2g}")} # Fig. S6: markers of six temperature columns; x ticks every 200 K from 5000 K, y ticks every 0.2 from r = 0 form = reader.pages[34]["/Resources"]["/XObject"]["/Im11"].get_object() ctm, stack, marks = [1, 0, 0, 1, 0, 0], [], [] for operands, op in ops(form, reader): if op == "q": stack.append(list(ctm)) elif op == "Q": ctm = stack.pop() elif op == "cm": a, _b, _c, d, e, f = (float(v) for v in operands) ctm = [ctm[0] * a, 0, 0, ctm[3] * d, ctm[4] + e * ctm[0], ctm[5] + f * ctm[3]] elif op == "Do" and str(operands[0]).startswith("/P"): marks.append((ctm[4], ctm[5])) s6 = strokes(reader, 35, "/Im11") xs6, _ = ticks(s6, 7.5) ys6 = sorted({round(p[0][1], 4) for rg, _w, p in s6 if len(p) == 2 and abs(p[0][1] - p[1][1]) < 1e-6 and abs(abs(p[0][0] - p[1][0]) - 7.5) < 1e-6}) TX = lambda x: 5000 + 200 * (x - xs6[0]) / ((xs6[-1] - xs6[0]) / (len(xs6) - 1)) RY = lambda y: 0.2 * (y - ys6[0]) / ((ys6[-1] - ys6[0]) / (len(ys6) - 1)) cols = {} for x, y in marks: cols.setdefault(round(TX(x)), []).append(round(RY(y), 4)) loge = [round(32.699 + 0.1875 * k, 3) for k in range(17)] res["figS6"] = {"text_p11": "Flares with energies E > 5 · 10^34 erg are detected on all stars (r = 1)", "r_by_teff": {str(t): {f"{loge[k]}": r for k, r in enumerate(v) if loge[k] >= 34.5} for t, v in sorted(cols.items())}} with open(out_path, "w", encoding="utf-8", newline="\n") as fh: json.dump(res, fh, indent=1) fh.write("\n") for k, v in res["fig3"].items(): print(f"Fig. 3, {k}: {v['rate_1e34_to_1e35']} at 1e34-1e35 erg, {v['rate_above_1e35']} above 1e35 erg") print("Table 1, full sample: 0.00863 and 0.00039; tighter cuts: 0.00931 and 0.000359") print(f"Fig. 3 uncorrected total x 220,000 star-years: {res['fig3']['full sample']['flares_at_220000_star_years']} flares " f"({res['fig3']['full sample']['share_of_2889_detected']} of 2,889)") for k, v in res["fig4"].items(): print(f"Fig. 4, {k}: {v['rate_above_1e34']} above 1e34 erg; differs from Fig. 3 cumulated by at most {v['max_difference_from_fig3_cumulated']}") for k, v in res["fig3"].items(): print(f"Fig. 3, {k}: bounds {v['rate_1e34_to_1e35_bounds']} at 1e34-1e35 erg, {v['rate_above_1e35_bounds']} above") print("Fig. 3 applied completeness:", [(a["bin_log_erg"], a["observed_over_corrected"]) for a in res["fig3_applied_completeness"]]) for t, v in res["figS6"]["r_by_teff"].items(): print(f"Fig. S6, {t} K:", v) print(f"wrote {out_path}") if __name__ == "__main__": main()