Geo Stars Library  0.9.3
Geodetic and Astrometry library
GeoStars Library
Author
Dean Nelson (DeanNelson at AOL period Com)
Version
0.9.2 (08/20/11) Added JavaScript routines that mirror C++ library
0.9.1 (07/18/11) Updating all 3rd party packages
0.9 (01/10/10) Updated World Magnetic Model to 2010
0.8 (02/06/08) Cleanup for final release
0.7 Functional Beta (11/26/05) Updates to EFG2LLH Routines
0.6 Functional Beta (01/08/05) First SourceForge Release - additional functions for easy Excel integration
0.5 Functional Beta (01/08/05) Updated World Magnetic Model to 2005
0.4 Functional beta (11/20/04) tuned for DLL and VB support
0.3 Functional beta (11/01/03) with geomagnetism routines
0.2 Functional beta (02/02/02) with test routines
0.1 functional beta (01/28/02)

Home

Geodetic Library routines augmented by astrometry routines

Table of Contents

Introduction

This geodetic library has functions for dealing with many geodesy-based problems found in positioning, pointing, and surveying situations. It is useful to determine absolute position on the earth, pointing vectors, coordinate transformations, and deg/min/sec conversions. It has the following features:

  • ANSI C code
  • Sun position
  • Accurate Azimuth, Elevation, and Range calculation
  • Cartesian to Polar conversions
  • Multiple geocentric to geodetic coordinate conversions
  • It has 23 ellipsoid definitions (can be used worldwide)
  • Links with the Naval Observatory's Novas library for astronomic calculations
  • Calculates the earth's magnetic declination of any location and time(2010-2015)

Compiling Instructions

Geo Library was designed to be as flexible as possible. It does not come as a linkable library. You compile the code on your command line or place it in your own library. Geo Library contains these files:

  • geoStars.h - Main Geo Library include file
  • geoEllips.c - Ellipsoid definitions and location descriptor routines
  • geoPoint.c - Coordinate transformations, polar and Cartesian conversions
  • geoAstro.c - Wrapper routines for the Novas library used in Sun Position
  • geoMag.c - Geo Magnetism routines for the World Magnetism Model of 2010
  • geoEfg2Llh.c - Geocentric to geodetic conversions (7 different methods)
  • geoEfg2Llh_fast.c - Fast and lower accuracy Efg2Llh routines
  • NOVAS Library - This is the Naval Observatory Vector Astrometry Subroutine library.

Test Code

  • geoTest.c - general test routines for the basic geodetic part of the GeoStars Library
  • geoMagTest.c - test routines for the World Magnetism Model of 2010 (WMM-2010)
  • geoTest.bas - VB test routines for geodetic and WMM-2010 functions
  • geoEfg2LlhTest.c - test routines for geocentric to geodetic conversions

Novas could be located in a subdirectory called novas in the main code tree. When using Novas, you will need to compile in the following files:

  • novas.c
  • novascon.c
  • readeph0.c
  • solsys3.c

There are also header files in that directory that you may want to make your compiler aware of. Novas source code is not supplied with GeoStars Library. It can be downloaded from http://aa.usno.navy.mil/software/novas/novas_info.html

Example compilation with GNU C++

with novas already in a library:

gcc geoTest.c geoAstro.c geoEllips.c geoPoint.c geoMag.c geoEfg2Llh.c novas\novas.a -lm -o geoTest

with Novas not in a library:

gcc geoTest.c geoAstro.c geoEllips.c geoPoint.c geoMag.c geoEfg2Llh.c \ novas\novas.c novas\novascon.c novas\readeph0.c novas\solsys3.c novas/nutation.c -lm -o geoTest

Example compilation with Borland C++

with novas already in a library:

bcc32 geoTest.c geoAstro.c geoEllips.c geoPoint.c geoMag.c geoEfg2Llh.c novas\novas.lib

with Novas not in a library:

bcc32 geoTest.c geoAstro.c geoEllips.c geoPoint.c geoMag.c geoEfg2Llh.c \ novas\novas.c novas\novascon.c novas\readeph0.c novas\solsys3.c

Example compilation with Microsoft C++

with novas already in a library:

cl geoTest.c geoAstro.c geoEllips.c geoPoint.c geoMag.c geoEfg2Llh.c novas\novas.lib

with Novas not in a library:

