#!/usr/bin/env python3
"""Recreate an approximate redraw of NASA/TP-2000-210387 Figure 2.

The NASA report contains plotted curves but no tabulated coordinates.  The
values below are manual visual readings from Fig. 2 (printed page 21), not
source data and not a reimplementation of the cycle equations.  Re-running
this script reproduces the CSV and two responsive SVG plots from those readings.
"""

from __future__ import annotations

import csv
import hashlib
import json
from pathlib import Path


OUT = Path(__file__).resolve().parent
SOURCE_FIGURE = "NASA/TP-2000-210387, Figure 2, printed p. 21 (PDF p. 27)"
POINTS = {
    0.00: [
        (1.00, .46), (1.25, .51), (1.50, .67), (1.75, .81), (2.00, .92),
        (2.25, 1.01), (2.50, 1.06), (2.75, 1.09), (3.00, 1.10),
        (3.25, 1.10), (3.50, 1.09), (3.75, 1.06), (4.00, 1.03),
        (4.25, .99), (4.50, .94), (4.75, .89), (5.00, .84), (5.25, .80),
        (5.50, .75), (5.75, .64), (6.00, .53), (6.25, .49), (6.50, .42),
        (6.75, .35), (7.00, .27), (7.25, .20), (7.50, .12), (7.75, .04),
        (7.90, 0.00),
    ],
    .25: [
        (1.75, .16), (2.00, .37), (2.25, .55), (2.50, .69), (2.75, .80),
        (3.00, .88), (3.25, .93), (3.50, .96), (3.75, .97), (4.00, .97),
        (4.25, .95), (4.50, .94), (4.75, .91), (5.00, .87), (5.25, .83),
        (5.50, .79), (5.75, .73), (6.00, .58), (6.25, .55), (6.50, .53),
        (6.75, .53), (7.00, .50), (7.25, .44), (7.50, .39), (7.75, .34),
        (8.00, .28), (8.25, .22), (8.50, .15), (8.75, .09), (9.00, .04),
        (9.10, 0.00),
    ],
    .50: [
        (2.90, 0.00), (3.00, .20), (3.25, .38), (3.50, .51), (3.75, .59),
        (4.00, .66), (4.25, .72), (4.50, .75), (4.75, .77), (5.00, .78),
        (5.25, .77), (5.50, .76), (5.75, .74), (6.00, .71), (6.25, .71),
        (6.50, .69), (6.75, .66), (7.00, .62), (7.25, .60), (7.50, .56),
        (7.75, .53), (8.00, .49), (8.25, .45), (8.50, .41), (8.75, .38),
        (9.00, .34), (9.25, .30), (9.50, .25), (9.75, .21), (10.00, .17),
        (10.25, .13), (10.50, .10), (10.75, .05), (11.00, .02), (11.10, 0.00),
    ],
}
SERIES = {
    0.00: ("#f2098c", "solid", "No MHD extraction"),
    .25: ("#80d9ff", "10 8", "Extraction ratio 0.25"),
    .50: ("#ffc857", "2 7", "Extraction ratio 0.50"),
}
REVIEW_TRACE = OUT / "nasa-mhd-cycle-fig2-independent-review.csv"
DIRECT_TRACE = OUT / "figure2-direct-raster-trace.csv"
if not DIRECT_TRACE.exists():
    project_trace = OUT.parents[4] / "analyses/nasa-2000-independent-blind/raster-adjudication/figure2-direct-raster-trace.csv"
    if not project_trace.exists():
        raise FileNotFoundError("Place figure2-direct-raster-trace.csv beside this script, or run it from the complete DMGB project checkout.")
    import shutil
    shutil.copyfile(project_trace, DIRECT_TRACE)
REVIEWED_SEGMENTS = []

with DIRECT_TRACE.open(newline="", encoding="utf-8-sig") as f:
    for row in csv.DictReader(f):
        if row["identity_status"] == "clear":
            REVIEWED_SEGMENTS.append((
                float(row["mach_from_pixel"]),
                float(row["ordinate_from_pixel_kN_s_per_kg"]),
                int(row["crop_pixel_x"]),
                int(row["crop_pixel_y"]),
                max(abs(float(row["ordinate_from_pixel_kN_s_per_kg"]) - float(row["ordinate_band_from_pixel_min"])),
                    abs(float(row["ordinate_band_from_pixel_max"]) - float(row["ordinate_from_pixel_kN_s_per_kg"]))),
            ))
