Last active
February 16, 2016 20:32
-
-
Save alexjohnj/eacbece12a6447b4defa to your computer and use it in GitHub Desktop.
Julia script to calculate length and azimuthal angles of a great circle joining two points on the Earth (ideally for seismology).
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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