"""Sum haversine segment lengths of every <trkpt> in a GPX file.
Usage: python3 -I route_length.py FILE.gpx [radius_m]   (default radius 6371008.8 m)"""
import math, sys
import xml.etree.ElementTree as ET
NS = {"g": "http://www.topografix.com/GPX/1/1"}
R = float(sys.argv[2]) if len(sys.argv) > 2 else 6371008.8
pts = [(float(p.get("lat")), float(p.get("lon"))) for p in ET.parse(sys.argv[1]).getroot().iterfind(".//g:trkpt", NS)]
def hav(a, b):
    la1, lo1, la2, lo2 = map(math.radians, (*a, *b))
    h = math.sin((la2 - la1) / 2) ** 2 + math.cos(la1) * math.cos(la2) * math.sin((lo2 - lo1) / 2) ** 2
    return 2 * R * math.asin(math.sqrt(h))
print(f"{len(pts)} points, {sum(hav(a, b) for a, b in zip(pts, pts[1:])):.1f} m (radius {R} m)")