cl geoTest.c geoAstro.c geoEllips.c geoPoint.c geoMag.c geoEfg2Llh.c \ novas\novas.c novas\novascon.c novas\readeph0.c novas\solsys3.c

Calling the geo.dll from Visual Basic

Example:

Private Declare Function MagD Lib "geoStarsLib.dll" Alias "_geoMagGetDecRet@36" _
(ByVal Latitude As Double, ByVal Longitude As Double, ByVal Altitude As Double, _
ByVal Month As Integer, ByVal Day As Integer, ByVal Year As Integer) As Double
' Although the path to "geoStarsLib.dll" could be included, if it is in the system path, the local
' directory, c:\windows or \windows\system32, the path is not required, as it will find it in going
' through the defaults.
Dim dlat As Double 'Latitude as decimal. Plus for Northern Latitudes
Dim dlong As Double 'Longitude as decimal. Minus for West Longitude's
Dim alt As Double 'Altitude in Meters An average value for a state can be used
Dim mo As Integer 'Month
Dim da As Integer 'Day
Dim yr As Integer 'Year
'_______________________________________________________________
Private Sub cmdDeclination_Click() 'They pushed the calculate button'
'Note, the following values would normally be generated by some program steps
'and checked to be sure they were numeric etc.
dlat = 37.1668611
dlong = -93.35138889
alt = 350
mo = 11
da = 2
yr = 2004
txtResult.Text = MagD(dlat, dlong, alt, mo, da, yr) 'The result is displayed in a text box txtResult
End Sub
'_______________________________________________________________

More examples and definintions are included in geo.bas and geoTest.bas.

\ref Home

Formulas

Basic Geodesy

Ellipsoids

Contrary to popular belief, the Earth is not round. It, like most people, is a bit bigger around the girth. That girth is called the equator. The equator is one of the two axis of the Earth. The other axis is the polar axis.

Ellipsoids contain several defining values:

  • Major Axis - This is the equatorial axis or radius $a$
  • Minor Axis - This is the polar axis or radius $b$
  • Flattening - This is a ratio describing the amount of flattening that is present on the ellipsoid. Flattening $f$ is represented by

    \[f = \frac{a - b}{a}\]

    Flattening is usually a small number (i.e. 0.0033900753) and is often expressed as a fraction (i.e. $\frac{1}{297}$). Inverse Flattening $f^{-1}$ is, of course, the denominator of the fraction.
  • Eccentricity - also called Eccentricity Squared $e^2$ is represented by

    \[e^2 = \frac{a^2 - b^2}{a^2}\]

  • Eccentricity Prime - also called Eccentricity Squared Prime $e^2_p$ is represented by

    \[e^2_p = \frac{a^2 - b^2}{b^2}\]

ellipse.gif
Ellipsoid Parameters
The Geo Library contains definitions of the following ellipsoids:

       Ellipsoid Name,                       ID,  Major Axis, Inverse Flattening
  ---------------------------------------------------------------------------
 Airy 1830,                                  AA, 6377563.396, 299.3249646
 Australian National,                        AN, 6378160.0  , 298.25
 Bessel 1841,                                BR, 6377397.155, 299.1528128
 Bessel 1841 (Namibia),                      BN, 6377483.865, 299.1528128
 Clarke 1866,                                CC, 6378206.4  , 294.9786982
 Clarke 1880,                                CD, 6378249.145, 293.465
 Everest (Brunei, E. Malaysia ),             EB, 6377298.556, 300.8017
 Everest 1830,                               EA, 6377276.345, 300.8017
 Everest 1956 (India and Nepal),             EC, 6377301.243, 300.8017
 Everest (Pakistan),                         EF, 6377309.613, 300.8017
 Everest 1948 (W. Malaysia and Singapore),   EE, 6377304.063, 300.8017
 Everest 1969 (W. Malaysia),                 ED, 6377295.664, 300.8017
 Geodetic Reference System 1980,             RF, 6378137.0  , 298.257222101
 Helmert 1906,                               HE, 6378200.0  , 298.3
 Hough 1960,                                 HO, 6378270.0  , 297.0
 Indonesian 1974,                            ID, 6378160.0  , 298.247
 International 1924,                         IN, 6378388.0  , 297.0
 Krassovsky 1940,                            KA, 6378245.0  , 298.3
 Modified Airy,                              AM, 6377340.189, 299.3249646
 Modified Fischer 1960,                      FA, 6378155.0  , 298.3
 South American 1969,                        SA, 6378160.0  , 298.25
 WGS 1972,                                   WD, 6378135.0  , 298.26
 WGS 1984,                                   WE, 6378137.0  , 298.257223563
 NAD 1983,                                   83, 6378137.0   ,298.257222101
