I'm brand new to geo processing/mapping, and I'm struggling with spherical trig, trying to understand what's going on without even knowing the terminology of what I'm looking for. After a lot of struggle, I've accumulated a bunch of code from snippets and other people's projects and assembled it into a rhumb line calculator (included below for reference) starting from a point and travelling in a direction, but (I think) ultimately I need to do this for a great circle instead.
The question is, if I'm at (lat, lon), and I travel $d$ distance in meters at $\theta$ angle (from due east), what point am I at now?
All of this is on the earth's surface, and while it isn't critical that it be pinpoint accurate, if I could at least understand the math enough to intuit what that inaccuracy would look like, that would be ideal.
Could you point me to resources or formulas for these kinds of operations?
I've included the below Frankenstein code. There are definitely pieces of it I don't know the reasoning behind, and I really want to understand the math involved! -- which is also why I posted here on the math stack exchange, rather than the GIS one.
import math
EARTH_RADIUS = 6378137.0
FLATTENING = 1/298.257223563
EARTH_MINOR = EARTH_RADIUS * (1-FLATTENING)
def rhumb(lat: float, lon: float, heading: float, dist: float) -> tuple[float, float]:
'''
calculate the point on the rhumb line at a given distance (in meters)
lat, lon - in decimal degrees
heading - decimal degrees (E of N)
dist - meters
Example:
>>> rhumb(33.0, -110.06, 45, 1000)
(33.006375586750394, -110.05243303029485)
'''
alpha1 = math.radians(heading)
sin_alpha1 = math.sin(alpha1)
cos_alpha1 = math.cos(alpha1)
tan_U1 = (1-FLATTENING) * math.tan(math.radians(lat))
cos_U1 = 1 / math.sqrt(1 + tan_U1**2)
sin_U1 = tan_U1 * cos_U1
sigma1 = math.atan2(tan_U1, cos_alpha1)
sin_alpha = cos_U1 * sin_alpha1
cosSq_alpha = 1 - sin_alpha**2
uSq = cosSq_alpha * (EARTH_RADIUS**2 - EARTH_MINOR**2) / (EARTH_MINOR**2)
# I have no idea what A and B represent here, but there are a lot of
# powers of two and uSq -> some kind error correction?
# A is slightly (0.001) larger than 1, B is slightly (0.001) larger than zero
A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
sigma = dist / (EARTH_MINOR * A)
sigmaP = 2 * math.pi
sin_sigma = math.sin(sigma)
# Error correction?
while abs(sigma-sigmaP) > 1e-12:
cos_2sigma_m = math.cos(2 * sigma1 + sigma)
sin_sigma = math.sin(sigma)
cos_sigma = math.cos(sigma)
delta_sigma = B * sin_sigma * (
cos_2sigma_m + B / 4 * (cos_sigma * (-1 + 2 * cos_2sigma_m**2) -
B / 6 * cos_2sigma_m * (-3 + 4 * sin_sigma**2) * (-3 + 4 * cos_2sigma_m**2)))
sigmaP = sigma
sigma = dist / (EARTH_MINOR * A) + delta_sigma
# final lat
tmp = sin_U1 * sin_sigma - cos_U1 * cos_sigma * cos_alpha1
final_lat = math.atan2(
sin_U1 * cos_sigma + cos_U1 * sin_sigma * cos_alpha1,
(1-FLATTENING) * math.sqrt(sin_alpha**2 + tmp**2))
final_lat = math.degrees(final_lat)
# final lon
lambda1 = math.atan2(sin_sigma * sin_alpha1,
cos_U1 * cos_sigma - sin_U1 * sin_sigma * cos_alpha1)
C = FLATTENING / 16 * cosSq_alpha * \
(4 + FLATTENING * (4 - 3 * cosSq_alpha))
L = lambda1 - (1-C) * FLATTENING * sin_alpha * (
sigma + C * sin_sigma * (
cos_2sigma_m + C * cos_sigma * (
-1 + 2 * cos_2sigma_m**2)))
final_lon = lon + math.degrees(L)
return final_lat, final_lon
```
After a ton of searching and trial and error, I found an answer here: http://www.edwilliams.org/avform147.htm#LL
It's a flight navigation aid!
As well, I found some information about the method used in the code I inherited, called the Vincenty method here: https://www.movable-type.co.uk/scripts/latlong-vincenty.html
it's for highly accurate (i.e. < mm) accuracy when calculating geos.