Like shown in the examples linked in my post above, you now do that in GMT with (for example). It uses PROJ for the calculations.
ortho = orthodrome([0 0; 70 60], step=1000, unit=:k);
hmm, better show also the result
using GMT
# Compute the great circle (orthodrome) line
ortho = orthodrome([0 0; 70 60], step=1000, unit=:k);
# Plot the orthodrome on an orthographic projection
coast(region=:global, proj=(name=:ortho, center=(0,45)), land=:peru, frame=:g)
plot!(ortho, lw=0.5, marker=:circ, ms=0.1, fill=:black, show=true)
