MathGroup Archive 2002

[Date Index] [Thread Index] [Author Index]

Search the Archive

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
> 
> ...................






  • Prev by Date: Re: Re: Functionality and Reliability
  • Next by Date: StoppingTest option for NDSolve
  • Previous by thread: Re: calculating the azimuth between two lat/lon's
  • Next by thread: Re: calculating the azimuth between two lat/lon's