\ref Home

Coordinate Systems

There are several coordinate systems used in the Geo Library.

  • Latitude ( $\phi$), Longitude ( $\lambda$), Height ( $h$) This coordinate system references a point on the Earth's ellipsoid. Latitude and longitude are measured in degrees. Height is measured in meters.
latlon2.gif
Latitude, Longitude, and height
  • Earth Fixed Geocentric Coordinates This Cartesian coordinate system uses the center of the earth as a reference point. The coordinates represent an X,Y,Z coordinate system. However, in Geo Library, these coordinates are not called X,Y,Z. They are called E,F,G. This is to help distinguish between the tangential plane X,Y,Z coordinate systems.
ecefxyz.jpg
Earth Fixed Geocentric Coordinates
  • Range ( $r$), Azimuth ( $\alpha$), Elevation ( $\epsilon$) This polar coordinate system originates from a point on the Earth's ellipsoid. Range ( $r$) is measured in meters and represents the slant range from the origin to the point. Azimuth ( $\alpha$) is measured clockwise from north (zero degrees) in degrees. $\epsilon$ is measured from the tangential plane upward. Also in degrees.
rae.gif
Range, Azimuth, and Elevation Coordinates
  • X, Y, Z (also known as ENU) This coordinate system also originates at point on the Earth's ellipsoid. X (also known as East) is measured in meters from the point eastward. Y (also known as North) is measured in meters from the point northward. Z (also known as Up) is measured in meters, in the tangential plane, in an upward direction.
xyz.gif
X,Y,Z Cartesian Coordinates
\ref Home

Conversion Formulas

  • Range $r$, Azimuth $\alpha$, and Elevation $\epsilon$ to X,Y,Z If the slant range $r$, azimuth $\alpha$, and elevation $\epsilon$ from the origin to a point are known, topocentric Cartesian coordinates can be computed in the following manner:

\begin{eqnarray*} X &=& r \; \cos \, \epsilon \; \sin \, \alpha \\ Y &=& r \; \cos \, \epsilon \; \cos \, \alpha \\ Z &=& r \; \sin \, \epsilon \\ \end{eqnarray*}

xyz2aer.gif
Conversion of Range, Azimuth, and Elevation to E,F,G Geocentric Coordinates
  • X,Y,Z to Range $r$, Azimuth $\alpha$, and Elevation $\epsilon$ If Cartesian coordinates X,Y,Z are known, then slant range $r$, azimuth $\alpha$, and elevation $\epsilon$ can be determined this way:

\begin{eqnarray*} r &=& \sqrt{X^2 + Y^2 + Z^2} \\ \alpha &=& \tan^{-1}\left\{\frac{X}{Y}\right\} \left(if \; \alpha < 0, \; then \; \alpha=\alpha +2\pi \right) \\ \epsilon &=& \tan^{-1}\left\{\frac{Z}{ \sqrt{X^2 + Y^2}}\right\} \\ \end{eqnarray*}

  • Latitude $\phi$, Longitude $\lambda$, Height $h$ to E,F,G Geodetic coordinates, Latitude, longitude, and height can be converted to geocentric coordinates through this method:

\begin{eqnarray*} Radius \; of \; curvature: \; N &=&\frac{a}{\sqrt{1-e^2\sin^2\phi}} \\ E &=& \left(N+h\right) \; \cos\phi \; \cos\lambda \\ F &=& \left(N+h\right) \; \cos\phi \; \sin\lambda \\ G &=& \left[N\left(1-e^2\right)+h\right] \; \sin\phi \\ \end{eqnarray*}

  • E,F,G to Latitude $\phi$, Longitude $\lambda$, Height $h$ Geocentric coordinates E,F,G can be converted to latitude $\phi$, longitude $\lambda$, and height $h$. The geoStarsLib has several methods to do this conversion.

Iterative Methods available:

  • Hirvonen & Moritz method
  • Torge method
  • Astronomical Almanac 2002 method
  • Bowring method

Closed Solutions available (non-iterative):

  • Borkowski method
  • Vermeille method
  • geoStarLib method

