Julia script to calculate length and azimuthal angles of a great circle joining two points on the Earth (ideally for seismology).

finddelta.jl

#!/usr/bin/env julia
#
# A script used to find the length of the great circle joining two points on the
# Earth as well as the azimuth and back-azimuth of the great circle.
#
# USAGE: finddelta LAT1 LONG1 LAT2 LONG2
# OUTPUT: Absolute distance, angular distance, eq azimuth and back azimuth
# EXAMPLE (Great circle between the poles): finddelta -90 0 90 0

const EARTH_RADIUS = 6371_000

immutable Coordinates{T<:AbstractFloat}
    lat::T
    lon::T
end

function parseargs{T<:AbstractString}(S::Array{T})
    if length(S) != 4
        @printf("USAGE: finddelta LAT1 LONG1 LAT2 LONG2\n")
        exit(1)
    end
    cs = map(S) do s
        c = tryparse(BigFloat, s)
        if isnull(c)
            error("Unable to parse coordinates")
        else
            return get(c)
        end
    end
    (Coordinates(cs[1], cs[2]), Coordinates(cs[3], cs[4]))
end

""" Calculate the central angle between two points on a sphere in radians. """
function calccangle(c1::Coordinates, c2::Coordinates)
    acos(sind(c1.lat) * sind(c2.lat) + cosd(c1.lat) * cosd(c2.lat) * cosd(abs(c2.lon - c1.lon)))
end

""" Calculate the length of the great circle joining two points on a sphere. """
function calcdist(c1::Coordinates, c2::Coordinates)
    EARTH_RADIUS * calccangle(c1, c2)
end

""" Calculate the angle the great circle from `c1` to `c2` makes with north
measured clockwise at `c1`. """
function calcazimuth(c1::Coordinates, c2::Coordinates)
    numer = cosd(c2.lat) * sind(c2.lon - c1.lon)
    denom = sind(c2.lat) * cosd(c1.lat) - cosd(c2.lat) * sind(c1.lat) * cosd(c2.lon - c1.lon)

    ξ = rad2deg((atan2(numer, denom) + 2π) % 2π)
end

function main()
    eqloc, statloc = parseargs(ARGS)
    Δ = rad2deg(calccangle(eqloc, statloc))
    dist = calcdist(eqloc, statloc)
    ξ = calcazimuth(eqloc, statloc)
    ξb = calcazimuth(statloc, eqloc)

    @printf("Source-Receiver Distance: %.2f km\nSource-Receiver Angle: %.2f°\n", dist/1000, Δ)
    @printf("Azimuth: %.2f°\nBack-Azimuth: %.2f°\n", ξ, ξb)
end

if !isinteractive()
    main()
end