import json
import numpy as np

# Harmonic oscillator H=(q^2+p^2)/2, dq=p, dp=-q.
def leapfrog(q, p, h):
    p = p - 0.5 * h * q
    q = q + h * p
    p = p - 0.5 * h * q
    return q, p

def map_jacobian(h):
    # Exact linear map from differentiating kick-drift-kick.
    return np.array([
        [1.0 - h*h/2.0, h],
        [-h + h**3/4.0, 1.0 - h*h/2.0],
    ])

def run(h, total_time=200.0):
    q, p = 1.0, 0.0
    values = []
    for _ in range(round(total_time / h)):
        values.append(0.5 * (q*q + p*p))
        q, p = leapfrog(q, p, h)
    values = np.asarray(values)
    return float(values.max() - values.min()), float(abs(values[-1] - values[0]))

def main():
    hs = np.array([0.4, 0.2, 0.1, 0.05, 0.025])
    ranges = np.array([run(float(h))[0] for h in hs])
    slope = float(np.polyfit(np.log(hs), np.log(ranges), 1)[0])
    determinants = [float(np.linalg.det(map_jacobian(float(h)))) for h in hs]
    output = {
        'step_sizes': hs.tolist(),
        'energy_ranges': ranges.tolist(),
        'loglog_energy_range_slope': slope,
        'jacobian_determinants': determinants,
        'symplectic_volume_preservation': bool(max(abs(d - 1.0) for d in determinants) < 1e-14),
        'second_order_energy_scaling_observed': bool(1.7 < slope < 2.3),
    }
    print(json.dumps(output, indent=2))
    with open('verification.json', 'w') as f:
        json.dump(output, f, indent=2)

if __name__ == '__main__':
    main()