#!/usr/bin/env python3
"""Recompute a three-row published-table comparison, not an NR simulation.

Source: Damour, Guercilena, Hinder, Hopper, Nagar, Rezzolla (2014),
Phys. Rev. D 89, 081503, doi:10.1103/PhysRevD.89.081503.
Version: https://arxiv.org/pdf/1402.7307v1, Tables I/II, printed pp. 2/3.
Angles and quoted uncertainties are in degrees. The source supplies error
estimates, not a joint Gaussian likelihood or raw convergence sequences.
Python 3 standard library only. Run: python3 published-scattering-comparison.py
"""
import json
import math

# b_NR/M, Ein/M, Jin/M^2, DeltaE/M, DeltaJ/M^2,
# chi_NR, quoted NR error, chi_EOBNR5PNlog, quoted EOB error,
# chi_PN3PN, quoted PN error. These are factual source-table transcriptions.
ROWS = [
    [10.0,1.0225791,1.145523,.010734,.1164,222.9,1.7,227.0,5.0,126.0,1.0],
    [12.0,1.0225907,1.374658,.001980,.04022,120.7,1.5,120.77,.06,104.34,.04],
    [16.0,1.0225932,1.832883,.0003338,.0162,70.7,1.9,70.666,.002,69.122,.002],
]


def main():
    output=[]
    for b,e,j,de,dj,nr,nr_error,eob,eob_error,pn,pn_error in ROWS:
        assert e>1 and j>0 and de>0 and dj>0
        output.append({
            'b_NR_over_M':b,
            'E_bar_over_M':e-de/2,
            'J_bar_over_M_squared':j-dj/2,
            'NR_angle_degrees':nr,'quoted_NR_error_degrees':nr_error,
            'EOB_angle_degrees':eob,'quoted_EOB_error_degrees':eob_error,
            'PN_angle_degrees':pn,'quoted_PN_error_degrees':pn_error,
            'EOB_minus_NR_degrees':eob-nr,
            'PN_minus_NR_degrees':pn-nr,
            'EOB_fractional_difference_percent':100*(eob-nr)/nr,
            'PN_fractional_difference_percent':100*(pn-nr)/nr,
            'EOB_NR_quadrature_error_degrees':math.hypot(eob_error,nr_error),
            'quadrature_assumption':'uncorrelated errors, illustrative only; not a confidence interval',
        })
    print(json.dumps({'source':'arXiv:1402.7307v1, Tables I/II',
        'scope':'Published 2014 model predictions; no refit, new Hamiltonian integration, or NR rerun',
        'rows':output},indent=2))


if __name__=='__main__':
    main()