print "Haversine distance: "; using( "####.###########", havDist( 36.12, -86.67, 33.94, -118.4)); " km."
end
function havDist( th1, ph1, th2, ph2)
  degtorad   = acs(-1)/180
  diameter   = 2 * 6372.8
    LgD      = degtorad  * (ph1 - ph2)
    th1      = degtorad  * th1
    th2      = degtorad  * th2
    dz       = sin( th1) - sin( th2)
    dx       = cos( LgD) * cos( th1) - cos( th2)
    dy       = sin( LgD) * cos( th1)
    havDist  = asn( ( dx^2 +dy^2 +dz^2)^0.5 /2) *diameter
end function