RE: Re: calculating the azimuth between two lat/lon's
- To: mathgroup at smc.vnet.net
- Subject: [mg35009] RE: [mg34984] Re: calculating the azimuth between two lat/lon's
- From: "DrBob" <majort at cox-internet.com>
- Date: Wed, 19 Jun 2002 05:52:45 -0400 (EDT)
- Reply-to: <drbob at bigfoot.com>
- Sender: owner-wri-mathgroup at wolfram.com
Here's an updated notebook expression for the azimuth problem. I've added calculation of distance between points on the Earth's surface, as well as azimuth (course bearing). A spherical approximation is used in both cases. Copy the expression into a notebook and then evaluate the notebook. Bobby Treat (* start copying here *) Notebook[{ Cell[BoxData[{ \(Needs["\<Miscellaneous`CityData`\>"]\), "\n", \(Needs["\<Calculus`VectorAnalysis`\>"]\), "\[IndentingNewLine]", \ \(Needs["\<LinearAlgebra`Orthogonalization`\>"]\)}], "Input", InitializationCell->True], Cell[CellGroupData[{ Cell["Derivation", "Subsubtitle"], Cell["\<\ Using the Law of Sines and Law of Cosines for spherical triangles, we \ can solve for the azimuth from ptA to ptB. First we define the a, b, c, A, B, and C as follows:\ \>", "Text"], Cell[BoxData[{ \(Clear[a, b, c, A, B]\), "\[IndentingNewLine]", \(info = TableForm[{{a, "\<arc angle, ptB to north \ (90\[Degree]-latitude)\>"}, {A, "\<angle at ptA (Azimuth)\>"}, {b, \ "\<arc angle, ptA to north (90\[Degree]-latitude)\>"}, {B, "\<angle \ at ptB\>"}, {c, "\<arc angle, ptA to ptB\>"}, {C, "\<angle at \ north\>"}}]\)}], "Input"], Cell["\<\ The Cosine of an arc angle is the dot product of vectors terminating \ the arc, with each angle in the interval {0, \[Pi]}:\ \>", "Text"], Cell[BoxData[ \(c \[Equal] ArcCos[ptA . ptB]\)], "Input", Evaluatable->False], Cell["Angle C is the difference in longitudes at ptA and ptB:", \ "Text"], Cell[BoxData[ \(C \[Equal] ArcCos[Cos[p[a] - p[b]]]\)], "Input", Evaluatable->False], Cell["The Law of Sines is:", "Text"], Cell[BoxData[ \(Sin[a]\/Sin[A] \[Equal] Sin[b]\/Sin[B] \[Equal] Sin[c]\/Sin[C]\)], "Input", Evaluatable->False], Cell["Hence the azimuth A satisfies", "Text"], Cell[BoxData[ \(\(Sin[A] \[Equal] \(Sin[a] Sin[C]\)\/Sin[c];\)\)], "Input", Evaluatable->False], Cell["The Law of Cosines at ptA tells us:", "Text"], Cell[BoxData[ \(Cos[a] \[Equal] Cos[b] Cos[c] + Sin[b] Sin[c] Cos[A]\)], "Input", Evaluatable->False], Cell["from which we can derive", "Text"], Cell[BoxData[ \(Cos[ A] \[Equal] \(Cos[a] - Cos[b] Cos[c]\)\/\(Sin[b] Sin[c]\)\)], \ "Input", Evaluatable->False], Cell["Hence A is determined by", "Text"], Cell[BoxData[{ \(A \[Equal] ArcTan[Cos[A], Sin[A]]\), "\[IndentingNewLine]", \(\ \ \ \ \(\(\[Equal]\)\(ArcTan[\(Cos[a] - Cos[b] \ Cos[c]\)\/\(Sin[b] Sin[c]\), \(Sin[a] Sin[C]\)\/Sin[c]]\)\)\), "\ \[IndentingNewLine]", \(\ \ \ \ \(\(\[Equal]\)\(ArcTan[Cos[a] - Cos[b] Cos[c], Sin[a] Sin[b], Sin[C]]\)\)\)}], "Input", Evaluatable->False], Cell[TextData[{ "The distance from ", StyleBox["ptA", FontWeight->"Bold"], " to ", StyleBox["ptB", FontWeight->"Bold"], " is the spherical radius ", StyleBox["r", FontWeight->"Bold"], " times the arc angle ", StyleBox["c", FontWeight->"Bold"], ":" }], "Text"], Cell[BoxData[ \(r\ ArcCos[ptA . ptB]\)], "Input", Evaluatable->False], Cell["\<\ For the Earth's surface, the radius in miles is approximately\ \>", "Text"], Cell[BoxData[ \(r = 3959.74\)], "Input", Evaluatable->False] }, Open ]], Cell[CellGroupData[{ Cell["Coordinate Functions", "Subsubtitle"], Cell["\<\ The following computes a vector in three dimensions from its \ spherical coordinates. t[a] is latitude and p[a] is longitude. All \ input angles are in radians.\ \>", "Text"], Cell[BoxData[ \(v[a_] := {Cos[t[a]]\ Cos[p[a]], Cos[t[a]]\ Sin[p[a]], Sin[t[a]]}\)], "Input", InitializationCell->True], Cell["\<\ The functions p and t can be computed from latitude and longitude \ data in the CityData format as follows:\ \>", "Text"], Cell[BoxData[{ \(ClearAll[p, t, north, s]\), "\[IndentingNewLine]", \(t[{{d_, m_, s_: 0}, {___}}] := \((d + m/60 + s/360)\) Pi/180.0\), "\[IndentingNewLine]", \(t[{{d_, m_: 0}, {___}}] := \((d + m/60)\) Pi/180.0\), "\[IndentingNewLine]", \(t[{{d_: 0}, {___}}] := d\ Pi/180.0\), "\[IndentingNewLine]", \(t[name_String] := t[CityData[name, CityPosition]]\), "\[IndentingNewLine]", \(t[{name__String}] := t[CityData[{name}, CityPosition]]\), "\[IndentingNewLine]", \(p[{{___}, {d_, m_, s_: 0}}] := \((d + m/60 + s/360)\) Pi/180.0\), "\[IndentingNewLine]", \(p[{{___}, {d_, m_: 0}}] := \((d + m/60)\) Pi/180.0\), "\[IndentingNewLine]", \(p[{{___}, {d_: 0}}] := d\ Pi/180.0\), "\[IndentingNewLine]", \(p[name_String] := p[CityData[name, CityPosition]]\), "\[IndentingNewLine]", \(p[{name__String}] := p[CityData[{name}, CityPosition]]\), "\[IndentingNewLine]", \(\(t["\<north\>"] = Pi/2.0;\)\), "\[IndentingNewLine]", \(\(p["\<north\>"] = 0.0;\)\), "\[IndentingNewLine]", \(\(v["\<north\>"] = \(north = {0, 0, 1}\);\)\), "\[IndentingNewLine]", \(\(north = v["\<north\>"];\)\), "\[IndentingNewLine]", \(p[{a_, b_, c_}] := N[ArcTan[a, b]]\), "\[IndentingNewLine]", \(t[{a_, b_, c_}] := N[ArcSin[c]]\), "\[IndentingNewLine]", \(s[name_String] := {1, Pi/2 - t[name], p[name]}\), "\[IndentingNewLine]", \(s[name : {__String}] := {1, Pi/2 - t[name], p[name]}\), "\[IndentingNewLine]", \(s[A_, B_, \[Theta]_] := v[A]\ Cos[\[Theta]] + basis[A, B] Sin[\[Theta]]\), "\[IndentingNewLine]", \(cityName[a_String] := a\), "\[IndentingNewLine]", \(cityName[{a__String}] := First[{a}]\)}], "Input", InitializationCell->True], Cell[BoxData[{ \(\(Unprotect[ArcTan];\)\), "\[IndentingNewLine]", \(\(ArcTan[0. , 0. ] = 0;\)\), "\[IndentingNewLine]", \(\(ArcTan[0, 0] = 0;\)\), "\[IndentingNewLine]", \(\(Protect[ArcTan];\)\)}], "Input", InitializationCell->True], Cell[BoxData[{ \(ClearAll[azimuth, distance]\), "\[IndentingNewLine]", \(azimuth[A_, B_] := Module[{ptA, ptB, a, b, c, CC}, \[IndentingNewLine]ptA = v[A]; \[IndentingNewLine]ptB = v[B]; \[IndentingNewLine]a = Pi/2 - t[B]; \[IndentingNewLine]b = Pi/2 - t[A]; \[IndentingNewLine]c = ArcCos[ptA . ptB]; \[IndentingNewLine]CC = Mod[p[B] - p[A], 2 Pi, \(-Pi\)]; \[IndentingNewLine]Mod[ ArcTan[\((north . ptB)\) - \((north . ptA)\) \((ptA . ptB)\), Sin[a] Sin[b] Sin[CC]], 2. Pi]/ N[Degree]\[IndentingNewLine]]\), "\[IndentingNewLine]", \(azimuth[A_, Indeterminate] := 0\), "\[IndentingNewLine]", \(azimuth[Indeterminate, A_] := 0\), "\[IndentingNewLine]", \(distance[A_, B_] := 3959.74\ ArcCos[v[A] . v[B]]\)}], "Input", InitializationCell->True] }, Open ]], Cell[CellGroupData[{ Cell["Examples from City Data", "Subsubtitle", InitializationCell->True], Cell[TextData[{ "From the ", StyleBox["VNR Concise Encyclopedia of Mathematics, 2nd Edition", FontWeight->"Bold"], ", page 273, we find that a course from Leningrad (St. Peterburg) \ to San Francisco begins on a heading of N 21.61\[Degree] W and ends \ on a heading of S 13.52\[Degree] W. The distance is given as 5511.5 \ miles." }], "Text"], Cell[BoxData[ \(360 - azimuth["\<St. Peterburg\>", "\<San Francisco\>"]\)], "Input"], Cell[BoxData[ \(azimuth["\<San Francisco\>", "\<St. Peterburg\>"]\)], "Input"], Cell[BoxData[ \(distance["\<St. Peterburg\>", "\<San Francisco\>"]\)], \ "Input"], Cell["\<\ There are 258 cities in the CityData database. To save time on \ slower computers, we'll take a sample of them:\ \>", "Text"], Cell[BoxData[{ \(\(origin = "\<Dallas\>";\)\), "\n", \(\(allCities = CityData[CityPosition];\)\), "\[IndentingNewLine]", \(count = Length[allCities]\), "\[IndentingNewLine]", \(indices = Union[\(Random[Integer, {1, 258}] &\) /@ Range[40], {First[ First[Position[allCities, origin]]]}]\)}], "Input", InitializationCell->True], Cell["\<\ The following give azimuth and distance from Dallas to various \ cities:\ \>", "Text"], Cell[BoxData[{ \(\(pts = Sort[allCities[\([indices]\)], OrderedQ[{azimuth[origin, #1], azimuth[origin, #2]}] &];\)\), "\n", \(TableForm[\({cityName[#], p[#]/N[Degree], t[#]/N[Degree], azimuth[origin, #], distance[origin, #]} &\) /@ pts, TableHeadings \[Rule] {None, {"\<City\>", "\<Longitude\>", \ "\<Latitude\>", "\<Azimuth\>", "\<Miles\>"}}]\)}], "Input", InitializationCell->True] }, Open ]] }, FrontEndVersion->"4.1 for Microsoft Windows", ScreenRectangle->{{0, 1024}, {0, 711}}, AutoGeneratedPackage->Automatic, WindowSize->{1016, 673}, WindowMargins->{{0, Automatic}, {Automatic, 0}}, PrintingCopies->1, PrintingPageRange->{Automatic, Automatic}, ShowCellLabel->False, Magnification->1.25 ] (* Stop copying here *) -----Original Message----- From: rainer gruber [mailto:rainer.gruber at gmx.at] To: mathgroup at smc.vnet.net Subject: [mg35009] [mg34984] Re: calculating the azimuth between two lat/lon's Sorry, my fault, I "redefined" azimuth. As I found on <http://aa.usno.navy.mil/data/docs/AltAz.html#Notes> "_Azimuth_ is the angle along the horizon, with zero degrees corresponding to North, and increasing in a clockwise fashion". I posted the angle between the two lines from both points to the center of the earth. Rainer Gruber DrBob wrote: > I think azimuth is far more complicated than rainer indicated. Copy the > following notebook expression into an empty notebook, for a discussion > of calculation of azimuth, plus examples using the CityData package. > > If I've missed the boat somewhere, I'm sure someone will let me know!! > > Bobby Treat > > ...................