assert len(REVIEWED_SEGMENTS) == 23, "Direct raster trace changed; re-audit before rebuilding."


def write_csv() -> None:
    with (OUT / "nasa-mhd-cycle-fig2.csv").open("w", newline="", encoding="utf-8") as f:
        writer = csv.writer(f)
        writer.writerow([
            "flight_mach", "generator_enthalpy_extraction_ratio",
            "specific_thrust_kN_s_per_kg", "visual_reading_uncertainty_kN_s_per_kg",
            "source_figure",
        ])
        for extraction, points in POINTS.items():
            for mach, thrust in points:
                uncertainty = .05 if 4.25 <= mach <= 5.75 else .03
                writer.writerow([f"{mach:.2f}", f"{extraction:.2f}", f"{thrust:.2f}", f"±{uncertainty:.2f}", SOURCE_FIGURE])


def write_methods() -> None:
    methods = {
        "title": "Approximate redraw of NASA/TP-2000-210387 Figure 2",
        "source": {
            "title": "Thermodynamic Cycle Analysis of Magnetohydrodynamic-Bypass Hypersonic Airbreathing Engines",
            "report": "NASA/TP-2000-210387",
            "ntrs_record": "https://ntrs.nasa.gov/citations/20000095934",
            "pdf": "https://ntrs.nasa.gov/api/citations/20000095934/downloads/20000095934.pdf",
            "figure": "Figure 2, printed p. 21 (PDF p. 27)",
            "local_primary_source_sha256": "12867527B3417F5F99A39C2BAE26B663A2D297157C565CCF3F6D4A1402F41BDF",
        },
        "status": "Approximate visual digitization of the published plot. The report publishes no numeric point table for these curves. These values are not raw source data and are not a recalculation of the cycle.",
        "digitization": {
            "method": "The plotted curves were visually read against the printed Mach and specific-thrust axes and sampled at convenient intervals. Rounded points are connected by straight segments for this redraw.",
            "estimated_reading_uncertainty": "About ±0.03 kN·s/kg for clear portions of the scan and about ±0.05 kN·s/kg near overlapping curves; these are reading estimates, not statistical confidence intervals.",
            "points": sum(len(v) for v in POINTS.values()),
            "source_render": {
                "renderer": "Poppler pdftoppm",
                "page_number_pdf": 27,
                "page_render_dpi": 150,
                "page_render_dimensions_pixels": [1271, 1650],
                "page_render_sha256": "C34D31F7CDC4AA1FF2EABB308596C3F52ED4E0ACDF1D64CC260D057D03270074",
                "crop_file": "nasa-mhd-cycle-fig2-source-crop.png",
                "crop_box_page_pixels_xyxy": [210, 130, 1000, 780],
                "crop_dimensions_pixels": [790, 650],
                "crop_sha256": "63CA9B3A2EE33FBA23BBA3E3A3C882C1CCD18AF278E957A942FDFFB63FFE745A",
                "note": "Faithful raster crop of the cited NASA report page for source inspection; distinct from the styled graph redraw.",
            },
            "axis_calibration": {
                "coordinate_origin": "top-left of the cropped PNG; pixel coordinates are zero-based",
                "approximate_plot_left_x_px": 52,
                "approximate_plot_right_x_px": 514,
                "mach_at_left_edge": 1.0,
                "mach_at_right_edge": 12.0,
                "approximate_plot_top_y_px": 55,
                "approximate_plot_bottom_y_px": 453,
                "specific_thrust_at_top_kN_s_per_kg": 1.2,
                "specific_thrust_at_bottom_kN_s_per_kg": 0.0,
                "calibration_coordinate_tolerance_px": "about ±1 px from rasterized axis-line centers",
                "mapping": "Mach = 1 + (x_px - 52) * 11 / 462; specific thrust = (453 - y_px) * 1.2 / 398.",
                "reading_protocol": "Visually follow each printed curve centerline at the listed Mach values, convert with the linear endpoint calibration, round readings to 0.01 kN·s/kg, and connect with straight segments only for display. Near overlapping strokes, centerline identity is ambiguous; listed reading tolerances are visual estimates, not statistical uncertainty.",
            },
        },
        "direct_raster_adjudication": {
            "status": "Provisional native-raster observations shown alongside the earlier approximate redraw. They audit our graph readings, not NASA's numeric source data or model.",
            "data_file": "figure2-direct-raster-trace.csv",
            "data_file_sha256": hashlib.sha256((OUT / "figure2-direct-raster-trace.csv").read_bytes()).hexdigest(),
            "review_record": "30 integer source-crop pixels: 23 clear and 7 explicitly ambiguous. The reproducible adjudication script and point-to-ink audit are in the project analysis archive.",
            "source_pdf_sha256": "12867527B3417F5F99A39C2BAE26B663A2D297157C565CCF3F6D4A1402F41BDF",
            "render": "Direct 300-dpi raster of PDF page 27; source crop and tick-fit calibration are retained with the project audit.",
            "clear_segment_points_displayed": len(REVIEWED_SEGMENTS),
            "ambiguous_points_retained": 7,
            "selection": "23 clear center-pixel readings are displayed as individual markers, not connected. Seven ambiguous observations remain in the downloadable project audit CSV, including crossings, alternatives near Mach 8, and the terminal point.",
            "marker_uncertainty": "Vertical bars show the thresholded dark-ink band at the same raster column; horizontal bars show ±0.5 source pixel. Axis tick-center/calibration uncertainty, antialiasing, original line width, and selection uncertainty are not fully propagated. Bars are not confidence intervals or total uncertainty.",
            "comparison_metric": "The former 23-point discrepancy metric against the old Site line is withdrawn because target-stroke identity was not established. Updated project-audit comparisons quantify manual reading differences only: C n=22 MAE 0.13375, RMSE 0.13848; D n=21 MAE 0.05376, RMSE 0.06151 kN·s/kg. These are not NASA source errors or physical model discrepancies.",
            "provenance_caveat": "Every selected direct coordinate is checked against a same-column dark-ink band in the frozen crop. This supports raster correspondence, but the clear labels, line-center picks, calibration, and comparisons remain manual analyst judgments; crossing/terminal readings are ambiguous.",
            "interpretation": "The green markers are an audit of project digitizations against the published raster. They are not tabulated NASA output, experimental measurements, a validation of the NASA model, or evidence the NASA report is wrong."
        },
        "source_cross_reference_note": "Figure 2 discussion in the final report (printed p. 20 / PDF p. 26) points to Eqs. 25 and 28, which identify diffuser relations in the final equation sequence. By content, the generator and accelerator pressure-ratio relations are Eqs. 36 and 53. The 1999 NASA preprint uses earlier numbering: Eqs. 26 and 30 give those component relations, and its Figure 2 discussion points to those numbers. This suggests stale cross-references in the final report, though the cause is not documented. The reconstruction uses Eqs. 36 and 53 by component definition. The comparison does not validate the source plot or this reconstruction.",
        "plotted_quantity": "Specific thrust F / m_dot_air in kN·s/kg, not total engine thrust and not aircraft acceleration.",
        "swept_parameter": "Generator enthalpy-extraction ratio η_N(g) = 0, 0.25, 0.5. This is a model input, not a measured device efficiency or electrical output fraction.",
        "reported_assumptions": {
            "freestream_static_temperature_K": 220,
            "gas_model": "Thermally and calorically perfect gas; detailed thermochemical effects are not resolved.",
            "fuel_heating_value_MJ_per_kg": 45,
            "combustion_efficiency": 1.0,
            "diffuser_total_pressure_ratio": 0.70,
            "burner_total_pressure_ratio": 0.95,
            "nozzle_total_pressure_ratio": 0.98,
            "conventional_component_pressure_ratios": "Held constant with flight Mach for simplicity; the authors note real ratios vary with Mach and are not independent.",
            "generator_isentropic_efficiency": 0.90,
            "accelerator_isentropic_efficiency": 0.90,
            "mhd_efficiencies": "Optimistic assumed values in the paper.",
            "preionizer_power_fraction": 0.05,
            "preionizer_power_fraction_note": "Chosen for demonstration because the required power was difficult to estimate.",
            "burner_exit_stagnation_temperature_limit_K": 3000,
        },
        "interpretation": "In this specific simplified model, MHD extraction reduces peak specific thrust but shifts the modeled thrust envelope to higher flight Mach numbers. This is a cycle-study result, not a flight test or an established net vehicle benefit. A separate Riggins 2004 peer-reviewed journal analysis reports less favorable thermally balanced simulations for its own configuration (abstract-level record; no numerical curve data available). Do not combine it with this 2000 curve.",
        "reproduce": "Run reproduce_nasa_mhd_cycle_fig2.py with Python 3. The script recreates the CSV and both SVG redraws from the rounded visually digitized points embedded in the script.",
    }
    (OUT / "nasa-mhd-cycle-fig2-methods.json").write_text(json.dumps(methods, indent=2) + "\n", encoding="utf-8")


