Well, the angle between 2 points in a cartesian axis (theta) is the angle counterclockwise from the x-axis, but bearing (azimuth) is the angle clockwise from North. So azimuth = 90 - theta
You can do all sorts of mapping with GMT.jl See for example this example and the corresponding julia script.
And as a simpler example, let’s plot an orthodrome between two points
ortho = project(start_pt=[-50.0 10], end_pt=[50 60], step=10, km=true);
coast(region=[-180 180 -90 90], proj=:equidistCylindrical, frame=:g, res=:crude, land=:navy, B=:a)
plot!(ortho, lw=3, lc=:green, name="ortho_line.png", show=true)
(to plot an loxodrome one could interpolate linearly between two points in a Mercator, convert back to gepgraphic (using mapproject) and add that line to the plot below)
