= Symbol("t")
lat, lon = Function("lat")(t), Function("lon")(t)

Physical constants for orbital mechanics

R_earth = 6371.0

def haversine_symbolic(lat1, lon1, lat2, lon2):
dlat = lat2 - lat1
dlon = lon2 - lon1
a = sin(dlat/2)2 + cos(lat1)cos(lat2)sin(dlon/2)2
return 2R_earthasin(sqrt(a))

Ground track velocity

v_ground = diff(haversine_symbolic(lat0, lon0, lat, lon), t)