def svg(mobile: bool) -> str:
    width, height = (840, 690) if mobile else (1200, 720)
    left, right, top, bottom = (84, 28, 110, 480) if mobile else (105, 355, 52, 610)
    x0, x1, y0, y1 = 1.0, 12.0, 0.0, 1.2
    plot_w, plot_h = width - left - right, bottom - top

    def X(x: float) -> float:
        return left + (x - x0) / (x1 - x0) * plot_w

    def Y(y: float) -> float:
        return bottom - (y - y0) / (y1 - y0) * plot_h

    e: list[str] = [
        f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {width} {height}" role="img" aria-labelledby="title desc">',
        '<title id="title">NASA MHD-bypass specific-thrust model with provisional raster audit points</title>',
        '<desc id="desc">Three approximate curves read from NASA Figure 2 are shown as lines. Green individual markers show 23 clear direct readings of dark pixels assigned to the one-quarter extraction stroke. They are not connected. Seven ambiguous readings are omitted from the plot and retained in the audit data. Vertical bars show same-column dark-ink bands and horizontal bars show half a source pixel; neither includes full calibration or selection uncertainty. These markers audit project readings, not NASA source data or model accuracy.</desc>',
        '<rect width="100%" height="100%" rx="18" fill="#08172d"/>',
    ]
    # Plot grid and axes.
    for tick in range(0, 7):
        v = tick * .2
        y = Y(v)
        e.append(f'<line x1="{left}" y1="{y:.1f}" x2="{width-right}" y2="{y:.1f}" stroke="#344968" stroke-width="1"/>')
        e.append(f'<text x="{left-14}" y="{y+6:.1f}" text-anchor="end" fill="#e8efff" font-size="18">{v:.1f}</text>')
    x_ticks = [1, 3, 5, 7, 9, 11, 12] if mobile else range(1, 13)
    for m in x_ticks:
        x = X(m)
        e.append(f'<line x1="{x:.1f}" y1="{top}" x2="{x:.1f}" y2="{bottom}" stroke="#263b5b" stroke-width="1"/>')
        e.append(f'<text x="{x:.1f}" y="{bottom+27}" text-anchor="middle" fill="#e8efff" font-size="16">{m}</text>')
    e.extend([
        f'<line x1="{left}" y1="{top}" x2="{left}" y2="{bottom}" stroke="#d6e5ff" stroke-width="2"/>',
        f'<line x1="{left}" y1="{bottom}" x2="{width-right}" y2="{bottom}" stroke="#d6e5ff" stroke-width="2"/>',
        f'<text x="{left+plot_w/2:.1f}" y="{bottom+61}" text-anchor="middle" fill="#f5f7ff" font-size="20" font-weight="700">Flight Mach number, Mₐ</text>',
        f'<text transform="translate(28 {top+plot_h/2:.1f}) rotate(-90)" text-anchor="middle" fill="#f5f7ff" font-size="19" font-weight="700">Specific thrust, F / ṁₐ (kN·s/kg)</text>',
    ])
    for extraction, points in POINTS.items():
        color, dash, _ = SERIES[extraction]
        coords = " ".join(f"{X(x):.1f},{Y(y):.1f}" for x, y in points)
        dash_attr = "" if dash == "solid" else f' stroke-dasharray="{dash}"'
        e.append(f'<polyline points="{coords}" fill="none" stroke="{color}" stroke-width="4" stroke-linejoin="round" stroke-linecap="round"{dash_attr}/>')
        # Show selected source-read points so the approximate digitization is visible.
        for x, y in points[::4]:
            e.append(f'<circle cx="{X(x):.1f}" cy="{Y(y):.1f}" r="3.2" fill="#08172d" stroke="{color}" stroke-width="2"/>')
    # Keep the approximate Site curves visible; direct raster readings remain
    # individual audit markers and are deliberately not connected into a line.
    review_color = "#b8f28b"
    for mach, thrust, _pixel_x, _pixel_y, y_px_allowance in REVIEWED_SEGMENTS:
        x, y = X(mach), Y(thrust)
        xerr = plot_w / (x1 - x0) * (0.5 / 82.89160839160839)
        yerr = plot_h / (y1 - y0) * y_px_allowance
        e.append(f'<line x1="{x-xerr:.2f}" y1="{y:.2f}" x2="{x+xerr:.2f}" y2="{y:.2f}" stroke="{review_color}" stroke-width="2"/>')
        e.append(f'<line x1="{x:.2f}" y1="{y-yerr:.2f}" x2="{x:.2f}" y2="{y+yerr:.2f}" stroke="{review_color}" stroke-width="2"/>')
        e.append(f'<circle cx="{x:.2f}" cy="{y:.2f}" r="5.5" fill="{review_color}" stroke="#08172d" stroke-width="1.5"/>')
    legend_y = 75 if mobile else 155
    if mobile:
        e.append(f'<text x="{width/2:.1f}" y="{height-144}" text-anchor="middle" fill="#cbd8ed" font-size="16">Generator enthalpy-extraction ratio, ηₙ(g)</text>')
        cols = [135, 420, 690]
        for (ratio, (color, dash, label)), lx in zip(SERIES.items(), cols):
            dash_attr = "" if dash == "solid" else f' stroke-dasharray="{dash}"'
            e.append(f'<line x1="{lx-48}" y1="{height-122}" x2="{lx-12}" y2="{height-122}" stroke="{color}" stroke-width="4"{dash_attr}/>')
            e.append(f'<text x="{lx-4}" y="{height-116}" fill="#f5f7ff" font-size="15">{ratio:g}</text>')
        e.append(f'<circle cx="{width/2-164}" cy="{height-81}" r="5" fill="{review_color}"/>')
        e.append(f'<text x="{width/2-150}" y="{height-76}" fill="#f5f7ff" font-size="15">Direct raster readings · provisional</text>')
    else:
        e.append(f'<text x="{width-right+30}" y="{legend_y-45}" fill="#d4def4" font-size="16" font-weight="700">GENERATOR EXTRACTION</text>')
        for i, (ratio, (color, dash, label)) in enumerate(SERIES.items()):
            y = legend_y + i * 40
            dash_attr = "" if dash == "solid" else f' stroke-dasharray="{dash}"'
            e.append(f'<line x1="{width-right+30}" y1="{y}" x2="{width-right+76}" y2="{y}" stroke="{color}" stroke-width="4"{dash_attr}/>')
            e.append(f'<text x="{width-right+90}" y="{y+6}" fill="#f5f7ff" font-size="18">{label}</text>')
        review_y = legend_y + len(SERIES) * 40
        e.append(f'<line x1="{width-right+30}" y1="{review_y}" x2="{width-right+76}" y2="{review_y}" stroke="{review_color}" stroke-width="2"/>')
        e.append(f'<circle cx="{width-right+53}" cy="{review_y}" r="5" fill="{review_color}" stroke="#08172d" stroke-width="1"/>')
        e.append(f'<text x="{width-right+90}" y="{review_y+6}" fill="#f5f7ff" font-size="16">Direct raster readings (n=23)</text>')
        e.append(f'<text x="{width-right+30}" y="{review_y+30}" fill="{review_color}" font-size="13">Audit of our readings, not NASA data</text>')
        e.append(f'<text x="{width-right+30}" y="{review_y+58}" fill="#cbd8ed" font-size="14">Model assumptions:</text>')
        e.append(f'<text x="{width-right+30}" y="{review_y+82}" fill="#cbd8ed" font-size="14">ηs,g = ηs,a = 0.9</text>')
        e.append(f'<text x="{width-right+30}" y="{review_y+104}" fill="#cbd8ed" font-size="14">T₀,lim = 3000 K</text>')
    note_y = height - 14 if mobile else height - 14
    e.append(f'<text x="{width/2:.1f}" y="{note_y}" text-anchor="middle" fill="#c4d0e6" font-size="13">Lines are approximate graph readings · markers are provisional raster observations, not NASA tabulated data</text>')
    e.append('</svg>')
    return "\n".join(e) + "\n"


