FW: Re: calculating the azimuth between two lat/lon's
- To: mathgroup at smc.vnet.net
- Subject: [mg34946] FW: [mg34790] Re: calculating the azimuth between two lat/lon's
- From: "DrBob" <majort at cox-internet.com>
- Date: Fri, 14 Jun 2002 02:38:50 -0400 (EDT)
- Reply-to: <drbob at bigfoot.com>
- Sender: owner-wri-mathgroup at wolfram.com
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
(* copy starting here --- with or without this comment *)
Notebook[{
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[
\(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 is 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]
}, Open ]],
Cell[CellGroupData[{
Cell["Coordinate Functions", "Subsubtitle",
InitializationCell->True],
Cell[BoxData[
\(<< Miscellaneous`CityData`\)], "Input",
InitializationCell->True],
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",
InitializationCell->True],
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",
InitializationCell->True],
Cell[BoxData[{
\(t[{{d_, m_, s_: 0}, {___}}] := \((d + m/60 + s/360)\)
Pi/180\), "\[IndentingNewLine]",
\(t[{{d_, m_: 0}, {___}}] := \((d + m/60)\)
Pi/180\), "\[IndentingNewLine]",
\(t[{{d_: 0}, {___}}] := d\ Pi/180\), "\[IndentingNewLine]",
\(t[name_String] :=
t[CityData[name, CityPosition]]\), "\[IndentingNewLine]",
\(t[name_List] :=
t[CityData[name, CityPosition]]\), "\[IndentingNewLine]",
\(p[{{___}, {d_, m_, s_: 0}}] := \((d + m/60 + s/360)\)
Pi/180\), "\[IndentingNewLine]",
\(p[{{___}, {d_, m_: 0}}] := \((d + m/60)\)
Pi/180\), "\[IndentingNewLine]",
\(p[{{___}, {d_: 0}}] := d\ Pi/180\), "\[IndentingNewLine]",
\(p[name_String] :=
p[CityData[name, CityPosition]]\), "\[IndentingNewLine]",
\(p[name_List] :=
p[CityData[name, CityPosition]]\), "\[IndentingNewLine]",
\(\(t["\<north\>"] = Pi/2;\)\), "\[IndentingNewLine]",
\(\(p["\<north\>"] = 0;\)\), "\[IndentingNewLine]",
\(\(north = v["\<north\>"];\)\)}], "Input",
InitializationCell->True],
Cell[BoxData[{
\(\(Unprotect[ArcTan];\)\), "\[IndentingNewLine]",
\(\(ArcTan[0. , 0. ] = 0;\)\), "\[IndentingNewLine]",
\(\(Protect[ArcTan];\)\)}], "Input"],
Cell[BoxData[{
\(ClearAll[azimuth]\), "\n",
\(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]]\)}], "Input",
InitializationCell->True]
}, Open ]],
Cell[CellGroupData[{
Cell["Examples from City Data", "Subsubtitle",
InitializationCell->True],
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",
\(count =
Length[CityData[CityPosition]]\), "\[IndentingNewLine]",
\(indices =
Union[\(Random[Integer, {1, 258}] &\) /@
Range[40]]\)}], "Input",
InitializationCell->True],
Cell[BoxData[{
\(\(pts =
Sort[\(CityData[CityPosition]\)[\([indices]\)],
OrderedQ[{azimuth[origin, #1],
azimuth[origin, #2]}] &];\)\), "\n",
\(TableForm[\({p[#]/N[Degree], t[#]/N[Degree],
azimuth[origin, #]} &\) /@ pts,
TableHeadings \[Rule] {pts, {"\<Longitude\>", "\<Latitude\>", "\
\<Azimuth\>"}}]\)}], "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}},
Magnification->1.25
]
(* copy ending here --- with or without this comment *)