@lazarus, please see code below that computes the geodesic line over the GRS80 ellipsoid given the (LAT, LONG) of 2 initial points.
This is a quick&dirty attempt to implement the differential equations in Panou et al. (2013) with some simplifications. In particular, the iterative process to fine tune one of the boundary conditions was simplified, hopefully not too much.
In any case, the results seem to match the authors’ numeric tables for different points tested, in terms of geodesic length and Clairaut constants. The results seem to match also those from this web calculator.
Not being an expert in any of these things, much less in Julia, any advice on how to improve the differential equations accuracy and the code are very welcome.
# Reference:
# Panou, G., Delikaraoglou, D. and Korakitis, R. (2013)
# Solving the geodesics on the ellipsoid as a boundary value problem.
# Journal of Geodetic Science, 3(1)
using DifferentialEquations
function geodesic!(du,u,p,λ)
du[1] = u[2]
du[2] = h1(u[1])*u[2]^2 + h2(u[1])
du[3] = u[4]
du[4] = (h3(u[1])*u[2]^2 + h4(u[1]))*u[3] + 2*h1(u[1])*u[2]*u[4]
end
trapez(t,A) = @views sum((A[1:end-1] + A[2:end]).*diff(t)/2) # trapezoidal integration
# GRS80 ellipsoid (Moritz, 1980)
a = 6378137; b = 6356752.3141 # m
e = sqrt(a^2 - b^2)/a; e2 = e^2
N(φ) = a / sqrt(1 - e2*sin(φ)^2) # radius of curvature in prime vertical normal section
E(φ) = a^2*(1 - e2)^2 / (1 - e2*sin(φ)^2)^3 # first fundamental coefficients E, F=0, G
G(φ) = a^2*cos(φ)^2 / (1 - e2*sin(φ)^2)
h1(φ) = -(2*(1 - e2)*tan(φ)+ 3*e2*sin(φ)*cos(φ)) / (1 - e2*sin(φ)^2)
h2(φ) = -(1 - e2*sin(φ)^2)*sin(φ)*cos(φ) / (1 - e2)
h3(φ) = -2*(1 - e2)*(1 - e2*sin(φ)^2 + 2*e2*sin(φ)^2 * cos(φ)^2) / ((1 - e2*sin(φ)^2)^2 * cos(φ)^2) +
- 3*e2*((1 - e2*sin(φ)^2)*(cos(φ)^2 - sin(φ)^2) + 2*e2*sin(φ)^2 * cos(φ)^2) / (1 - e2*sin(φ)^2)^2
h4(φ) = -((1 - e2*sin(φ)^2)*(cos(φ)^2 - sin(φ)^2) - 2*e2*sin(φ)^2 * cos(φ)^2) / (1 - e2)
# INPUT geodetic endpoints (LAT, LONG), convert from degrees to radians
# NOTE: geodetic LAT: –π/2 < φ < π/2; geodetic LONG: –π < λ < π
d2r = π/180
φ0 = 60*d2r ; λ0 = 0*d2r # λ0 < λ1
φ1 = -30*d2r ; λ1 = 120*d2r # λ0 = λ1 (meridian) excluded as trival case
# initial conditions
u0 = zeros(4)
u0[1] = φ0
u0[2] = (φ1 - φ0)/(λ1 - λ0) # need to iterate to find this initial condition so that geodetic ends at φ1(λ1)
u0[3] = 0.
u0[4] = 1.
λspan = (λ0, λ1) # longitude-span for ODE solving (takes the role of time variable)
prob = ODEProblem(geodesic!, u0, λspan) # define problem
sol = solve(prob, alg=RK4())
φ10 = sol[1,end]
n = 8000
λ = LinRange(λ0,λ1,n)
dϕ = 0.05*u0[2] # perturb dφ/dλ by 5% (empirical, can be improved)
edφ = +Inf
while abs(edφ) > 1e-12
u0[2] += dϕ
prob = ODEProblem(geodesic!, u0, λspan) # define problem
sol = solve(prob, alg=RK4(),saveat=λ,reltol=1e-12)
φ11 = sol[1,end]
edφ = φ1 - φ11
dφ1dφ0 = (φ11 - φ10)/dϕ # simplification of Panou et al., 2013, Eq.33)
dϕ = edφ/dφ1dφ0
φ10 = φ11
end
λ = sol.t
φ = sol[1,:]
dφ = sol[2,:]
# geodetic distance : s
A = @. sqrt(E(φ)*(dφ)^2 + G(φ))
s = trapez(λ, A) # geodetic distance
# Clairaut constant : C
C = mean(@. G(φ)/sqrt(E(φ)*(dφ)^2 + G(φ)))
# Plot computed geodesic line on ellipsoid using GMT
using Printf, GMT
str = @sprintf("%i m, over (LAT,LONG) = (%i,%i) -> (%i,%i) deg",s,φ0/d2r,λ0/d2r,φ1/d2r,λ1/d2r);
coast(region=:global, proj=(name=:ortho, center=((λ0+λ1)/d2r/2,(φ0+φ1)/d2r/2)), land=:lightblue, frame=:g,
title="Geodesic L=$str", par=(FONT_TITLE="14p,Helvetica,black",))
GMT.plot!([λ/d2r φ/d2r], lw=1.5, lc=:red, show=true)
# plot computed geodesic line on ellipsoid using Plots.jl pyplot()
using Printf, Plots; pyplot()
n = 100; u = collect(LinRange(0,2*π,n)); v = collect(LinRange(0,π,n));
x = a*cos.(u)*sin.(v)'; y = a*sin.(u)*sin.(v)'; z = b*ones(n)*cos.(v)';
Plots.surface(x,y,z, colorbar=false)
x0 = @. N(φ)*cos(φ)*cos(λ);
y0 = @. N(φ)*cos(φ)*sin(λ);
z0 = @. N(φ)*(1 - e2)*sin(φ);
str = @sprintf("%i m, over (LAT,LONG) = (%i,%i) -> (%i,%i) deg",s,φ0/d2r,λ0/d2r,φ1/d2r,λ1/d2r);
plot!(x0,y0,z0, lc=:red, label="geodesic L=$str")
Off course. I tried to see if there was an option to generate loxodromes that I forgot the obvious.
)