def main() -> None:
    for extraction, points in POINTS.items():
        assert len(points) >= 20 and all(1 <= x <= 12 and 0 <= y <= 1.2 for x, y in points)
        assert all(points[i][0] < points[i + 1][0] for i in range(len(points) - 1))
        peak = max(range(len(points)), key=lambda i: points[i][1])
        assert all(points[i][1] >= points[i + 1][1] for i in range(peak, len(points) - 1))
    assert max(y for _, y in POINTS[0.00]) > max(y for _, y in POINTS[0.25]) > max(y for _, y in POINTS[0.50])
    assert 7.7 <= POINTS[0.00][-1][0] <= 8.1
    assert 8.9 <= POINTS[0.25][-1][0] <= 9.3
    assert 10.8 <= POINTS[0.50][-1][0] <= 11.3
    write_csv()
    write_methods()
    audit_methods = {
        "title": "NASA Figure 2 direct-raster reading audit",
        "scope": "Checks DMGB project figure readings against visible pixels in NASA/TP-2000-210387 Figure 2; it does not validate NASA's source numbers, model, or conclusions.",
        "source": "NASA/TP-2000-210387, Figure 2, printed p. 21 / PDF p. 27",
        "source_pdf_sha256": "12867527b3417f5f99a39c2bae26b663a2d297157c565ccf3f6d4a1402f41bdf",
        "data_file": "figure2-direct-raster-trace.csv",
        "data_file_sha256": hashlib.sha256((OUT / "figure2-direct-raster-trace.csv").read_bytes()).hexdigest(),
        "selection": "30 manually selected integer pixels in a frozen 300-dpi source crop; 23 clear center-pixel reads and 7 explicit ambiguities. Values are calculated from recorded pixels after selection using a least-squares fit to major axis ticks.",
        "pixel_check": "Each selected point is checked for grayscale below 100 and membership in a contiguous same-column dark-ink band. This establishes visible-ink correspondence at the chosen threshold, not unique curve identity or statistical validity.",
        "omitted_uncertainty": "Vertical bars show only the same-column thresholded ink band and horizontal bars show ±0.5 source pixel. Axis tick-center/calibration error, antialiasing, original stroke width, and manual selection uncertainty are not fully propagated.",
        "ambiguities": "Four crossing/crowding stations at Mach 4.50–5.25; two alternative adjacent bands at Mach 8; one uncertain near-baseline terminal point at Mach 9.",
        "withdrawn_result": "The earlier 23-point comparison to the approximate Site line (MAE 0.108884, RMSE 0.117331, maximum 0.167021 near Mach 7) is withdrawn as a validated discrepancy because stroke identity was not established.",
        "current_comparison": "Comparisons of traces C and D with this manual raster center trace are descriptive graph-reading differences only. They are not NASA source errors, model-error estimates, or ground truth.",
        "reproduction": "The underlying project script analyses/nasa-2000-independent-blind/raster-adjudication/trace_and_compare.py regenerates the audit data and overlay; see the project analysis report for commands and tests."
    }
    (OUT / "figure2-raster-adjudication-methods.json").write_text(json.dumps(audit_methods, indent=2) + "\n", encoding="utf-8")
    (OUT / "nasa-mhd-cycle-fig2.svg").write_text(svg(False), encoding="utf-8")
    (OUT / "nasa-mhd-cycle-fig2-mobile.svg").write_text(svg(True), encoding="utf-8")
    print(f"Wrote {sum(len(v) for v in POINTS.values())} visually digitized points and two SVG redraws; source data are not tabulated in NASA/TP-2000-210387.")


if __name__ == "__main__":
    main()