>>105 Maxima
f(x1,y1,x2,y2):=(r:%pi/180,
6378.137*acos(sin(y1*r)*sin(y2*r)+cos(y1*r)*cos(y2*r)*cos(x2*r-x1*r)));

(%i60) float(f(139.766086,35.681377,-58.381595,-34.603722));
(%o60) 18383.55983532794