Vincenty Inverse Part 2

After further investigation, I found that, given locations near antipodal points, the Vincenty geodetic inverse calculations fail miserably. So, I am looking into alternative calculations that have been modified to account for antipodal situations. Vincenty has a paper where certain formulae are rearranged to account for these types of situations. But first, for my own clarification, I will depict the standard calculations, hopefully in a manner clear to me! We begin with converting all values to radians before we start (i.e., \(degrees / (180/\pi)\). Then we start calculations:

\[\begin{align} & U_{1} = \tan^{-1}{(1-f) \tan{\phi_{1}}} \\ & U_{2} = \tan^{-1}{(1-f) \tan{\phi_{2}}} \\ & L = L_{2} - L_{1} \\ & b = (1-f) a \\ \end{align}\]

Then we set an initial value for lambda, \(\lambda = L\), before we start iterating until \(\lambda\) converges:

\[\begin{align} & \sin{\sigma} = \sqrt{(\cos{U_{2} \sin{\lambda})^{2} + (\cos{U_{1}} \sin{U_{2}} - \sin{U_{1}} \cos{U_{2}} \cos{\lambda})^{2}}} \\ & \cos{\sigma} = \sin{U_{1}} \sin{U_{2}} + \cos{U_{1}} \cos{U_{2}} \cos{\lambda} \\ & \sigma = \arctan{\frac{\sin{\sigma}}{\cos{\sigma}}} \\ &\sin{\alpha} = \frac{\cos{U_{1}} \cos{U_{2}} \sin{\lambda}}{\sin{\sigma}} \\ & \cos^{2}{\alpha} = 1 - \sin^{2}{\alpha} \\ & \cos^{2}{(2\sigma_{m})} = \cos{\sigma} - \frac{2 \sin{U_{1}} \sin{U_{2}}}{\cos^{2}{\alpha}} \\ & C = \frac{f}{16} \cos^{2}{\alpha} [4 + f (4 - 3 \cos^{2}{\alpha})] \\ & \lambda = L + (1 - C) f \sin{\alpha} \lbrace \sigma + C \sin{\sigma} [\cos{(2\sigma_{m})} + C \cos{\sigma} (-1+2 \cos^{2}{(2 \sigma_{m}))}]\rbrace \\ \end{align}\] Once \(\lambda\) has converged to the desired accuracy (<1e-12), we continue with the below calculations:

\[\begin{align} & \upsilon^{2} = (\cos^{2}{\alpha}) \biggl(\frac{a^{2} - b^{2}}{b^{2}}\biggr) \\ & A = 1 + \frac{\upsilon^{2}}{16384} \lbrace 4096 + \upsilon^{2} [-768 + \upsilon^{2} (320 - 175 \upsilon^{2})]\rbrace \\ & B = \frac{\upsilon^{2}}{1024} \lbrace 256 + \upsilon^{2} [-128 + \upsilon^{2} (74-47 \upsilon^{2})]\rbrace \\ & \Delta\sigma = B \sin{\sigma} \biggl\lbrace \cos{(2\sigma_{m})} + \frac{1}{4} B \biggl[ \cos{\sigma} (-1+2 \cos^{2}{(2\sigma_{m})}) - \frac{1}{6} B \cos{(2\sigma_{m})} (-3+4 \sin^{2}{\sigma}) (-3+4 \cos^{2}{(2\sigma_{m})}) \biggr] \biggr\rbrace \\ \end{align}\] Then we have the distance calculated:

\[s = b A ( \sigma - \Delta \sigma )\]

Finally we determine the starting and ending bearing:

\[\begin{align} & \alpha_{1} = \arctan{\biggl(\frac{\cos{U_{2}} \sin{\lambda}}{\cos{U_{1}} \sin{U_{2}} - \sin{U_{1}} \cos{U_{2}} \cos{\lambda}}\biggr)} \\ & \alpha_{2} = \arctan{\biggl(\frac{\cos{U_{1}} \sin{\lambda}}{-\sin{U_{1}} \cos{U_{2}} + \cos{U_{1}} \sin{U_{2}} \cos{\lambda}}\biggr)} \\ \end{align}\]

It seems that when the \(\lambda\) is greater than \(\pi\), different methods must be used. The solution I will attempt, also from Vincenty, is used if L\(\mathcal{'} = \pi - L\) if L > 0, or L\(\mathcal{'} = - \pi - L\) if L < 0. To start we set \(\lambda\mathcal{'} = 0, \cos^{2}{\alpha} = 1/2, \cos{2\sigma} = 0\), and \(\sigma = \pi - |U_{1} + U_{2}|\). Here we use L\(\mathcal{'}\) and \(\lambda\mathcal{'}\) to avoid loosing resolution through rounding errors. Then, as for the standard formulae above, we iterate through the equations:

\[\begin{align} & C = \frac{1}{16} f \cos^{2}{\alpha} [4 + f (4-3 \cos^{2}{\alpha}) ] \\ & \cos{2\sigma_{m}} = \cos {\sigma} - 2 \sin{U_{1}} \sin{U_{2}} / \cos^{2}{\alpha} \\ & D = (1-C) f \lbrace \sigma + C \sin{\sigma} [\cos{2\sigma_{m}} + C \cos{\sigma} (-1+2 \cos^{2}{2\sigma_{m}} ] \rbrace \\ & \sin{\alpha} = (L\mathcal{'} - \lambda\mathcal{'}) / D \\ & \sin{\lambda\mathcal{'}} = \sin{\alpha} \sin{\sigma} / (\cos{U_{1}} \cos{U_{2}}) \\ & \sin^{2}{\sigma} = ( \cos{U_{2}} \sin{\lambda\mathcal{'}})^{2} + (\cos{U_{1}} \sin{U_{2}} + \sin{U_{1}} \cos{U_{2}} \cos{\lambda\mathcal{'}})^{2} \\ \end{align}\] The iterations stop when the change in \(\sin{\alpha}\) has converged to the desired accuracy (<1e-12). We then continue:

\[\begin{align} & \sin{\alpha_{1}} = \sin{\alpha} / \cos{U_{1}} \\ & \cos{\alpha_{1}} = \pm (1 - \sin^{2}{\alpha_{1}})^{1/2} \\ \end{align}\] The sign of \(\cos{\alpha_{1}}\) being negative if: \[\begin{align} & \cos{U_{1}} \sin{U_{2}} + \sin{U_{1}} \cos{U_{2}} \cos{\lambda\mathcal{'}} \, < \, 0 \\ \end{align}\] The azimuth at the end point (\(P_{2}\)) is obtained by: \[\tan{\alpha_{2}} = \frac{\sin{\alpha}}{-\sin{U_{1}} \sin{\sigma} + \cos{U_{1}} \cos{\sigma} \cos{\alpha_{1}}}\] To determine the geodesic distance, the below formulas are used: \[\begin{align} & E = (1 + \epsilon \cos^{2}{\alpha} )^{1/2} \\ & F = (E-1) / (E+1) \\ & A = (1 + \frac{1}{4} F^{2}) / (1-F) \\ & B = F (1 - \frac{3}{8} F^{2}) \\ & \Delta\sigma = B \sin{\sigma} \biggl\lbrace \cos{(2\sigma_{m})} + \frac{1}{4} B \biggl[ \cos{\sigma} (-1+2 \cos^{2}{(2\sigma_{m})}) - \frac{1}{6} B \cos{(2\sigma_{m})} (-3+4 \sin^{2}{\sigma}) (-3+4 \cos^{2}{(2\sigma_{m})}) \biggr] \biggr\rbrace \\ \end{align}\] And, finally we calculate the distance: \[s = b A (\sigma - \Delta\sigma = (1 - f) a A (\sigma - \Delta\sigma)\] The last two formulas (\(\Delta\sigma, s\)) are identical to the standard method for distance. At this juncture, I have not attempted to implement these formulae, but will make an effort to incorporate them into the Vincenty inverse program described in the previous post.

So, it seems I have my work cut out for me. That’s all for now. God Bless you and yours. If you don’t know Jesus as your Lord and Savior, there’s no time like the present. And time could be something that may be in short supply, if you wait too long!