sympy hangs when solving for component of Haversine formula

Question

I am using the Haversine formula to calculate a distance in miles between two (lat, lng) coordinate pairs. (Note: I am aware and okay with limitations in the formula related to the non-spheroidal (ellipsoidal) shape of the Earth.)

I would like to use this formula to solve for either a single latitude or single longitude that is due north, east, west, or south of a given coordinate. This is maybe best illustrated through a diagram; I have the central red point as given and am trying to solve for the 4 outer red points below:

From the central coordinate of (38.0, -77.0), I want to solve (individually) for the 4 missing points at each side of the circle pictured, assuming a distance of 5 miles. So in each equation, I am given a distance and 3 coordinates, and want to solve for the 4th coordinate.

What I have tried is to use sympy, but the calculation seems to time out, unless I have a symbol wrong somewhere. To use the top point (lat2, -77.0) as an example:

import sympy as s

lat1 = s.rad(38.0)
lat2 = s.Symbol('lat2')
lon1 = s.rad(-77.0)
lon2 = s.rad(-77.0)
d = 5.0  # Given distance
R = 3950.  # Radius of earth in miles
dlon = lon2 - lon1
dlat = lat2 - lat1
a = (s.sin(dlat/2))**2 + s.cos(lat1) * s.cos(lat2) * (s.sin(dlon/2))**2
c = 2 * s.atan2( s.sqrt(a), s.sqrt(1-a) )
s.solve(3950 * c - d, lat2)  # HANGS

Here I'm trying to solve for lat2, but the .solve() call hangs up indefinitely.

Sorry, I knew that @MichaelButscher ;) Copied the wrong formula here. I've corrected that but am seeing the same problem. — Brad Solomon, Jul 18 '20 at 23:27
Does this answer your question? https://stackoverflow.com/questions/15258078/latitude-longitude-and-meters — Joni, Jul 18 '20 at 23:37
It very well may @Joni, but let me give that a more thorough look before we can close as duplicate — Brad Solomon, Jul 18 '20 at 23:39

score 0 · Accepted Answer · answered Jul 19 '20 at 00:19

Since you have floats in the equation and no symbols you aren't going to get an exact analytic solution so you can just use nsolve. Note that I see different solutions for different initial estimates and they both appear to be solutions. I'm not sure how to interpret that in the context of your problem (maybe it's two of the 4 points?).

In [25]: nsolve(eq, lat2, 0)                                                                                           
Out[25]: 0.661959292973035

In [26]: nsolve(eq, lat2, 1)                                                                                           
Out[26]: 0.664490938542655

In [27]: v1 = nsolve(eq, lat2, 0)                                                                                      

In [28]: v2 = nsolve(eq, lat2, 2)                                                                                      

In [29]: eq.subs(lat2, v1).n()                                                                                         
Out[29]: -7.37924650738106e-14

In [30]: eq.subs(lat2, v2).n()                                                                                         
Out[30]: 1.27083170255818e-13

If you use atan instead of arctan2:

c = 2 * s.atan( s.sqrt(a) / s.sqrt(1-a) )

then you can get those both directly from solve:

In [43]: s.solve(3950 * c - d, lat2, check=False, simplify=False)                                                      
Out[43]: [0.661959292973035, 0.664490938542655]

I'm not sure how it goes wrong with atan2...

sympy hangs when solving for component of Haversine formula

1 Answers1