<html><body name="Mail Message Editor"><div><br></div><div>Mikæl:</div><div><br></div><div>The usual spherical formulæ for computing a point at a known distance and azimuth are:<br></div><div><br></div><div><span class="Apple-style-span" style="font-family: Monaco;">φ₂ = arcsin (sin φ₁ cos c + cos φ₁ sin c cos Az)<br></span></div><div><span class="Apple-style-span" style="font-family: Monaco;">λ₂ = λ₁ + atan2 [sin c sin Az, (cos φ₁ cos c - sin φ₁ sin c cos Az)]</span></div><div><br></div><div>where the start point is subscripted with 1, the end point is subscripted with 2, and the arc distance is c. We know φ₂, λ₂, Az, and c. Therefore these reduce to</div><div><br></div><div><div><span class="Apple-style-span" style="font-family: Monaco;">(1) n₄ = n₁ sin φ₁ + n₂ cos φ₁</span></div><div><span class="Apple-style-span" style="font-family: Monaco;">(2) λ₂ = λ₁ + atan2 [n₃, (n₁ cos φ₁ - n₂ sin φ₁)]</span></div><div><br></div><div>where we have made substitutions for symbolic clarity:</div><div><br></div><div><span class="Apple-style-span" style="font-family: Monaco;">n₁ = cos c</span></div><div><span class="Apple-style-span" style="font-family: Monaco;">n₂ = sin c cos Az</span></div><div><span class="Apple-style-span" style="font-family: Monaco;">n₃ = sin c sin Az</span></div><div><span class="Apple-style-span" style="font-family: Monaco; ">n₄ = sin φ₂</span><br></div><div><br></div><div>Since equation 1 <span class="Apple-style-span" style="font-family: Monaco; "><span class="Apple-style-span" style="font-family: Helvetica;"><span class="Apple-style-span" style="background-color: transparent;">ha</span></span><span class="Apple-style-span" style="font-family: Helvetica; ">s only φ₁ unknown, we solve with trigonometric and algebraic manipulations:</span></span></div></div><div><br></div><div><span class="Apple-style-span" style="font-family: Monaco;">tan φ₁ = A / B</span></div><div><br></div><div>with φ₁ resolvable using the usual atan2 function in order to preserve quadrants. Here</div><div><br></div><div><span class="Apple-style-span" style="font-family: Monaco;"> n₄ - n₂ B</span></div><div><span class="Apple-style-span" style="font-family: Monaco;">A = _________</span></div><div><span class="Apple-style-span" style="font-family: Monaco;"> n₁</span></div><div><span class="Apple-style-span" style="font-family: Monaco;"><br></span></div><div><span class="Apple-style-span" style="font-family: Monaco;"> n₂ n₄ ± √(n₂² n₁² + n₁⁴ - n₁² n₄²)</span></div><div><span class="Apple-style-span" style="font-family: Monaco;">B = ----------------------------------</span></div><div><span class="Apple-style-span" style="font-family: Monaco;"> n₁² + n₂²</span></div><div><br></div><div><br></div><div>yielding two solutions, depending on the sign chosen for the radical in B. Both are valid solutions.</div><div><br></div><div>Thus armed with φ₁, we can easily compute<span class="Apple-style-span" style="background-color: transparent;"> </span><span class="Apple-style-span" style="font-family: Monaco; "><span class="Apple-style-span" style="font-family: Helvetica;"><span class="Apple-style-span" style="background-color: transparent;">λ₁ </span></span><span class="Apple-style-span" style="font-family: Helvetica; ">as</span></span></div><div><br></div><div><div><span class="Apple-style-span" style="font-family: Monaco; ">λ₁ = λ₂ - atan2 [n₃, (n₁ cos φ₁ - n₂ sin φ₁)]</span></div><div><span class="Apple-style-span" style="font-family: Monaco;"><br></span></div></div><div>The ellipsoidal case is much harder. Your iterative process will work, of course, though it will be terribly inefficient and likely comes with problematic regions. In principle it seems like the problem you are trying to solve is intermediate between the forward and the inverse problems. The reason it is not as simple as the forward problem is that, in the forward case, you know both the starting point and the azimuth. This wholly characterizes which geodesic you are dealing with, since a geodesic is completely characterized by its constant and any point along it:</div><div><span class="Apple-style-span" style="font-family: Monaco;"><br></span></div><div><span class="Apple-style-span" style="font-family: Monaco;"> cos φ sin Az</span></div><div><span class="Apple-style-span" style="font-family: Monaco;">K = ________________</span></div><div><span class="Apple-style-span" style="font-family: Monaco;"> √(1 - e² sin² φ)</span></div><div><br></div><div>In your case, you cannot compute the geodesic's constant because you don't know both the latitude and azimuth from any known point; instead you know the azimuth from an unknown point. On the other hand, you already know the distance, the calculation of which is what makes the inverse problem so difficult.</div><div><br></div><div>Warning: Do NOT use Pearson's analysis of geodesics (Frederick Pearson II 1984: "Map Projection Methods", pp. 81–82, and 1990: "Map Projections: Theory and Applications", pp. 80–81). He makes an error early on in the derivation, making the final, easily computed integral nothing but a fiction. Sadly, Donald Fenna has repeated this error in his 2007 volume, "A Compendium of Map Projections, with Derivations", pp. 394–395. </div><div><br></div><div>Regards,</div><div>-- daan Strebe</div><div><br></div><br>On Oct 20, 2008, at 12:36:56 AM, "Mikael Rittri" <Mikael.Rittri@carmenta.com> wrote:<br><blockquote style="padding-left: 5px; margin-left: 5px; border-left-width: 2px; border-left-style: solid; border-left-color: blue; color: blue; "><span class="Apple-style-span" style="border-collapse: separate; color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: 2; text-align: auto; text-indent: 0px; text-transform: none; white-space: normal; widows: 2; word-spacing: 0px; -webkit-border-horizontal-spacing: 0px; -webkit-border-vertical-spacing: 0px; -webkit-text-decorations-in-effect: none; -webkit-text-size-adjust: auto; -webkit-text-stroke-width: 0; "><div style="width: 100%; "><div id="felix-mail-header-block" style="color: black; background-color: white; border-bottom-width: 1px; border-bottom-style: solid; border-bottom-color: silver; padding-bottom: 1em; margin-bottom: 1em; width: 100%; "><table border="0" cellpadding="1" cellspacing="1" width="100%"><tbody><tr><td width="70px" style="font-family: 'Lucida Grande'; font-size: 8pt; color: gray; text-align: right; vertical-align: top; font-weight: bold; "><span>From:</span></td><td style="font-family: 'Lucida Grande'; font-size: 8pt; color: black; text-align: left; vertical-align: top; padding-left: 5px; "><span title=""Mikael Rittri" <Mikael.Rittri@carmenta.com>">"Mikael Rittri" <Mikael.Rittri@carmenta.com></span></td></tr><tr><td width="70px" style="font-family: 'Lucida Grande'; font-size: 8pt; color: gray; text-align: right; vertical-align: top; font-weight: bold; "><span>Subject:</span></td><td style="font-family: 'Lucida Grande'; font-size: 8pt; color: black; text-align: left; vertical-align: top; padding-left: 5px; "><span style="font-weight: bold; ">[Proj] Locate a point from distance and backwards azimuth?</span></td></tr><tr><td width="70px" style="font-family: 'Lucida Grande'; font-size: 8pt; color: gray; text-align: right; vertical-align: top; font-weight: bold; "><span>Date:</span></td><td style="font-family: 'Lucida Grande'; font-size: 8pt; color: black; text-align: left; vertical-align: top; padding-left: 5px; "><span>October 20, 2008 12:36:56 AM PDT</span></td></tr><tr><td width="70px" style="font-family: 'Lucida Grande'; font-size: 8pt; color: gray; text-align: right; vertical-align: top; font-weight: bold; "><span>To:</span></td><td style="font-family: 'Lucida Grande'; font-size: 8pt; color: black; text-align: left; vertical-align: top; padding-left: 5px; "><span title="proj@lists.maptools.org">proj@lists.maptools.org</span></td></tr></tbody></table></div><div id="felix-mail-content-block" style="color: black; background-color: white; width: 100%; "><div style="font-family: monospace; color: black; background-color: white; font-size: 8pt; ">Hello,<br>Here is a variation of the first/principal/forward geodetic problem:<br><br>Known:<span class="Apple-converted-space"> </span><br>Position of point C.<br>Distance from point C to an unknown point A.<span class="Apple-converted-space"> </span><br>Azimuth, at A, of the great-circle arc between A and C.<span class="Apple-converted-space"> </span><br>Where is A?<span class="Apple-converted-space"> </span><br><br>Note the difference from the first/principal/forward problem:<span class="Apple-converted-space"> </span><br>the azimuth is known at the unknown point A, instead of at the<span class="Apple-converted-space"> </span><br>known point C.<span class="Apple-converted-space"> </span><br><br>I am not really sure when this is useful, but I have been asked<span class="Apple-converted-space"> </span><br>about it twice by different people, so I feel I ought to solve it<span class="Apple-converted-space"> </span><br>some day.<span class="Apple-converted-space"> </span><br><br>Does this problem has a name, and are there detailed published solutions<br>anywhere?<span class="Apple-converted-space"> </span><br><br>Actually, I should say that I am not completely lost:<span class="Apple-converted-space"> </span><br>For ellipsoid accuracy, I would use Vincenty's formulas for the<span class="Apple-converted-space"> </span><br>forward problem starting from C, and then use some numerical root-finding<span class="Apple-converted-space"> </span><br>algorithm to find the azimuth at C that gives the desired azimuth at A.<span class="Apple-converted-space"> </span><br>For spherical accuracy, one can add the north pole as a third point B,<span class="Apple-converted-space"> </span><br>and then ABC is a spherical triangle where two sides (a and b) and one<span class="Apple-converted-space"> </span><br>angle (A) is known. For this case, textbooks usually suggest that one uses<span class="Apple-converted-space"> </span><br>the spherical law of sines to find sin(B) (where B is the difference in<span class="Apple-converted-space"> </span><br>longitude between point A and point C), and then one of Napier's analogies to<span class="Apple-converted-space"> </span><br>find the angle C and the side c. The most detailed formulas that I have found<br>are at http://encarta.msn.com/encyclopedia_761572350_2/Trigonometry.html#s7 ,<span class="Apple-converted-space"> </span><br>but the formula font makes it impossible to distinguish the upper case C<span class="Apple-converted-space"> </span><br>from the lower case c. I think it must be a lower case c in the 2nd formula<br>and a an upper case C in the 3rd formula, though. I also think there is a<span class="Apple-converted-space"> </span><br>typo in the 3rd formula: (a+b) and (a-b) should both be multiplied by 1/2,<span class="Apple-converted-space"> </span><br>before the sine is taken.<br><br>But I am not quite sure about the fine details. If I know sin(B),<span class="Apple-converted-space"> </span><br>I still don't known whether the angle B is acute or obtuse, or if both<span class="Apple-converted-space"> </span><br>are possible solutions. I think that if the distance between A and C is<span class="Apple-converted-space"> </span><br>larger than the distance from C to the nearest pole, there would be either<span class="Apple-converted-space"> </span><br>two solutions or none. (Well, of course, if C is the north pole and the<span class="Apple-converted-space"> </span><br>azimuth at A should be zero, there would be infinitely many solutions.)<br><br>So, if anyone recognizes this problem, I'd be grateful for any advice<span class="Apple-converted-space"> </span><br>or literature references.<span class="Apple-converted-space"> </span><br><br>Best regards,<br><br>--<br>Mikael Rittri<br>Carmenta AB<br>Box 11354<br>SE-404 28 Göteborg<br>Visitors: Sankt Eriksgatan 5<br>SWEDEN<br>Tel: +46-31-775 57 37<br>Mob: +46-703-60 34 07<span class="Apple-converted-space"> </span><br>mikael.rittri@carmenta.com<br>www.carmenta.com<br><br>_______________________________________________<br>Proj mailing list<br>Proj@lists.maptools.org<br>http://lists.maptools.org/mailman/listinfo/proj<br><br></div></div></div></span></blockquote><br><div><br></div><div class="aol_ad_footer" id="u6A4A039FB93E49A38C55DE23FBA50F2B"></div></body></html>