\[\mathbf{\underline{geoStarsLib \; Method}}\]

\begin{eqnarray*} \lambda &=& \tan^{-1}\left\{\frac{F}{E}\right\} \\ p &=&\sqrt{E^2+F^2} \\ u &=& \tan^{-1}\left\{ \frac{G}{p} \; \frac{a}{b} \right\} \\ \phi &=& \tan^{-1} \left\{ \frac{G+e^2_p \; b \; \sin^3u}{p - e^2 \; a \; \cos^3u} \right\} \\ u_p &=& \tan^{-1} (( 1-f)\tan \; \phi) \\ sign &=& \frac {|p-a \; \cos \; u_p|}{p-a \; \cos \; u_p} \\ h &=& sign \sqrt{(p-a \; \cos \; u_p)^2 + (G-b \; \sin \; u_p)^2} \\ \end{eqnarray*}

\[\mathbf{\underline{Hirvonen \; and \;Moritz \; Method}} \]

\begin{eqnarray*} \lambda &=& \tan^{-1} \left\{ \frac{F}{E} \right\} \\ p &=& \sqrt{E^2+F^2} \\ \phi_1 &=& \tan^{-1} \left\{ \frac{G}{p} \;\; \left(\frac{1}{1-e^2} \right) \right\} \;first \; approximation \\ iterate \; &\phi_1& \; until \; the \; \Delta \phi \; is \;insignificant. \\ N &=& \frac{a}{\sqrt{1-e^2\sin^2\phi_1}} \\ \phi &=& \tan^{-1} \left\{ \frac{G}{p} \left({1 + \frac{e^2 \;N \;\sin(\phi_1)}{G}}\right) \right\} \\ \phi_1 &=& \phi \\ h &=& \frac{p}{\cos\phi_1}-N \\ \end{eqnarray*}

\[\mathbf{\underline{Torge \; Method \;/\; Heiskanen \; and \; Moritz \;Method}} \]

\begin{eqnarray*} \lambda &=& \tan^{-1} \left\{ \frac{F}{E} \right\} \\ p &=& \sqrt{E^2+F^2} \\ prepare \; initial \; estimate \; \; \phi_1 &=& \tan^{-1} \left\{\frac{G}{p} \left(\frac{1}{1-e^2}\right)\right\} \\ iterate \; &\phi_1& \; until \; the \; \Delta \phi \; is \;insignificant. \\ N &=& \frac{a}{\sqrt{1-e^2 \sin^2\phi_1}} \\ h &=& \frac{p}{\cos\phi_1}-N \\ \phi_1 &=& \tan^{-1} \left\{\frac{G}{p} \left(\frac{1}{1-e^2\frac{N}{N+h}}\right)\right\}\\ \phi &=& \phi_1 \\ \end{eqnarray*}

\[\mathbf{\underline{Astronomical \; Almanac \; 2002 \; Method}} \]

\begin{eqnarray*} \lambda &=& \tan^{-1} \left\{ \frac{F}{E} \right\} \\ p &=& \sqrt{E^2+F^2} \\ prepare \; &initial& \; estimate \; \; \\ \phi_1 &=& \tan^{-1} \left\{\frac{G}{p} \right\} \\ iterate \; &\phi_1& \; until \; the \; \Delta \phi \; is \;insignificant. \\ C &=& \frac{1}{\sqrt{1 - e^2 \sin^2 \phi_1}} \\ \phi_1 &=& \tan^{-1} \left\{\frac{G+aCe^2\sin\phi_1}{p}\right\} \\ \phi &=& \phi_1 \\ \end{eqnarray*}

\[\mathbf{\underline{Borkowski \; Method}} \]

\begin{eqnarray*} \lambda &=& \tan^{-1} \left\{ \frac{F}{E} \right\} \\ r &=& \frac {E}{\cos \lambda} \\ e_b &=& \frac{b G - \left( a^2 - b^2\right)}{a r} \\ f_b &=& \frac{b G + \left( a^2 - b^2\right)}{a r} \\ q &=& 2 \left(e_b^2 - f_b^2) \right) \\ p &=& \frac{4}{3} \left(e_b f_b + 1 \right) \\ d &=& p^3 + q^2 \\ v &=& \sqrt[3]{\sqrt{d} - q} - \sqrt[3]{\sqrt{d} + q} \\ g_b &=& \frac{\sqrt{e^2_b + v} + e_b}{2} \\ t &=& \sqrt{g_b^2 + \frac{f_b-vg_b}{2g_b-e_b} - g_b} \\ \phi &=& \tan^-1\left(\frac{a\left(1-t^2\right)}{2bt} \right) \\ h &=& \left(r-at\right)\cos\phi + \left(G-b \right)\sin\phi \\ \end{eqnarray*}

