Skip to content

Instantly share code, notes, and snippets.

@alexjohnj alexjohnj/finddelta.jl
Last active Feb 16, 2016

Embed
What would you like to do?
Julia script to calculate length and azimuthal angles of a great circle joining two points on the Earth (ideally for seismology).
#!/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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.