"""Reproduce an ideal axial dipole teaching example; not measured data.

Run: python reproduce.py (requires numpy and matplotlib).
Outputs are written beside this script. No sensor or object model is fitted.
"""
from pathlib import Path
import csv
import json
import math
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter

OUT = Path(__file__).resolve().parent
MU0 = 4 * math.pi * 1e-7  # H/m, a teaching approximation, not an exact SI constant.
MOMENT = 0.01  # A m^2, assigned rather than measured.

def axial_nt(moment, distance):
    """Positive axial component in nT for an aligned ideal point dipole."""
    return MU0 * moment / (2 * math.pi * np.asarray(distance)**3) * 1e9

def main():
    r = np.linspace(0.1, 1.0, 181)
    b = axial_nt(MOMENT, r)
    pair = [(0.01, 0.2), (0.08, 0.4)]
    assert math.isclose(axial_nt(MOMENT, 0.1), 2000, rel_tol=1e-12)
    assert math.isclose(axial_nt(MOMENT, 1.0), 2, rel_tol=1e-12)
    for distance in [0.1, 0.2, 0.3, 0.5]:
        assert math.isclose(axial_nt(MOMENT, distance) / axial_nt(MOMENT, 2*distance), 8, rel_tol=1e-12)
    assert all(math.isclose(axial_nt(m, d), 250, rel_tol=1e-12) for m, d in pair)
    with (OUT / 'ideal-dipole.csv').open('w', newline='', encoding='utf-8') as f:
        writer = csv.writer(f)
        writer.writerow(['distance_m', 'assigned_moment_A_m2', 'ideal_axial_B_nT'])
        writer.writerows((f'{d:.3f}', MOMENT, f'{field:.12g}') for d, field in zip(r, b))
    with (OUT / 'equal-reading-pair.csv').open('w', newline='', encoding='utf-8') as f:
        writer = csv.writer(f)
        writer.writerow(['scenario', 'distance_m', 'assigned_moment_A_m2', 'ideal_axial_B_nT'])
        writer.writerows((name, d, m, f'{axial_nt(m,d):.12g}') for name, (m,d) in zip(['A','B'], pair))

    plt.rcParams.update({'font.family': 'DejaVu Sans', 'font.size': 12, 'svg.fonttype': 'none', 'pdf.fonttype': 42})
    for mobile in [False, True]:
        fig, axes = plt.subplots(2 if mobile else 1, 1 if mobile else 2,
            figsize=(5.2, 9.8) if mobile else (12.2, 6.5))
        fig.subplots_adjust(left=0.18 if mobile else 0.085, right=0.96, bottom=0.22 if mobile else 0.34,
                            top=0.86 if mobile else 0.79, hspace=0.58, wspace=0.32)
        fig.suptitle('An ideal magnetic dipole\nDistance matters; one reading is ambiguous', fontsize=15, y=0.97)
        ax = axes[0]
        ax.plot(r, b, color='#14638c', lw=2.6)
        ax.set_yscale('log')
        ax.set_yticks([2, 10, 100, 1000, 2000])
        ax.yaxis.set_major_formatter(ScalarFormatter())
        ax.set_xlim(0.1, 1)
        ax.set_ylim(1.6, 2600)
        ax.set_xlabel('Axial distance r (m)')
        ax.set_ylabel('Ideal axial B (nT) — logarithmic scale')
        ax.set_title('Assigned m = 0.01 A m²', fontsize=13, pad=12)
        ax.grid(True, which='major', alpha=0.25)
        ax.text(0.95, 0.92, 'Double distance →\none-eighth of B', transform=ax.transAxes, ha='right', va='top', fontsize=12)
        ax = axes[1]
        ax.bar([0, 1], [250, 250], width=0.55, color=['#14638c', '#a35016'])
        ax.set_xticks([0, 1], ['A\n0.01 A m²\nat 0.20 m', 'B\n0.08 A m²\nat 0.40 m'])
        ax.set_ylim(0, 330)
        ax.set_ylabel('Ideal axial B (nT) — linear scale')
        ax.set_title('Different scenarios, same B', fontsize=13, pad=12)
        ax.grid(True, axis='y', alpha=0.25)
        ax.set_axisbelow(True)
        for i in [0, 1]:
            ax.text(i, 260, '250 nT', ha='center', fontsize=12)
        footer = ('Teaching model, not a measurement or object identification.\n'
                  'B = μ₀m/(2πr³); μ₀ ≈ 4π × 10⁻⁷ H/m.\n'
                  'Ideal point dipole / far-field approximation, on its axis.\n'
                  'No finite geometry, background, noise, off-axis\n'
                  'or induced-magnetization effects are included.')
        fig.text(0.07 if mobile else 0.085, 0.018, footer, fontsize=10.5 if mobile else 11, va='bottom', linespacing=1.5)
        stem = 'ideal-dipole-mobile' if mobile else 'ideal-dipole'
        fig.savefig(OUT / (stem + '.svg'))
        fig.savefig(OUT / (stem + '.pdf'))
        plt.close(fig)
    methods = {
        'title': 'Ideal axial dipole reading: a civilian teaching model',
        'data_kind': 'constructed values, not measurements',
        'formula': 'B_axial = mu0*m/(2*pi*r^3)',
        'units': {'B': 'tesla; plotted and tabulated in nanotesla', 'm': 'A m^2', 'r': 'm', 'mu0': 'H/m'},
        'mu0_used': MU0,
        'mu0_note': '4*pi*1e-7 H/m is an adequate rounded teaching approximation. In the revised SI, mu0 is experimentally determined rather than fixed exactly to this value. These ideal constructed outputs do not justify comparable measurement precision.',
        'assigned_moment_A_m2': MOMENT,
        'range_m': [0.1, 1.0], 'sample_count': len(r), 'step_m': 0.005,
        'orientation': 'Positive axial component for an aligned dipole; no off-axis direction is modeled.',
        'applicability': 'An ideal point dipole is used. For a finite source, the dipole approximation requires distance much larger than its characteristic dimensions. No source size is assigned, so applicability to a real source over 0.1–1 m is not established.',
        'equal_pair': [{'moment_A_m2': m, 'distance_m': d, 'B_nT': 250} for m,d in pair],
        'interpretation': 'Within this scalar ideal model, differing moment and distance can yield the same reading; a single scalar value does not uniquely identify a source.',
        'excluded': ['finite source geometry', 'background fields', 'noise', 'off-axis effects', 'induced magnetization'],
        'not_claimed': ['actual measurement', 'object identification', 'sensor performance', 'orbital design'],
        'figure_axes': {'left_or_top': 'linear distance; logarithmic B', 'right_or_bottom': 'scenario categories; linear B'},
        'validation': {'endpoint_nT': [2000, 2], 'doubling_distance_B_ratio': 8, 'both_pair_values_nT': 250, 'relative_tolerance': 1e-12, 'passed': True},
        'reproduction': 'python reproduce.py; dependencies: numpy, matplotlib'
    }
    (OUT / 'methods.json').write_text(json.dumps(methods, indent=2) + '\n', encoding='utf-8')
    print('PASS: inverse-cube scaling, endpoints, equal-reading pair; 181 curve rows and 2 pair rows.')

if __name__ == '__main__':
    main()