\[\mathbf{\underline{Vermeille \; Method}} \]

\begin{eqnarray*} p &=& \frac{E^2 +F^2}{a2} \\ q &=& G^2\frac{1-e^2}{a^2} \\ r &=& \frac{p+q-e^4}{6} \\ s &=& e^4 \frac{pq}{4r^3} \\ t &=& \sqrt[3]{1+s+\sqrt{s\left(2+s\right)}} \\ u &=& r\left(1+t+\frac{1}{t} \right) \\ v &=& \sqrt{u^2+e^4q} \\ w &=& e^2\frac{u+v-q}{2v} \\ k &=& \sqrt{u+v+w^2}-w \\ D &=& \frac{k\sqrt{E^2 +G^2}}{k+e^2} \\ \lambda &=& 2 \;tan^{-1} \frac{F}{E+\sqrt{E^2+F^2}} \\ \phi &=& 2 \;tan^{-1} \frac{G}{D+\sqrt{D^2+G^2}} \\ h &=& \frac{k+e^2-1}{k}\sqrt{D^2+G^2} \\ \end{eqnarray*}

Geomagnetic Routines

The modeling of the earth's magnetic field is called geomagnetism. The most current application of this model is magnetic declination - deflection of the compass needle off true north. The earth's magnetic field is in constant motion. Thus, a stable model needed to be developed so that one might be able to correct for the changes of the field in time.

By using the geomagnetic routines in this library, true north can be determined quite easily.

The geoStarsLib library uses the World Magnetic Model of 2010.

The following figure shows the magentic declination in all areas of the world.

WMM2010_D_MERC_reduced.png
WMM 2010 declination

Since the magnetic field is always in flux, the WMM2010 discusses the annual rate of change and charts it as follows:

WMM2010_D_SV_MERC_reduced.png
WMM 2010 rate of change in declination

Using Visual Basic and Excel with the geoStarslib DLL

geoStarsLib.dll can be used with VB, VBA, and other programs that can access external DLL's. The DLL must be located in the path to be used with Excel. The file geoTest.bas is an example of how the DLL functions are called from VB. It tests the sun position at the current time. You must make sure that geo.bas is loaded in your project to provide the interface to the DLL.

All of the geoStarsLib routines can be called from Visual Basic using the Function and Call methods.

Follow this link for more information:

\ref Home

How do I ....

Here are some common questions:

  1. How do I convert Range, Azimuth, and Elevation to XYZ coordinates?
  2. How do I convert XYZ coordinates to Range, Azimuth, and Elevation?
  3. How do I get a Range, Azimuth, and Elevation from one point to another and all I have are their geodetic coordinates?
    • use geoInitLocation() to establish a base location and then again to establish the point location
    • use geoEfg2XyzDiff() to get the XYZ coordinates of the point (with the base as the origin).
    • use geoXyz2Rae() to get your polar coordinates

  4. How do I get X,Y,Z from one point to another and all I have are their geodetic coordinates?
    • use geoInitLocation() to establish a base location and then again to establish the point location
    • use geoEfg2XyzDiff() to get the XYZ coordinates of the point (with the base as the origin).

  5. How do I get earth centered coordinates (EFG) from Range, Azimuth, and Elevation?
  6. How do I get the Latitude, Longitude and height of the point indicated by Range, Azimuth, and Elevation?
  7. How do I convert radians to degrees, minutes, and seconds?
  8. How do I convert degrees, minutes, and seconds to radians?
  9. How do I find out the magnetic declination given a Latitude and Longitude?
    • use geoMag() to get every blessed geomagnetic variable under the sun, or...
    • use geoMagGetDec() to get just the magnetic declination, or ...
    • use geoMagGetDecRet() to return a double precision value declination (also used in VB)

  10. How do I get the position of the sun right now?
  11. How do I get the sun position for any day?
  12. How do I get the sun position for any day when I only have a Julian date?
  13. How do I get a date with this Julian?
    • get a life!
\ref Home