mgrs util

.Rbuildignore

@@ -0,0 +1,10 @@
+^.*\.Rproj$
+^\.Rproj\.user$
+^\.travis\.yml$
+^README\.*Rmd$
+^README\.*html$
+^NOTES\.*Rmd$
+^NOTES\.*html$
+^\.codecov\.yml$
+^README_files$
+^doc$

.codecov.yml

@@ -0,0 +1 @@
+comment: false

.gitignore

@@ -0,0 +1,8 @@
+.DS_Store
+.Rproj.user
+.Rhistory
+.RData
+.Rproj
+src/*.o
+src/*.so
+src/*.dll

.travis.yml

@@ -0,0 +1,31 @@
+language: r
+
+warnings_are_errors: true
+
+sudo: required
+
+cache: packages
+
+r:
+ - oldrel
+ - release
+ - devel
+
+apt_packages:
+  - libv8-dev
+  - xclip
+
+env:
+ global:
+   - CRAN: http://cran.rstudio.com
+
+after_success:
+  - Rscript -e 'covr::codecov()'
+
+notifications:
+  email:
+    - bob@rud.is
+  irc:
+    channels:
+      - "104.236.112.222#builds"
+    nick: travisci

DESCRIPTION

@@ -0,0 +1,20 @@
+Package: mgrs
+Type: Package
+Title: Convert 'MGRS' Coordinates To and From Other Coordiante Systems
+Version: 0.1.0
+Date: 2017-04-09
+Author: Bob Rudis (bob@rud.is), Howard Butler
+Maintainer: Bob Rudis <bob@rud.is>
+Description: Convert 'MGRS' Coordinates To and From Other Coordiante Systems
+URL: https://github.com/hrbrmstr/mgrs
+BugReports: https://github.com/hrbrmstr/mgrs/issues
+License: AGPL
+Suggests:
+    testthat,
+    covr
+Depends:
+    R (>= 3.2.0)
+Imports:
+    Rcpp
+LinkingTo: Rcpp
+RoxygenNote: 6.0.1

NAMESPACE

@@ -0,0 +1,6 @@
+# Generated by roxygen2: do not edit by hand
+
+export(latlng_to_mgrs)
+export(mgrs_to_latlng)
+importFrom(Rcpp,sourceCpp)
+useDynLib(mgrs)

NEWS.md

@@ -0,0 +1,2 @@
+0.1.0 
+* Initial release

R/RcppExports.R

@@ -0,0 +1,28 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#' Convert an MGRS string to latitude/longitude
+#'
+#' @md
+#' @param MGRS an MGRS string
+#' @param degrees convert to degrees? Default: `TRUE`
+#' @export
+#' @examples
+#' mgrs_to_latlng("15TWG0000049776")
+mgrs_to_latlng <- function(MGRS, degrees = TRUE) {
+    .Call('mgrs_mgrs_to_latlng', PACKAGE = 'mgrs', MGRS, degrees)
+}
+
+#' Convert latitude/longitude to MGRS string
+#'
+#' @md
+#' @param latitude,longitude coordinates
+#' @param degrees are latitude/longitude in degrees? Default: `TRUE`
+#' @param precision level of precision for the conversion. Default `5`
+#' @export
+#' @examples
+#' latlng_to_mgrs(42, -93)
+latlng_to_mgrs <- function(latitude, longitude, degrees = TRUE, precision = 5L) {
+    .Call('mgrs_latlng_to_mgrs', PACKAGE = 'mgrs', latitude, longitude, degrees, precision)
+}
+

R/mgrs-package.R

@@ -0,0 +1,8 @@
+#' Convert 'MGRS' Coordinates To and From Other Coordiante Systems
+#'
+#' @name mgrs
+#' @docType package
+#' @author Bob Rudis (bob@@rud.is)
+#' @useDynLib mgrs
+#' @importFrom Rcpp sourceCpp
+NULL

README.Rmd

@@ -0,0 +1,43 @@
+---
+output: rmarkdown::github_document
+---
+
+`mgrs` : Convert 'MGRS' Coordinates To and From Other Coordiante Systems
+
+NOTE: These will be vectorized soon.
+
+The following functions are implemented:
+
+- `latlng_to_mgrs`:	Convert latitude/longitude to MGRS string
+- `mgrs_to_latlng`:	Convert an MGRS string to latitude/longitude
+
+### Installation
+
+```{r eval=FALSE}
+devtools::install_github("hrbrmstr/mgrs")
+```
+
+```{r message=FALSE, warning=FALSE, error=FALSE, include=FALSE}
+options(width=120)
+```
+
+### Usage
+
+```{r message=FALSE, warning=FALSE, error=FALSE}
+library(mgrs)
+
+# current verison
+packageVersion("mgrs")
+
+mgrs_to_latlng("33UXP04")
+
+latlng_to_mgrs(48.20535, 16.34593)
+
+mgrs_to_latlng("33UXP0500444996")
+
+latlng_to_mgrs(48.24948, 16.41449)
+
+mgrs_to_latlng("24XWT783908")
+
+latlng_to_mgrs(83.62738, -32.66879)
+```

README.md

@@ -0,0 +1,65 @@
+
+`mgrs` : Convert 'MGRS' Coordinates To and From Other Coordiante Systems
+
+NOTE: These will be vectorized soon.
+
+The following functions are implemented:
+
+-   `latlng_to_mgrs`: Convert latitude/longitude to MGRS string
+-   `mgrs_to_latlng`: Convert an MGRS string to latitude/longitude
+
+### Installation
+
+``` r
+devtools::install_github("hrbrmstr/mgrs")
+```
+
+### Usage
+
+``` r
+library(mgrs)
+
+# current verison
+packageVersion("mgrs")
+```
+
+    ## [1] '0.1.0'
+
+``` r
+mgrs_to_latlng("33UXP04")
+```
+
+    ##      lat      lng 
+    ## 48.20535 16.34593
+
+``` r
+latlng_to_mgrs(48.20535, 16.34593)
+```
+
+    ## [1] "33UXP0000040000"
+
+``` r
+mgrs_to_latlng("33UXP0500444996")
+```
+
+    ##      lat      lng 
+    ## 48.24947 16.41449
+
+``` r
+latlng_to_mgrs(48.24948, 16.41449)
+```
+
+    ## [1] "33UXP0500344996"
+
+``` r
+mgrs_to_latlng("24XWT783908")
+```
+
+    ##       lat       lng 
+    ##  83.62738 -32.66879
+
+``` r
+latlng_to_mgrs(83.62738, -32.66879)
+```
+
+    ## [1] "25XEN0410486507"

man/latlng_to_mgrs.Rd

@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{latlng_to_mgrs}
+\alias{latlng_to_mgrs}
+\title{Convert latitude/longitude to MGRS string}
+\usage{
+latlng_to_mgrs(latitude, longitude, degrees = TRUE, precision = 5L)
+}
+\arguments{
+\item{latitude, longitude}{coordinates}
+
+\item{degrees}{are latitude/longitude in degrees? Default: \code{TRUE}}
+
+\item{precision}{level of precision for the conversion. Default \code{5}}
+}
+\description{
+Convert latitude/longitude to MGRS string
+}
+\examples{
+latlng_to_mgrs(42, -93)
+}

man/mgrs.Rd

@@ -0,0 +1,13 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mgrs-package.R
+\docType{package}
+\name{mgrs}
+\alias{mgrs}
+\alias{mgrs-package}
+\title{Convert 'MGRS' Coordinates To and From Other Coordiante Systems}
+\description{
+Convert 'MGRS' Coordinates To and From Other Coordiante Systems
+}
+\author{
+Bob Rudis (bob@rud.is)
+}

man/mgrs_to_latlng.Rd

@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{mgrs_to_latlng}
+\alias{mgrs_to_latlng}
+\title{Convert an MGRS string to latitude/longitude}
+\usage{
+mgrs_to_latlng(MGRS, degrees = TRUE)
+}
+\arguments{
+\item{MGRS}{an MGRS string}
+
+\item{degrees}{convert to degrees? Default: \code{TRUE}}
+}
+\description{
+Convert an MGRS string to latitude/longitude
+}
+\examples{
+mgrs_to_latlng("15TWG0000049776")
+}

mgrs.Rproj

@@ -0,0 +1,21 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+StripTrailingWhitespace: Yes
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source
+PackageBuildArgs: --resave-data
+PackageRoxygenize: rd,collate,namespace

src/.gitignore

@@ -0,0 +1,3 @@
+*.o
+*.so
+*.dll

src/RcppExports.cpp

@@ -0,0 +1,33 @@
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#include <Rcpp.h>
+
+using namespace Rcpp;
+
+// mgrs_to_latlng
+NumericVector mgrs_to_latlng(std::string MGRS, bool degrees);
+RcppExport SEXP mgrs_mgrs_to_latlng(SEXP MGRSSEXP, SEXP degreesSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< std::string >::type MGRS(MGRSSEXP);
+    Rcpp::traits::input_parameter< bool >::type degrees(degreesSEXP);
+    rcpp_result_gen = Rcpp::wrap(mgrs_to_latlng(MGRS, degrees));
+    return rcpp_result_gen;
+END_RCPP
+}
+// latlng_to_mgrs
+std::string latlng_to_mgrs(double latitude, double longitude, bool degrees, int precision);
+RcppExport SEXP mgrs_latlng_to_mgrs(SEXP latitudeSEXP, SEXP longitudeSEXP, SEXP degreesSEXP, SEXP precisionSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< double >::type latitude(latitudeSEXP);
+    Rcpp::traits::input_parameter< double >::type longitude(longitudeSEXP);
+    Rcpp::traits::input_parameter< bool >::type degrees(degreesSEXP);
+    Rcpp::traits::input_parameter< int >::type precision(precisionSEXP);
+    rcpp_result_gen = Rcpp::wrap(latlng_to_mgrs(latitude, longitude, degrees, precision));
+    return rcpp_result_gen;
+END_RCPP
+}

src/main.cpp

@@ -0,0 +1,57 @@
+#include <Rcpp.h>
+
+#include "mgrs.h"
+
+using namespace Rcpp;
+
+//' Convert an MGRS string to latitude/longitude
+//'
+//' @md
+//' @param MGRS an MGRS string
+//' @param degrees convert to degrees? Default: `TRUE`
+//' @export
+//' @examples
+//' mgrs_to_latlng("15TWG0000049776")
+// [[Rcpp::export]]
+NumericVector mgrs_to_latlng(std::string MGRS, bool degrees = true) {
+
+  double lat, lng;
+  long ret;
+
+  ret = Convert_MGRS_To_Geodetic((char *)MGRS.c_str(), &lat, &lng);
+
+  NumericVector coords = NumericVector::create(
+    _["lat"] = degrees ? lat * 180.0/PI : lat,
+    _["lng"] = degrees ? lng * 180.0/PI : lng
+  );
+
+  return(coords);
+
+}
+
+//' Convert latitude/longitude to MGRS string
+//'
+//' @md
+//' @param latitude,longitude coordinates
+//' @param degrees are latitude/longitude in degrees? Default: `TRUE`
+//' @param precision level of precision for the conversion. Default `5`
+//' @export
+//' @examples
+//' latlng_to_mgrs(42, -93)
+// [[Rcpp::export]]
+std::string latlng_to_mgrs(double latitude, double longitude,
+                           bool degrees = true, int precision = 5) {
+
+  if (degrees) {
+    latitude *= PI / 180.0;
+    longitude *= PI / 180.0;
+  }
+
+  char buf[80];
+  int ret;
+
+  ret = Convert_Geodetic_To_MGRS(latitude, longitude, precision, buf);
+
+  return(std::string(buf));
+
+}

src/mgrs.h

@@ -0,0 +1,37 @@
+#ifndef MGRS_H
+#define MGRS_H
+
+#define MGRS_NO_ERROR                0x0000
+#define MGRS_LAT_ERROR               0x0001
+#define MGRS_LON_ERROR               0x0002
+#define MGRS_STRING_ERROR            0x0004
+#define MGRS_PRECISION_ERROR         0x0008
+#define MGRS_A_ERROR                 0x0010
+#define MGRS_INV_F_ERROR             0x0020
+#define MGRS_EASTING_ERROR           0x0040
+#define MGRS_NORTHING_ERROR          0x0080
+#define MGRS_ZONE_ERROR              0x0100
+#define MGRS_HEMISPHERE_ERROR        0x0200
+#define MGRS_LAT_WARNING             0x0400
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+long Set_MGRS_Parameters(double a, double f, char *Ellipsoid_Code);
+void Get_MGRS_Parameters(double *a, double *f, char *Ellipsoid_Code);
+long Convert_Geodetic_To_MGRS (double Latitude, double Longitude, long Precision, char *MGRS);
+long Convert_MGRS_To_Geodetic (char *MGRS, double *Latitude, double *Longitude);
+long Convert_UTM_To_MGRS (long Zone, char Hemisphere, double Easting,  double Northing,
+                          long Precision, char *MGRS);
+long Convert_MGRS_To_UTM (char *MGRS, long *Zone, char *Hemisphere, double *Easting,
+                          double *Northing);
+long Convert_UPS_To_MGRS (char Hemisphere, double Easting, double Northing, long Precision,
+                          char *MGRS);
+long Convert_MGRS_To_UPS (char *MGRS, char *Hemisphere, double *Easting, double *Northing);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* MGRS_H */

src/polarst.c

@@ -0,0 +1,523 @@
+/***************************************************************************/
+/* RSC IDENTIFIER: POLAR STEREOGRAPHIC 
+ *
+ *
+ * ABSTRACT
+ *
+ *    This component provides conversions between geodetic (latitude and
+ *    longitude) coordinates and Polar Stereographic (easting and northing) 
+ *    coordinates.
+ *
+ * ERROR HANDLING
+ *
+ *    This component checks parameters for valid values.  If an invalid 
+ *    value is found the error code is combined with the current error code 
+ *    using the bitwise or.  This combining allows multiple error codes to 
+ *    be returned. The possible error codes are:
+ *
+ *          POLAR_NO_ERROR           : No errors occurred in function
+ *          POLAR_LAT_ERROR          : Latitude outside of valid range
+ *                                      (-90 to 90 degrees)
+ *          POLAR_LON_ERROR          : Longitude outside of valid range
+ *                                      (-180 to 360 degrees) 
+ *          POLAR_ORIGIN_LAT_ERROR   : Latitude of true scale outside of valid
+ *                                      range (-90 to 90 degrees)
+ *          POLAR_ORIGIN_LON_ERROR   : Longitude down from pole outside of valid
+ *                                      range (-180 to 360 degrees)
+ *          POLAR_EASTING_ERROR      : Easting outside of valid range,
+ *                                      depending on ellipsoid and
+ *                                      projection parameters
+ *          POLAR_NORTHING_ERROR     : Northing outside of valid range,
+ *                                      depending on ellipsoid and
+ *                                      projection parameters
+ *          POLAR_RADIUS_ERROR       : Coordinates too far from pole,
+ *                                      depending on ellipsoid and
+ *                                      projection parameters
+ *          POLAR_A_ERROR            : Semi-major axis less than or equal to zero
+ *          POLAR_INV_F_ERROR        : Inverse flattening outside of valid range
+ *								  	                  (250 to 350)
+ *
+ *
+ * REUSE NOTES
+ *
+ *    POLAR STEREOGRAPHIC is intended for reuse by any application that  
+ *    performs a Polar Stereographic projection.
+ *
+ *
+ * REFERENCES
+ *
+ *    Further information on POLAR STEREOGRAPHIC can be found in the
+ *    Reuse Manual.
+ *
+ *
+ *    POLAR STEREOGRAPHIC originated from :
+ *                                U.S. Army Topographic Engineering Center
+ *                                Geospatial Information Division
+ *                                7701 Telegraph Road
+ *                                Alexandria, VA  22310-3864
+ *
+ *
+ * LICENSES
+ *
+ *    None apply to this component.
+ *
+ *
+ * RESTRICTIONS
+ *
+ *    POLAR STEREOGRAPHIC has no restrictions.
+ *
+ *
+ * ENVIRONMENT
+ *
+ *    POLAR STEREOGRAPHIC was tested and certified in the following
+ *    environments:
+ *
+ *    1. Solaris 2.5 with GCC, version 2.8.1
+ *    2. Window 95 with MS Visual C++, version 6
+ *
+ *
+ * MODIFICATIONS
+ *
+ *    Date              Description
+ *    ----              -----------
+ *    06-11-95          Original Code
+ *    03-01-97          Original Code
+ *
+ *
+ */
+
+
+/************************************************************************/
+/*
+ *                               INCLUDES
+ */
+
+#include <math.h>
+#include "polarst.h"
+
+/*
+ *    math.h     - Standard C math library
+ *    polarst.h  - Is for prototype error checking
+ */
+
+
+/************************************************************************/
+/*                               DEFINES
+ *
+ */
+
+
+#define PI           3.14159265358979323e0       /* PI     */
+#define PI_OVER_2    (PI / 2.0)           
+#define TWO_PI       (2.0 * PI)
+#define POLAR_POW(EsSin)     pow((1.0 - EsSin) / (1.0 + EsSin), es_OVER_2)
+
+/************************************************************************/
+/*                           GLOBAL DECLARATIONS
+ *
+ */
+
+const double PI_Over_4 = (PI / 4.0);
+
+/* Ellipsoid Parameters, default to WGS 84  */
+static double Polar_a = 6378137.0;                    /* Semi-major axis of ellipsoid in meters  */
+static double Polar_f = 1 / 298.257223563;            /* Flattening of ellipsoid  */
+static double es = 0.08181919084262188000;            /* Eccentricity of ellipsoid    */
+static double es_OVER_2 = .040909595421311;           /* es / 2.0 */
+static double Southern_Hemisphere = 0;                /* Flag variable */
+static double tc = 1.0;
+static double e4 = 1.0033565552493;
+static double Polar_a_mc = 6378137.0;                 /* Polar_a * mc */
+static double two_Polar_a = 12756274.0;               /* 2.0 * Polar_a */
+
+/* Polar Stereographic projection Parameters */
+static double Polar_Origin_Lat = ((PI * 90) / 180);   /* Latitude of origin in radians */
+static double Polar_Origin_Long = 0.0;                /* Longitude of origin in radians */
+static double Polar_False_Easting = 0.0;              /* False easting in meters */
+static double Polar_False_Northing = 0.0;             /* False northing in meters */
+
+/* Maximum variance for easting and northing values for WGS 84. */
+static double Polar_Delta_Easting = 12713601.0;
+static double Polar_Delta_Northing = 12713601.0;
+
+/* These state variables are for optimization purposes. The only function
+ * that should modify them is Set_Polar_Stereographic_Parameters.         
+ */
+
+
+/************************************************************************/
+/*                              FUNCTIONS
+ *
+ */
+
+
+long Set_Polar_Stereographic_Parameters (double a,
+                                         double f,
+                                         double Latitude_of_True_Scale,
+                                         double Longitude_Down_from_Pole,
+                                         double False_Easting,
+                                         double False_Northing)
+
+{  /* BEGIN Set_Polar_Stereographic_Parameters   */
+/*  
+ *  The function Set_Polar_Stereographic_Parameters receives the ellipsoid
+ *  parameters and Polar Stereograpic projection parameters as inputs, and
+ *  sets the corresponding state variables.  If any errors occur, error
+ *  code(s) are returned by the function, otherwise POLAR_NO_ERROR is returned.
+ *
+ *  a                : Semi-major axis of ellipsoid, in meters         (input)
+ *  f                : Flattening of ellipsoid					               (input)
+ *  Latitude_of_True_Scale  : Latitude of true scale, in radians       (input)
+ *  Longitude_Down_from_Pole : Longitude down from pole, in radians    (input)
+ *  False_Easting    : Easting (X) at center of projection, in meters  (input)
+ *  False_Northing   : Northing (Y) at center of projection, in meters (input)
+ */
+
+  double es2;
+  double slat, clat;
+  double essin;
+  double one_PLUS_es, one_MINUS_es;
+  double pow_es;
+  double temp, temp_northing;
+  double inv_f = 1 / f;
+  double mc;                    
+//  const double  epsilon = 1.0e-2;
+  long Error_Code = POLAR_NO_ERROR;
+
+  if (a <= 0.0)
+  { /* Semi-major axis must be greater than zero */
+    Error_Code |= POLAR_A_ERROR;
+  }
+  if ((inv_f < 250) || (inv_f > 350))
+  { /* Inverse flattening must be between 250 and 350 */
+    Error_Code |= POLAR_INV_F_ERROR;
+  }
+  if ((Latitude_of_True_Scale < -PI_OVER_2) || (Latitude_of_True_Scale > PI_OVER_2))
+  { /* Origin Latitude out of range */
+    Error_Code |= POLAR_ORIGIN_LAT_ERROR;
+  }
+  if ((Longitude_Down_from_Pole < -PI) || (Longitude_Down_from_Pole > TWO_PI))
+  { /* Origin Longitude out of range */
+    Error_Code |= POLAR_ORIGIN_LON_ERROR;
+  }
+
+  if (!Error_Code)
+  { /* no errors */
+
+    Polar_a = a;
+    two_Polar_a = 2.0 * Polar_a;
+    Polar_f = f;
+
+    if (Longitude_Down_from_Pole > PI)
+      Longitude_Down_from_Pole -= TWO_PI;
+    if (Latitude_of_True_Scale < 0)
+    {
+      Southern_Hemisphere = 1;
+      Polar_Origin_Lat = -Latitude_of_True_Scale;
+      Polar_Origin_Long = -Longitude_Down_from_Pole;
+    }
+    else
+    {
+      Southern_Hemisphere = 0;
+      Polar_Origin_Lat = Latitude_of_True_Scale;
+      Polar_Origin_Long = Longitude_Down_from_Pole;
+    }
+    Polar_False_Easting = False_Easting;
+    Polar_False_Northing = False_Northing;
+
+    es2 = 2 * Polar_f - Polar_f * Polar_f;
+    es = sqrt(es2);
+    es_OVER_2 = es / 2.0;
+
+    if (fabs(fabs(Polar_Origin_Lat) - PI_OVER_2) > 1.0e-10)
+    {
+      slat = sin(Polar_Origin_Lat);
+      essin = es * slat;
+      pow_es = POLAR_POW(essin);
+      clat = cos(Polar_Origin_Lat);
+      mc = clat / sqrt(1.0 - essin * essin);
+      Polar_a_mc = Polar_a * mc;
+      tc = tan(PI_Over_4 - Polar_Origin_Lat / 2.0) / pow_es;
+    }
+    else
+    {
+      one_PLUS_es = 1.0 + es;
+      one_MINUS_es = 1.0 - es;
+      e4 = sqrt(pow(one_PLUS_es, one_PLUS_es) * pow(one_MINUS_es, one_MINUS_es));
+    }
+
+    /* Calculate Radius */
+    Convert_Geodetic_To_Polar_Stereographic(0, Longitude_Down_from_Pole, 
+                                            &temp, &temp_northing);
+
+    Polar_Delta_Northing = temp_northing;
+    if(Polar_False_Northing)
+      Polar_Delta_Northing -= Polar_False_Northing;
+    if (Polar_Delta_Northing < 0)
+      Polar_Delta_Northing = -Polar_Delta_Northing;
+    Polar_Delta_Northing *= 1.01;
+
+    Polar_Delta_Easting = Polar_Delta_Northing;
+
+  /*  Polar_Delta_Easting = temp_northing;
+    if(Polar_False_Easting)
+      Polar_Delta_Easting -= Polar_False_Easting;
+    if (Polar_Delta_Easting < 0)
+      Polar_Delta_Easting = -Polar_Delta_Easting;
+    Polar_Delta_Easting *= 1.01;*/
+  }
+
+  return (Error_Code);
+} /* END OF Set_Polar_Stereographic_Parameters */
+
+
+
+void Get_Polar_Stereographic_Parameters (double *a,
+                                         double *f,
+                                         double *Latitude_of_True_Scale,
+                                         double *Longitude_Down_from_Pole,
+                                         double *False_Easting,
+                                         double *False_Northing)
+
+{ /* BEGIN Get_Polar_Stereographic_Parameters  */
+/*
+ * The function Get_Polar_Stereographic_Parameters returns the current
+ * ellipsoid parameters and Polar projection parameters.
+ *
+ *  a                : Semi-major axis of ellipsoid, in meters         (output)
+ *  f                : Flattening of ellipsoid					               (output)
+ *  Latitude_of_True_Scale  : Latitude of true scale, in radians       (output)
+ *  Longitude_Down_from_Pole : Longitude down from pole, in radians    (output)
+ *  False_Easting    : Easting (X) at center of projection, in meters  (output)
+ *  False_Northing   : Northing (Y) at center of projection, in meters (output)
+ */
+
+  *a = Polar_a;
+  *f = Polar_f;
+  *Latitude_of_True_Scale = Polar_Origin_Lat;
+  *Longitude_Down_from_Pole = Polar_Origin_Long;
+  *False_Easting = Polar_False_Easting;
+  *False_Northing = Polar_False_Northing;
+  return;
+} /* END OF Get_Polar_Stereographic_Parameters */
+
+
+long Convert_Geodetic_To_Polar_Stereographic (double Latitude,
+                                              double Longitude,
+                                              double *Easting,
+                                              double *Northing)
+
+{  /* BEGIN Convert_Geodetic_To_Polar_Stereographic */
+
+/*
+ * The function Convert_Geodetic_To_Polar_Stereographic converts geodetic
+ * coordinates (latitude and longitude) to Polar Stereographic coordinates
+ * (easting and northing), according to the current ellipsoid
+ * and Polar Stereographic projection parameters. If any errors occur, error
+ * code(s) are returned by the function, otherwise POLAR_NO_ERROR is returned.
+ *
+ *    Latitude   :  Latitude, in radians                      (input)
+ *    Longitude  :  Longitude, in radians                     (input)
+ *    Easting    :  Easting (X), in meters                    (output)
+ *    Northing   :  Northing (Y), in meters                   (output)
+ */
+
+  double dlam;
+  double slat;
+  double essin;
+  double t;
+  double rho;
+  double pow_es;
+  long Error_Code = POLAR_NO_ERROR;
+
+  if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
+  {   /* Latitude out of range */
+    Error_Code |= POLAR_LAT_ERROR;
+  }
+  if ((Latitude < 0) && (Southern_Hemisphere == 0))
+  {   /* Latitude and Origin Latitude in different hemispheres */
+    Error_Code |= POLAR_LAT_ERROR;
+  }
+  if ((Latitude > 0) && (Southern_Hemisphere == 1))
+  {   /* Latitude and Origin Latitude in different hemispheres */
+    Error_Code |= POLAR_LAT_ERROR;
+  }
+  if ((Longitude < -PI) || (Longitude > TWO_PI))
+  {  /* Longitude out of range */
+    Error_Code |= POLAR_LON_ERROR;
+  }
+
+
+  if (!Error_Code)
+  {  /* no errors */
+
+    if (fabs(fabs(Latitude) - PI_OVER_2) < 1.0e-10)
+    {
+      *Easting = Polar_False_Easting;
+      *Northing = Polar_False_Northing;
+    }
+    else
+    {
+      if (Southern_Hemisphere != 0)
+      {
+        Longitude *= -1.0;
+        Latitude *= -1.0;
+      }
+      dlam = Longitude - Polar_Origin_Long;
+      if (dlam > PI)
+      {
+        dlam -= TWO_PI;
+      }
+      if (dlam < -PI)
+      {
+        dlam += TWO_PI;
+      }
+      slat = sin(Latitude);
+      essin = es * slat;
+      pow_es = POLAR_POW(essin);
+      t = tan(PI_Over_4 - Latitude / 2.0) / pow_es;
+
+      if (fabs(fabs(Polar_Origin_Lat) - PI_OVER_2) > 1.0e-10)
+        rho = Polar_a_mc * t / tc;
+      else
+        rho = two_Polar_a * t / e4;
+
+
+      if (Southern_Hemisphere != 0)
+      {
+        *Easting = -(rho * sin(dlam) - Polar_False_Easting);
+     //   *Easting *= -1.0;
+        *Northing = rho * cos(dlam) + Polar_False_Northing;
+      }
+      else
+      {
+        *Easting = rho * sin(dlam) + Polar_False_Easting;
+        *Northing = -rho * cos(dlam) + Polar_False_Northing;
+      }
+
+    }
+  }
+  return (Error_Code);
+} /* END OF Convert_Geodetic_To_Polar_Stereographic */
+
+
+long Convert_Polar_Stereographic_To_Geodetic (double Easting,
+                                              double Northing,
+                                              double *Latitude,
+                                              double *Longitude)
+
+{ /*  BEGIN Convert_Polar_Stereographic_To_Geodetic  */
+/*
+ *  The function Convert_Polar_Stereographic_To_Geodetic converts Polar
+ *  Stereographic coordinates (easting and northing) to geodetic
+ *  coordinates (latitude and longitude) according to the current ellipsoid
+ *  and Polar Stereographic projection Parameters. If any errors occur, the
+ *  code(s) are returned by the function, otherwise POLAR_NO_ERROR
+ *  is returned.
+ *
+ *  Easting          : Easting (X), in meters                   (input)
+ *  Northing         : Northing (Y), in meters                  (input)
+ *  Latitude         : Latitude, in radians                     (output)
+ *  Longitude        : Longitude, in radians                    (output)
+ *
+ */
+
+  double dy = 0, dx = 0;
+  double rho = 0;
+  double t;
+  double PHI, sin_PHI;
+  double tempPHI = 0.0;
+  double essin;
+  double pow_es;
+  double delta_radius;
+  long Error_Code = POLAR_NO_ERROR;
+  double min_easting = Polar_False_Easting - Polar_Delta_Easting;
+  double max_easting = Polar_False_Easting + Polar_Delta_Easting;
+  double min_northing = Polar_False_Northing - Polar_Delta_Northing;
+  double max_northing = Polar_False_Northing + Polar_Delta_Northing;
+
+  if (Easting > max_easting || Easting < min_easting)
+  { /* Easting out of range */
+    Error_Code |= POLAR_EASTING_ERROR;
+  }
+  if (Northing > max_northing || Northing < min_northing)
+  { /* Northing out of range */
+    Error_Code |= POLAR_NORTHING_ERROR;
+  }
+
+  if (!Error_Code)
+  {
+    dy = Northing - Polar_False_Northing;
+    dx = Easting - Polar_False_Easting;
+
+    /* Radius of point with origin of false easting, false northing */
+    rho = sqrt(dx * dx + dy * dy);   
+    
+    delta_radius = sqrt(Polar_Delta_Easting * Polar_Delta_Easting + Polar_Delta_Northing * Polar_Delta_Northing);
+
+    if(rho > delta_radius)
+    { /* Point is outside of projection area */
+      Error_Code |= POLAR_RADIUS_ERROR;
+    }
+
+    if (!Error_Code)
+    { /* no errors */
+      if ((dy == 0.0) && (dx == 0.0))
+      {
+        *Latitude = PI_OVER_2;
+        *Longitude = Polar_Origin_Long;
+
+      }
+      else
+      {
+        if (Southern_Hemisphere != 0)
+        {
+          dy *= -1.0;
+          dx *= -1.0;
+        }
+
+        if (fabs(fabs(Polar_Origin_Lat) - PI_OVER_2) > 1.0e-10)
+          t = rho * tc / (Polar_a_mc);
+        else
+          t = rho * e4 / (two_Polar_a);
+        PHI = PI_OVER_2 - 2.0 * atan(t);
+        while (fabs(PHI - tempPHI) > 1.0e-10)
+        {
+          tempPHI = PHI;
+          sin_PHI = sin(PHI);
+          essin =  es * sin_PHI;
+          pow_es = POLAR_POW(essin);
+          PHI = PI_OVER_2 - 2.0 * atan(t * pow_es);
+        }
+        *Latitude = PHI;
+        *Longitude = Polar_Origin_Long + atan2(dx, -dy);
+
+        if (*Longitude > PI)
+          *Longitude -= TWO_PI;
+        else if (*Longitude < -PI)
+          *Longitude += TWO_PI;
+
+
+        if (*Latitude > PI_OVER_2)  /* force distorted values to 90, -90 degrees */
+          *Latitude = PI_OVER_2;
+        else if (*Latitude < -PI_OVER_2)
+          *Latitude = -PI_OVER_2;
+
+        if (*Longitude > PI)  /* force distorted values to 180, -180 degrees */
+          *Longitude = PI;
+        else if (*Longitude < -PI)
+          *Longitude = -PI;
+
+      }
+      if (Southern_Hemisphere != 0)
+      {
+        *Latitude *= -1.0;
+        *Longitude *= -1.0;
+      }
+    }
+  }
+  return (Error_Code);
+} /* END OF Convert_Polar_Stereographic_To_Geodetic */
+
+
+

src/polarst.h

@@ -0,0 +1,202 @@
+#ifndef POLARST_H
+  #define POLARST_H
+/***************************************************************************/
+/* RSC IDENTIFIER: POLAR STEREOGRAPHIC 
+ *
+ *
+ * ABSTRACT
+ *
+ *    This component provides conversions between geodetic (latitude and
+ *    longitude) coordinates and Polar Stereographic (easting and northing) 
+ *    coordinates.
+ *
+ * ERROR HANDLING
+ *
+ *    This component checks parameters for valid values.  If an invalid 
+ *    value is found the error code is combined with the current error code 
+ *    using the bitwise or.  This combining allows multiple error codes to 
+ *    be returned. The possible error codes are:
+ *
+ *          POLAR_NO_ERROR           : No errors occurred in function
+ *          POLAR_LAT_ERROR          : Latitude outside of valid range
+ *                                      (-90 to 90 degrees)
+ *          POLAR_LON_ERROR          : Longitude outside of valid range
+ *                                      (-180 to 360 degrees) 
+ *          POLAR_ORIGIN_LAT_ERROR   : Latitude of true scale outside of valid
+ *                                      range (-90 to 90 degrees)
+ *          POLAR_ORIGIN_LON_ERROR   : Longitude down from pole outside of valid
+ *                                      range (-180 to 360 degrees)
+ *          POLAR_EASTING_ERROR      : Easting outside of valid range,
+ *                                      depending on ellipsoid and
+ *                                      projection parameters
+ *          POLAR_NORTHING_ERROR     : Northing outside of valid range,
+ *                                      depending on ellipsoid and
+ *                                      projection parameters
+ *          POLAR_RADIUS_ERROR       : Coordinates too far from pole,
+ *                                      depending on ellipsoid and
+ *                                      projection parameters
+ *          POLAR_A_ERROR            : Semi-major axis less than or equal to zero
+ *          POLAR_INV_F_ERROR        : Inverse flattening outside of valid range
+ *								  	                  (250 to 350)
+ *
+ *
+ * REUSE NOTES
+ *
+ *    POLAR STEREOGRAPHIC is intended for reuse by any application that  
+ *    performs a Polar Stereographic projection.
+ *
+ *
+ * REFERENCES
+ *
+ *    Further information on POLAR STEREOGRAPHIC can be found in the
+ *    Reuse Manual.
+ *
+ *
+ *    POLAR STEREOGRAPHIC originated from :
+ *                                U.S. Army Topographic Engineering Center
+ *                                Geospatial Information Division
+ *                                7701 Telegraph Road
+ *                                Alexandria, VA  22310-3864
+ *
+ *
+ * LICENSES
+ *
+ *    None apply to this component.
+ *
+ *
+ * RESTRICTIONS
+ *
+ *    POLAR STEREOGRAPHIC has no restrictions.
+ *
+ *
+ * ENVIRONMENT
+ *
+ *    POLAR STEREOGRAPHIC was tested and certified in the following
+ *    environments:
+ *
+ *    1. Solaris 2.5 with GCC, version 2.8.1
+ *    2. Window 95 with MS Visual C++, version 6
+ *
+ *
+ * MODIFICATIONS
+ *
+ *    Date              Description
+ *    ----              -----------
+ *    06-11-95          Original Code
+ *    03-01-97          Original Code
+ *
+ *
+ */
+
+
+/**********************************************************************/
+/*
+ *                        DEFINES
+ */
+
+  #define POLAR_NO_ERROR                0x0000
+  #define POLAR_LAT_ERROR               0x0001
+  #define POLAR_LON_ERROR               0x0002
+  #define POLAR_ORIGIN_LAT_ERROR        0x0004
+  #define POLAR_ORIGIN_LON_ERROR        0x0008
+  #define POLAR_EASTING_ERROR			  0x0010
+  #define POLAR_NORTHING_ERROR		  0x0020
+  #define POLAR_A_ERROR                 0x0040
+  #define POLAR_INV_F_ERROR             0x0080
+  #define POLAR_RADIUS_ERROR            0x0100
+
+/**********************************************************************/
+/*
+ *                        FUNCTION PROTOTYPES
+ */
+
+/* ensure proper linkage to c++ programs */
+  #ifdef __cplusplus
+extern "C" {
+  #endif
+
+  long Set_Polar_Stereographic_Parameters (double a,
+                                           double f,
+                                           double Latitude_of_True_Scale,
+                                           double Longitude_Down_from_Pole,
+                                           double False_Easting,
+                                           double False_Northing);
+/*  
+ *  The function Set_Polar_Stereographic_Parameters receives the ellipsoid
+ *  parameters and Polar Stereograpic projection parameters as inputs, and
+ *  sets the corresponding state variables.  If any errors occur, error
+ *  code(s) are returned by the function, otherwise POLAR_NO_ERROR is returned.
+ *
+ *  a                : Semi-major axis of ellipsoid, in meters         (input)
+ *  f                : Flattening of ellipsoid                         (input)
+ *  Latitude_of_True_Scale  : Latitude of true scale, in radians       (input)
+ *  Longitude_Down_from_Pole : Longitude down from pole, in radians    (input)
+ *  False_Easting    : Easting (X) at center of projection, in meters  (input)
+ *  False_Northing   : Northing (Y) at center of projection, in meters (input)
+ */
+
+
+  void Get_Polar_Stereographic_Parameters (double *a,
+                                           double *f,
+                                           double *Latitude_of_True_Scale,
+                                           double *Longitude_Down_from_Pole,
+                                           double *False_Easting,
+                                           double *False_Northing);
+/*
+ * The function Get_Polar_Stereographic_Parameters returns the current
+ * ellipsoid parameters and Polar projection parameters.
+ *
+ *  a                : Semi-major axis of ellipsoid, in meters         (output)
+ *  f                : Flattening of ellipsoid                         (output)
+ *  Latitude_of_True_Scale  : Latitude of true scale, in radians       (output)
+ *  Longitude_Down_from_Pole : Longitude down from pole, in radians    (output)
+ *  False_Easting    : Easting (X) at center of projection, in meters  (output)
+ *  False_Northing   : Northing (Y) at center of projection, in meters (output)
+ */
+
+
+  long Convert_Geodetic_To_Polar_Stereographic (double Latitude,
+                                                double Longitude,
+                                                double *Easting,
+                                                double *Northing);
+/*
+ * The function Convert_Geodetic_To_Polar_Stereographic converts geodetic
+ * coordinates (latitude and longitude) to Polar Stereographic coordinates
+ * (easting and northing), according to the current ellipsoid
+ * and Polar Stereographic projection parameters. If any errors occur, error
+ * code(s) are returned by the function, otherwise POLAR_NO_ERROR is returned.
+ *
+ *    Latitude   :  Latitude, in radians                      (input)
+ *    Longitude  :  Longitude, in radians                     (input)
+ *    Easting    :  Easting (X), in meters                    (output)
+ *    Northing   :  Northing (Y), in meters                   (output)
+ */
+
+
+
+  long Convert_Polar_Stereographic_To_Geodetic (double Easting,
+                                                double Northing,
+                                                double *Latitude,
+                                                double *Longitude);
+
+/*
+ *  The function Convert_Polar_Stereographic_To_Geodetic converts Polar
+ *  Stereographic coordinates (easting and northing) to geodetic
+ *  coordinates (latitude and longitude) according to the current ellipsoid
+ *  and Polar Stereographic projection Parameters. If any errors occur, the
+ *  code(s) are returned by the function, otherwise POLAR_NO_ERROR
+ *  is returned.
+ *
+ *  Easting          : Easting (X), in meters                   (input)
+ *  Northing         : Northing (Y), in meters                  (input)
+ *  Latitude         : Latitude, in radians                     (output)
+ *  Longitude        : Longitude, in radians                    (output)
+ *
+ */
+
+  #ifdef __cplusplus
+}
+  #endif
+
+#endif  /* POLARST_H  */
+

src/tranmerc.c

@@ -0,0 +1,618 @@
+/***************************************************************************/
+/* RSC IDENTIFIER: TRANSVERSE MERCATOR
+ *
+ * ABSTRACT
+ *
+ *    This component provides conversions between Geodetic coordinates 
+ *    (latitude and longitude) and Transverse Mercator projection coordinates
+ *    (easting and northing).
+ *
+ * ERROR HANDLING
+ *
+ *    This component checks parameters for valid values.  If an invalid value
+ *    is found the error code is combined with the current error code using 
+ *    the bitwise or.  This combining allows multiple error codes to be
+ *    returned. The possible error codes are:
+ *
+ *       TRANMERC_NO_ERROR           : No errors occurred in function
+ *       TRANMERC_LAT_ERROR          : Latitude outside of valid range
+ *                                      (-90 to 90 degrees)
+ *       TRANMERC_LON_ERROR          : Longitude outside of valid range
+ *                                      (-180 to 360 degrees, and within
+ *                                        +/-90 of Central Meridian)
+ *       TRANMERC_EASTING_ERROR      : Easting outside of valid range
+ *                                      (depending on ellipsoid and
+ *                                       projection parameters)
+ *       TRANMERC_NORTHING_ERROR     : Northing outside of valid range
+ *                                      (depending on ellipsoid and
+ *                                       projection parameters)
+ *       TRANMERC_ORIGIN_LAT_ERROR   : Origin latitude outside of valid range
+ *                                      (-90 to 90 degrees)
+ *       TRANMERC_CENT_MER_ERROR     : Central meridian outside of valid range
+ *                                      (-180 to 360 degrees)
+ *       TRANMERC_A_ERROR            : Semi-major axis less than or equal to zero
+ *       TRANMERC_INV_F_ERROR        : Inverse flattening outside of valid range
+ *								  	                  (250 to 350)
+ *       TRANMERC_SCALE_FACTOR_ERROR : Scale factor outside of valid
+ *                                     range (0.3 to 3.0)
+ *		 TM_LON_WARNING              : Distortion will result if longitude is more
+ *                                       than 9 degrees from the Central Meridian
+ *
+ * REUSE NOTES
+ *
+ *    TRANSVERSE MERCATOR is intended for reuse by any application that 
+ *    performs a Transverse Mercator projection or its inverse.
+ *    
+ * REFERENCES
+ *
+ *    Further information on TRANSVERSE MERCATOR can be found in the 
+ *    Reuse Manual.
+ *
+ *    TRANSVERSE MERCATOR originated from :  
+ *                      U.S. Army Topographic Engineering Center
+ *                      Geospatial Information Division
+ *                      7701 Telegraph Road
+ *                      Alexandria, VA  22310-3864
+ *
+ * LICENSES
+ *
+ *    None apply to this component.
+ *
+ * RESTRICTIONS
+ *
+ *    TRANSVERSE MERCATOR has no restrictions.
+ *
+ * ENVIRONMENT
+ *
+ *    TRANSVERSE MERCATOR was tested and certified in the following 
+ *    environments:
+ *
+ *    1. Solaris 2.5 with GCC, version 2.8.1
+ *    2. Windows 95 with MS Visual C++, version 6
+ *
+ * MODIFICATIONS
+ *
+ *    Date              Description
+ *    ----              -----------
+ *    10-02-97          Original Code
+ *    03-02-97          Re-engineered Code
+ *
+ */
+
+
+/***************************************************************************/
+/*
+ *                               INCLUDES
+ */
+
+#include <math.h>
+#include "tranmerc.h"
+
+/*
+ *    math.h      - Standard C math library
+ *    tranmerc.h  - Is for prototype error checking
+ */
+
+
+/***************************************************************************/
+/*                               DEFINES 
+ *
+ */
+
+#define PI              3.14159265358979323e0   /* PI     */
+#define PI_OVER_2         (PI/2.0e0)            /* PI over 2 */
+#define MAX_LAT         ((PI * 89.99)/180.0)    /* 89.99 degrees in radians */
+#define MAX_DELTA_LONG  ((PI * 90)/180.0)       /* 90 degrees in radians */
+#define MIN_SCALE_FACTOR  0.3
+#define MAX_SCALE_FACTOR  3.0
+
+#define SPHTMD(Latitude) ((double) (TranMerc_ap * Latitude \
+      - TranMerc_bp * sin(2.e0 * Latitude) + TranMerc_cp * sin(4.e0 * Latitude) \
+      - TranMerc_dp * sin(6.e0 * Latitude) + TranMerc_ep * sin(8.e0 * Latitude) ) )
+
+#define SPHSN(Latitude) ((double) (TranMerc_a / sqrt( 1.e0 - TranMerc_es * \
+      pow(sin(Latitude), 2))))
+
+#define SPHSR(Latitude) ((double) (TranMerc_a * (1.e0 - TranMerc_es) / \
+    pow(DENOM(Latitude), 3)))
+
+#define DENOM(Latitude) ((double) (sqrt(1.e0 - TranMerc_es * pow(sin(Latitude),2))))
+
+
+/**************************************************************************/
+/*                               GLOBAL DECLARATIONS
+ *
+ */
+
+/* Ellipsoid Parameters, default to WGS 84  */
+static double TranMerc_a = 6378137.0;              /* Semi-major axis of ellipsoid in meters */
+static double TranMerc_f = 1 / 298.257223563;      /* Flattening of ellipsoid  */
+static double TranMerc_es = 0.0066943799901413800; /* Eccentricity (0.08181919084262188000) squared */
+static double TranMerc_ebs = 0.0067394967565869;   /* Second Eccentricity squared */
+
+/* Transverse_Mercator projection Parameters */
+static double TranMerc_Origin_Lat = 0.0;           /* Latitude of origin in radians */
+static double TranMerc_Origin_Long = 0.0;          /* Longitude of origin in radians */
+static double TranMerc_False_Northing = 0.0;       /* False northing in meters */
+static double TranMerc_False_Easting = 0.0;        /* False easting in meters */
+static double TranMerc_Scale_Factor = 1.0;         /* Scale factor  */
+
+/* Isometeric to geodetic latitude parameters, default to WGS 84 */
+static double TranMerc_ap = 6367449.1458008;
+static double TranMerc_bp = 16038.508696861;
+static double TranMerc_cp = 16.832613334334;
+static double TranMerc_dp = 0.021984404273757;
+static double TranMerc_ep = 3.1148371319283e-005;
+
+/* Maximum variance for easting and northing values for WGS 84. */
+static double TranMerc_Delta_Easting = 40000000.0;
+static double TranMerc_Delta_Northing = 40000000.0;
+
+/* These state variables are for optimization purposes. The only function
+ * that should modify them is Set_Tranverse_Mercator_Parameters.         */
+
+
+/************************************************************************/
+/*                              FUNCTIONS     
+ *
+ */
+
+
+long Set_Transverse_Mercator_Parameters(double a,
+                                        double f,
+                                        double Origin_Latitude,
+                                        double Central_Meridian,
+                                        double False_Easting,
+                                        double False_Northing,
+                                        double Scale_Factor)
+
+{ /* BEGIN Set_Tranverse_Mercator_Parameters */
+  /*
+   * The function Set_Tranverse_Mercator_Parameters receives the ellipsoid
+   * parameters and Tranverse Mercator projection parameters as inputs, and
+   * sets the corresponding state variables. If any errors occur, the error
+   * code(s) are returned by the function, otherwise TRANMERC_NO_ERROR is
+   * returned.
+   *
+   *    a                 : Semi-major axis of ellipsoid, in meters    (input)
+   *    f                 : Flattening of ellipsoid						         (input)
+   *    Origin_Latitude   : Latitude in radians at the origin of the   (input)
+   *                         projection
+   *    Central_Meridian  : Longitude in radians at the center of the  (input)
+   *                         projection
+   *    False_Easting     : Easting/X at the center of the projection  (input)
+   *    False_Northing    : Northing/Y at the center of the projection (input)
+   *    Scale_Factor      : Projection scale factor                    (input) 
+   */
+
+  double tn;        /* True Meridianal distance constant  */
+  double tn2;
+  double tn3;
+  double tn4;
+  double tn5;
+  double dummy_northing;
+  double TranMerc_b; /* Semi-minor axis of ellipsoid, in meters */
+  double inv_f = 1 / f;
+  long Error_Code = TRANMERC_NO_ERROR;
+
+  if (a <= 0.0)
+  { /* Semi-major axis must be greater than zero */
+    Error_Code |= TRANMERC_A_ERROR;
+  }
+  if ((inv_f < 250) || (inv_f > 350))
+  { /* Inverse flattening must be between 250 and 350 */
+    Error_Code |= TRANMERC_INV_F_ERROR;
+  }
+  if ((Origin_Latitude < -PI_OVER_2) || (Origin_Latitude > PI_OVER_2))
+  { /* origin latitude out of range */
+    Error_Code |= TRANMERC_ORIGIN_LAT_ERROR;
+  }
+  if ((Central_Meridian < -PI) || (Central_Meridian > (2*PI)))
+  { /* origin longitude out of range */
+    Error_Code |= TRANMERC_CENT_MER_ERROR;
+  }
+  if ((Scale_Factor < MIN_SCALE_FACTOR) || (Scale_Factor > MAX_SCALE_FACTOR))
+  {
+    Error_Code |= TRANMERC_SCALE_FACTOR_ERROR;
+  }
+  if (!Error_Code)
+  { /* no errors */
+    TranMerc_a = a;
+    TranMerc_f = f;
+    TranMerc_Origin_Lat = Origin_Latitude;
+    if (Central_Meridian > PI)
+      Central_Meridian -= (2*PI);
+    TranMerc_Origin_Long = Central_Meridian;
+    TranMerc_False_Northing = False_Northing;
+    TranMerc_False_Easting = False_Easting; 
+    TranMerc_Scale_Factor = Scale_Factor;
+
+    /* Eccentricity Squared */
+    TranMerc_es = 2 * TranMerc_f - TranMerc_f * TranMerc_f;
+    /* Second Eccentricity Squared */
+    TranMerc_ebs = (1 / (1 - TranMerc_es)) - 1;
+
+    TranMerc_b = TranMerc_a * (1 - TranMerc_f);    
+    /*True meridianal constants  */
+    tn = (TranMerc_a - TranMerc_b) / (TranMerc_a + TranMerc_b);
+    tn2 = tn * tn;
+    tn3 = tn2 * tn;
+    tn4 = tn3 * tn;
+    tn5 = tn4 * tn;
+
+    TranMerc_ap = TranMerc_a * (1.e0 - tn + 5.e0 * (tn2 - tn3)/4.e0
+                                + 81.e0 * (tn4 - tn5)/64.e0 );
+    TranMerc_bp = 3.e0 * TranMerc_a * (tn - tn2 + 7.e0 * (tn3 - tn4)
+                                       /8.e0 + 55.e0 * tn5/64.e0 )/2.e0;
+    TranMerc_cp = 15.e0 * TranMerc_a * (tn2 - tn3 + 3.e0 * (tn4 - tn5 )/4.e0) /16.0;
+    TranMerc_dp = 35.e0 * TranMerc_a * (tn3 - tn4 + 11.e0 * tn5 / 16.e0) / 48.e0;
+    TranMerc_ep = 315.e0 * TranMerc_a * (tn4 - tn5) / 512.e0;
+    Convert_Geodetic_To_Transverse_Mercator(MAX_LAT,
+                                            MAX_DELTA_LONG + Central_Meridian,
+                                            &TranMerc_Delta_Easting,
+                                            &TranMerc_Delta_Northing);
+    Convert_Geodetic_To_Transverse_Mercator(0,
+                                            MAX_DELTA_LONG + Central_Meridian,
+                                            &TranMerc_Delta_Easting,
+                                            &dummy_northing);
+    TranMerc_Delta_Northing++;
+    TranMerc_Delta_Easting++;
+
+  } /* END OF if(!Error_Code) */
+  return (Error_Code);
+}  /* END of Set_Transverse_Mercator_Parameters  */
+
+
+void Get_Transverse_Mercator_Parameters(double *a,
+                                        double *f,
+                                        double *Origin_Latitude,
+                                        double *Central_Meridian,
+                                        double *False_Easting,
+                                        double *False_Northing,
+                                        double *Scale_Factor)
+
+{ /* BEGIN Get_Tranverse_Mercator_Parameters  */
+  /*
+   * The function Get_Transverse_Mercator_Parameters returns the current
+   * ellipsoid and Transverse Mercator projection parameters.
+   *
+   *    a                 : Semi-major axis of ellipsoid, in meters    (output)
+   *    f                 : Flattening of ellipsoid						         (output)
+   *    Origin_Latitude   : Latitude in radians at the origin of the   (output)
+   *                         projection
+   *    Central_Meridian  : Longitude in radians at the center of the  (output)
+   *                         projection
+   *    False_Easting     : Easting/X at the center of the projection  (output)
+   *    False_Northing    : Northing/Y at the center of the projection (output)
+   *    Scale_Factor      : Projection scale factor                    (output) 
+   */
+
+  *a = TranMerc_a;
+  *f = TranMerc_f;
+  *Origin_Latitude = TranMerc_Origin_Lat;
+  *Central_Meridian = TranMerc_Origin_Long;
+  *False_Easting = TranMerc_False_Easting;
+  *False_Northing = TranMerc_False_Northing;
+  *Scale_Factor = TranMerc_Scale_Factor;
+  return;
+} /* END OF Get_Tranverse_Mercator_Parameters */
+
+
+
+long Convert_Geodetic_To_Transverse_Mercator (double Latitude,
+                                              double Longitude,
+                                              double *Easting,
+                                              double *Northing)
+
+{      /* BEGIN Convert_Geodetic_To_Transverse_Mercator */
+
+  /*
+   * The function Convert_Geodetic_To_Transverse_Mercator converts geodetic
+   * (latitude and longitude) coordinates to Transverse Mercator projection
+   * (easting and northing) coordinates, according to the current ellipsoid
+   * and Transverse Mercator projection coordinates.  If any errors occur, the
+   * error code(s) are returned by the function, otherwise TRANMERC_NO_ERROR is
+   * returned.
+   *
+   *    Latitude      : Latitude in radians                         (input)
+   *    Longitude     : Longitude in radians                        (input)
+   *    Easting       : Easting/X in meters                         (output)
+   *    Northing      : Northing/Y in meters                        (output)
+   */
+
+  double c;       /* Cosine of latitude                          */
+  double c2;
+  double c3;
+  double c5;
+  double c7;
+  double dlam;    /* Delta longitude - Difference in Longitude       */
+  double eta;     /* constant - TranMerc_ebs *c *c                   */
+  double eta2;
+  double eta3;
+  double eta4;
+  double s;       /* Sine of latitude                        */
+  double sn;      /* Radius of curvature in the prime vertical       */
+  double t;       /* Tangent of latitude                             */
+  double tan2;
+  double tan3;
+  double tan4;
+  double tan5;
+  double tan6;
+  double t1;      /* Term in coordinate conversion formula - GP to Y */
+  double t2;      /* Term in coordinate conversion formula - GP to Y */
+  double t3;      /* Term in coordinate conversion formula - GP to Y */
+  double t4;      /* Term in coordinate conversion formula - GP to Y */
+  double t5;      /* Term in coordinate conversion formula - GP to Y */
+  double t6;      /* Term in coordinate conversion formula - GP to Y */
+  double t7;      /* Term in coordinate conversion formula - GP to Y */
+  double t8;      /* Term in coordinate conversion formula - GP to Y */
+  double t9;      /* Term in coordinate conversion formula - GP to Y */
+  double tmd;     /* True Meridional distance                        */
+  double tmdo;    /* True Meridional distance for latitude of origin */
+  long    Error_Code = TRANMERC_NO_ERROR;
+  double temp_Origin;
+  double temp_Long;
+
+  if ((Latitude < -MAX_LAT) || (Latitude > MAX_LAT))
+  {  /* Latitude out of range */
+    Error_Code|= TRANMERC_LAT_ERROR;
+  }
+  if (Longitude > PI)
+    Longitude -= (2 * PI);
+  if ((Longitude < (TranMerc_Origin_Long - MAX_DELTA_LONG))
+      || (Longitude > (TranMerc_Origin_Long + MAX_DELTA_LONG)))
+  {
+    if (Longitude < 0)
+      temp_Long = Longitude + 2 * PI;
+    else
+      temp_Long = Longitude;
+    if (TranMerc_Origin_Long < 0)
+      temp_Origin = TranMerc_Origin_Long + 2 * PI;
+    else
+      temp_Origin = TranMerc_Origin_Long;
+    if ((temp_Long < (temp_Origin - MAX_DELTA_LONG))
+        || (temp_Long > (temp_Origin + MAX_DELTA_LONG)))
+      Error_Code|= TRANMERC_LON_ERROR;
+  }
+  if (!Error_Code)
+  { /* no errors */
+
+    /* 
+     *  Delta Longitude
+     */
+    dlam = Longitude - TranMerc_Origin_Long;
+
+    if (fabs(dlam) > (9.0 * PI / 180))
+    { /* Distortion will result if Longitude is more than 9 degrees from the Central Meridian */
+      Error_Code |= TRANMERC_LON_WARNING;
+    }
+
+    if (dlam > PI)
+      dlam -= (2 * PI);
+    if (dlam < -PI)
+      dlam += (2 * PI);
+    if (fabs(dlam) < 2.e-10)
+      dlam = 0.0;
+
+    s = sin(Latitude);
+    c = cos(Latitude);
+    c2 = c * c;
+    c3 = c2 * c;
+    c5 = c3 * c2;
+    c7 = c5 * c2;
+    t = tan (Latitude);
+    tan2 = t * t;
+    tan3 = tan2 * t;
+    tan4 = tan3 * t;
+    tan5 = tan4 * t;
+    tan6 = tan5 * t;
+    eta = TranMerc_ebs * c2;
+    eta2 = eta * eta;
+    eta3 = eta2 * eta;
+    eta4 = eta3 * eta;
+
+    /* radius of curvature in prime vertical */
+    sn = SPHSN(Latitude);
+
+    /* True Meridianal Distances */
+    tmd = SPHTMD(Latitude);
+
+    /*  Origin  */
+    tmdo = SPHTMD (TranMerc_Origin_Lat);
+
+    /* northing */
+    t1 = (tmd - tmdo) * TranMerc_Scale_Factor;
+    t2 = sn * s * c * TranMerc_Scale_Factor/ 2.e0;
+    t3 = sn * s * c3 * TranMerc_Scale_Factor * (5.e0 - tan2 + 9.e0 * eta 
+                                                + 4.e0 * eta2) /24.e0; 
+
+    t4 = sn * s * c5 * TranMerc_Scale_Factor * (61.e0 - 58.e0 * tan2
+                                                + tan4 + 270.e0 * eta - 330.e0 * tan2 * eta + 445.e0 * eta2
+                                                + 324.e0 * eta3 -680.e0 * tan2 * eta2 + 88.e0 * eta4 
+                                                -600.e0 * tan2 * eta3 - 192.e0 * tan2 * eta4) / 720.e0;
+
+    t5 = sn * s * c7 * TranMerc_Scale_Factor * (1385.e0 - 3111.e0 * 
+                                                tan2 + 543.e0 * tan4 - tan6) / 40320.e0;
+
+    *Northing = TranMerc_False_Northing + t1 + pow(dlam,2.e0) * t2
+                + pow(dlam,4.e0) * t3 + pow(dlam,6.e0) * t4
+                + pow(dlam,8.e0) * t5; 
+
+    /* Easting */
+    t6 = sn * c * TranMerc_Scale_Factor;
+    t7 = sn * c3 * TranMerc_Scale_Factor * (1.e0 - tan2 + eta ) /6.e0;
+    t8 = sn * c5 * TranMerc_Scale_Factor * (5.e0 - 18.e0 * tan2 + tan4
+                                            + 14.e0 * eta - 58.e0 * tan2 * eta + 13.e0 * eta2 + 4.e0 * eta3 
+                                            - 64.e0 * tan2 * eta2 - 24.e0 * tan2 * eta3 )/ 120.e0;
+    t9 = sn * c7 * TranMerc_Scale_Factor * ( 61.e0 - 479.e0 * tan2
+                                             + 179.e0 * tan4 - tan6 ) /5040.e0;
+
+    *Easting = TranMerc_False_Easting + dlam * t6 + pow(dlam,3.e0) * t7 
+               + pow(dlam,5.e0) * t8 + pow(dlam,7.e0) * t9;
+  }
+  return (Error_Code);
+} /* END OF Convert_Geodetic_To_Transverse_Mercator */
+
+
+long Convert_Transverse_Mercator_To_Geodetic (
+                                             double Easting,
+                                             double Northing,
+                                             double *Latitude,
+                                             double *Longitude)
+{      /* BEGIN Convert_Transverse_Mercator_To_Geodetic */
+
+  /*
+   * The function Convert_Transverse_Mercator_To_Geodetic converts Transverse
+   * Mercator projection (easting and northing) coordinates to geodetic
+   * (latitude and longitude) coordinates, according to the current ellipsoid
+   * and Transverse Mercator projection parameters.  If any errors occur, the
+   * error code(s) are returned by the function, otherwise TRANMERC_NO_ERROR is
+   * returned.
+   *
+   *    Easting       : Easting/X in meters                         (input)
+   *    Northing      : Northing/Y in meters                        (input)
+   *    Latitude      : Latitude in radians                         (output)
+   *    Longitude     : Longitude in radians                        (output)
+   */
+
+  double c;       /* Cosine of latitude                          */
+  double de;      /* Delta easting - Difference in Easting (Easting-Fe)    */
+  double dlam;    /* Delta longitude - Difference in Longitude       */
+  double eta;     /* constant - TranMerc_ebs *c *c                   */
+  double eta2;
+  double eta3;
+  double eta4;
+  double ftphi;   /* Footpoint latitude                              */
+  int    i;       /* Loop iterator                   */
+  double s;       /* Sine of latitude                        */
+  double sn;      /* Radius of curvature in the prime vertical       */
+  double sr;      /* Radius of curvature in the meridian             */
+  double t;       /* Tangent of latitude                             */
+  double tan2;
+  double tan4;
+  double t10;     /* Term in coordinate conversion formula - GP to Y */
+  double t11;     /* Term in coordinate conversion formula - GP to Y */
+  double t12;     /* Term in coordinate conversion formula - GP to Y */
+  double t13;     /* Term in coordinate conversion formula - GP to Y */
+  double t14;     /* Term in coordinate conversion formula - GP to Y */
+  double t15;     /* Term in coordinate conversion formula - GP to Y */
+  double t16;     /* Term in coordinate conversion formula - GP to Y */
+  double t17;     /* Term in coordinate conversion formula - GP to Y */
+  double tmd;     /* True Meridional distance                        */
+  double tmdo;    /* True Meridional distance for latitude of origin */
+  long Error_Code = TRANMERC_NO_ERROR;
+
+  if ((Easting < (TranMerc_False_Easting - TranMerc_Delta_Easting))
+      ||(Easting > (TranMerc_False_Easting + TranMerc_Delta_Easting)))
+  { /* Easting out of range  */
+    Error_Code |= TRANMERC_EASTING_ERROR;
+  }
+  if ((Northing < (TranMerc_False_Northing - TranMerc_Delta_Northing))
+      || (Northing > (TranMerc_False_Northing + TranMerc_Delta_Northing)))
+  { /* Northing out of range */
+    Error_Code |= TRANMERC_NORTHING_ERROR;
+  }
+
+  if (!Error_Code)
+  {
+    /* True Meridional Distances for latitude of origin */
+    tmdo = SPHTMD(TranMerc_Origin_Lat);
+
+    /*  Origin  */
+    tmd = tmdo +  (Northing - TranMerc_False_Northing) / TranMerc_Scale_Factor; 
+
+    /* First Estimate */
+    sr = SPHSR(0.e0);
+    ftphi = tmd/sr;
+
+    for (i = 0; i < 5 ; i++)
+    {
+      t10 = SPHTMD (ftphi);
+      sr = SPHSR(ftphi);
+      ftphi = ftphi + (tmd - t10) / sr;
+    }
+
+    /* Radius of Curvature in the meridian */
+    sr = SPHSR(ftphi);
+
+    /* Radius of Curvature in the meridian */
+    sn = SPHSN(ftphi);
+
+    /* Sine Cosine terms */
+    s = sin(ftphi);
+    c = cos(ftphi);
+
+    /* Tangent Value  */
+    t = tan(ftphi);
+    tan2 = t * t;
+    tan4 = tan2 * tan2;
+    eta = TranMerc_ebs * pow(c,2);
+    eta2 = eta * eta;
+    eta3 = eta2 * eta;
+    eta4 = eta3 * eta;
+    de = Easting - TranMerc_False_Easting;
+    if (fabs(de) < 0.0001)
+      de = 0.0;
+
+    /* Latitude */
+    t10 = t / (2.e0 * sr * sn * pow(TranMerc_Scale_Factor, 2));
+    t11 = t * (5.e0  + 3.e0 * tan2 + eta - 4.e0 * pow(eta,2)
+               - 9.e0 * tan2 * eta) / (24.e0 * sr * pow(sn,3) 
+                                       * pow(TranMerc_Scale_Factor,4));
+    t12 = t * (61.e0 + 90.e0 * tan2 + 46.e0 * eta + 45.E0 * tan4
+               - 252.e0 * tan2 * eta  - 3.e0 * eta2 + 100.e0 
+               * eta3 - 66.e0 * tan2 * eta2 - 90.e0 * tan4
+               * eta + 88.e0 * eta4 + 225.e0 * tan4 * eta2
+               + 84.e0 * tan2* eta3 - 192.e0 * tan2 * eta4)
+          / ( 720.e0 * sr * pow(sn,5) * pow(TranMerc_Scale_Factor, 6) );
+    t13 = t * ( 1385.e0 + 3633.e0 * tan2 + 4095.e0 * tan4 + 1575.e0 
+                * pow(t,6))/ (40320.e0 * sr * pow(sn,7) * pow(TranMerc_Scale_Factor,8));
+    *Latitude = ftphi - pow(de,2) * t10 + pow(de,4) * t11 - pow(de,6) * t12 
+                + pow(de,8) * t13;
+
+    t14 = 1.e0 / (sn * c * TranMerc_Scale_Factor);
+
+    t15 = (1.e0 + 2.e0 * tan2 + eta) / (6.e0 * pow(sn,3) * c * 
+                                        pow(TranMerc_Scale_Factor,3));
+
+    t16 = (5.e0 + 6.e0 * eta + 28.e0 * tan2 - 3.e0 * eta2
+           + 8.e0 * tan2 * eta + 24.e0 * tan4 - 4.e0 
+           * eta3 + 4.e0 * tan2 * eta2 + 24.e0 
+           * tan2 * eta3) / (120.e0 * pow(sn,5) * c  
+                             * pow(TranMerc_Scale_Factor,5));
+
+    t17 = (61.e0 +  662.e0 * tan2 + 1320.e0 * tan4 + 720.e0 
+           * pow(t,6)) / (5040.e0 * pow(sn,7) * c 
+                          * pow(TranMerc_Scale_Factor,7));
+
+    /* Difference in Longitude */
+    dlam = de * t14 - pow(de,3) * t15 + pow(de,5) * t16 - pow(de,7) * t17;
+
+    /* Longitude */
+    (*Longitude) = TranMerc_Origin_Long + dlam;
+
+    if((fabs)(*Latitude) > (90.0 * PI / 180.0))
+      Error_Code |= TRANMERC_NORTHING_ERROR;
+
+    if((*Longitude) > (PI))
+    {
+      *Longitude -= (2 * PI);
+      if((fabs)(*Longitude) > PI)
+        Error_Code |= TRANMERC_EASTING_ERROR;
+    }
+    else if((*Longitude) < (-PI))
+    {
+      *Longitude += (2 * PI);
+      if((fabs)(*Longitude) > PI)
+        Error_Code |= TRANMERC_EASTING_ERROR;
+    }
+
+    if (fabs(dlam) > (9.0 * PI / 180) * cos(*Latitude))
+    { /* Distortion will result if Longitude is more than 9 degrees from the Central Meridian at the equator */
+      /* and decreases to 0 degrees at the poles */
+      /* As you move towards the poles, distortion will become more significant */
+      Error_Code |= TRANMERC_LON_WARNING;
+    }
+  }
+  return (Error_Code);
+} /* END OF Convert_Transverse_Mercator_To_Geodetic */

src/tranmerc.h

@@ -0,0 +1,209 @@
+#ifndef TRANMERC_H
+  #define TRANMERC_H
+
+/***************************************************************************/
+/* RSC IDENTIFIER: TRANSVERSE MERCATOR
+ *
+ * ABSTRACT
+ *
+ *    This component provides conversions between Geodetic coordinates 
+ *    (latitude and longitude) and Transverse Mercator projection coordinates
+ *    (easting and northing).
+ *
+ * ERROR HANDLING
+ *
+ *    This component checks parameters for valid values.  If an invalid value
+ *    is found the error code is combined with the current error code using 
+ *    the bitwise or.  This combining allows multiple error codes to be
+ *    returned. The possible error codes are:
+ *
+ *       TRANMERC_NO_ERROR           : No errors occurred in function
+ *       TRANMERC_LAT_ERROR          : Latitude outside of valid range
+ *                                      (-90 to 90 degrees)
+ *       TRANMERC_LON_ERROR          : Longitude outside of valid range
+ *                                      (-180 to 360 degrees, and within
+ *                                        +/-90 of Central Meridian)
+ *       TRANMERC_EASTING_ERROR      : Easting outside of valid range
+ *                                      (depending on ellipsoid and
+ *                                       projection parameters)
+ *       TRANMERC_NORTHING_ERROR     : Northing outside of valid range
+ *                                      (depending on ellipsoid and
+ *                                       projection parameters)
+ *       TRANMERC_ORIGIN_LAT_ERROR   : Origin latitude outside of valid range
+ *                                      (-90 to 90 degrees)
+ *       TRANMERC_CENT_MER_ERROR     : Central meridian outside of valid range
+ *                                      (-180 to 360 degrees)
+ *       TRANMERC_A_ERROR            : Semi-major axis less than or equal to zero
+ *       TRANMERC_INV_F_ERROR        : Inverse flattening outside of valid range
+ *								  	                  (250 to 350)
+ *       TRANMERC_SCALE_FACTOR_ERROR : Scale factor outside of valid
+ *                                     range (0.3 to 3.0)
+ *		 TM_LON_WARNING              : Distortion will result if longitude is more
+ *                                      than 9 degrees from the Central Meridian
+ *
+ * REUSE NOTES
+ *
+ *    TRANSVERSE MERCATOR is intended for reuse by any application that 
+ *    performs a Transverse Mercator projection or its inverse.
+ *    
+ * REFERENCES
+ *
+ *    Further information on TRANSVERSE MERCATOR can be found in the 
+ *    Reuse Manual.
+ *
+ *    TRANSVERSE MERCATOR originated from :  
+ *                      U.S. Army Topographic Engineering Center
+ *                      Geospatial Information Division
+ *                      7701 Telegraph Road
+ *                      Alexandria, VA  22310-3864
+ *
+ * LICENSES
+ *
+ *    None apply to this component.
+ *
+ * RESTRICTIONS
+ *
+ *    TRANSVERSE MERCATOR has no restrictions.
+ *
+ * ENVIRONMENT
+ *
+ *    TRANSVERSE MERCATOR was tested and certified in the following 
+ *    environments:
+ *
+ *    1. Solaris 2.5 with GCC, version 2.8.1
+ *    2. Windows 95 with MS Visual C++, version 6
+ *
+ * MODIFICATIONS
+ *
+ *    Date              Description
+ *    ----              -----------
+ *    10-02-97          Original Code
+ *    03-02-97          Re-engineered Code
+ *
+ */
+
+
+/***************************************************************************/
+/*
+ *                              DEFINES
+ */
+
+  #define TRANMERC_NO_ERROR           0x0000
+  #define TRANMERC_LAT_ERROR          0x0001
+  #define TRANMERC_LON_ERROR          0x0002
+  #define TRANMERC_EASTING_ERROR      0x0004
+  #define TRANMERC_NORTHING_ERROR     0x0008
+  #define TRANMERC_ORIGIN_LAT_ERROR   0x0010
+  #define TRANMERC_CENT_MER_ERROR     0x0020
+  #define TRANMERC_A_ERROR            0x0040
+  #define TRANMERC_INV_F_ERROR        0x0080
+  #define TRANMERC_SCALE_FACTOR_ERROR 0x0100
+  #define TRANMERC_LON_WARNING        0x0200
+
+
+/***************************************************************************/
+/*
+ *                              FUNCTION PROTOTYPES
+ *                                for TRANMERC.C
+ */
+
+/* ensure proper linkage to c++ programs */
+  #ifdef __cplusplus
+extern "C" {
+  #endif
+
+
+  long Set_Transverse_Mercator_Parameters(double a,      
+                                          double f,
+                                          double Origin_Latitude,
+                                          double Central_Meridian,
+                                          double False_Easting,
+                                          double False_Northing,
+                                          double Scale_Factor);
+/*
+ * The function Set_Tranverse_Mercator_Parameters receives the ellipsoid
+ * parameters and Tranverse Mercator projection parameters as inputs, and
+ * sets the corresponding state variables. If any errors occur, the error
+ * code(s) are returned by the function, otherwise TRANMERC_NO_ERROR is
+ * returned.
+ *
+ *    a                 : Semi-major axis of ellipsoid, in meters    (input)
+ *    f                 : Flattening of ellipsoid                    (input)
+ *    Origin_Latitude   : Latitude in radians at the origin of the   (input)
+ *                         projection
+ *    Central_Meridian  : Longitude in radians at the center of the  (input)
+ *                         projection
+ *    False_Easting     : Easting/X at the center of the projection  (input)
+ *    False_Northing    : Northing/Y at the center of the projection (input)
+ *    Scale_Factor      : Projection scale factor                    (input) 
+ */
+
+
+  void Get_Transverse_Mercator_Parameters(double *a,
+                                          double *f,
+                                          double *Origin_Latitude,
+                                          double *Central_Meridian,
+                                          double *False_Easting,
+                                          double *False_Northing,
+                                          double *Scale_Factor);
+/*
+ * The function Get_Transverse_Mercator_Parameters returns the current
+ * ellipsoid and Transverse Mercator projection parameters.
+ *
+ *    a                 : Semi-major axis of ellipsoid, in meters    (output)
+ *    f                 : Flattening of ellipsoid                    (output)
+ *    Origin_Latitude   : Latitude in radians at the origin of the   (output)
+ *                         projection
+ *    Central_Meridian  : Longitude in radians at the center of the  (output)
+ *                         projection
+ *    False_Easting     : Easting/X at the center of the projection  (output)
+ *    False_Northing    : Northing/Y at the center of the projection (output)
+ *    Scale_Factor      : Projection scale factor                    (output) 
+ */
+
+
+  long Convert_Geodetic_To_Transverse_Mercator (double Latitude,
+                                                double Longitude,
+                                                double *Easting,
+                                                double *Northing);
+
+/*
+ * The function Convert_Geodetic_To_Transverse_Mercator converts geodetic
+ * (latitude and longitude) coordinates to Transverse Mercator projection
+ * (easting and northing) coordinates, according to the current ellipsoid
+ * and Transverse Mercator projection coordinates.  If any errors occur, the
+ * error code(s) are returned by the function, otherwise TRANMERC_NO_ERROR is
+ * returned.
+ *
+ *    Latitude      : Latitude in radians                         (input)
+ *    Longitude     : Longitude in radians                        (input)
+ *    Easting       : Easting/X in meters                         (output)
+ *    Northing      : Northing/Y in meters                        (output)
+ */
+
+
+  long Convert_Transverse_Mercator_To_Geodetic (double Easting,
+                                                double Northing,
+                                                double *Latitude,
+                                                double *Longitude);
+
+/*
+ * The function Convert_Transverse_Mercator_To_Geodetic converts Transverse
+ * Mercator projection (easting and northing) coordinates to geodetic
+ * (latitude and longitude) coordinates, according to the current ellipsoid
+ * and Transverse Mercator projection parameters.  If any errors occur, the
+ * error code(s) are returned by the function, otherwise TRANMERC_NO_ERROR is
+ * returned.
+ *
+ *    Easting       : Easting/X in meters                         (input)
+ *    Northing      : Northing/Y in meters                        (input)
+ *    Latitude      : Latitude in radians                         (output)
+ *    Longitude     : Longitude in radians                        (output)
+ */
+
+
+  #ifdef __cplusplus
+}
+  #endif
+
+#endif /* TRANMERC_H */

src/ups.c

@@ -0,0 +1,302 @@
+/********************************************************************/
+/* RSC IDENTIFIER: UPS
+ *
+ *
+ * ABSTRACT
+ *
+ *    This component provides conversions between geodetic (latitude
+ *    and longitude) coordinates and Universal Polar Stereographic (UPS)
+ *    projection (hemisphere, easting, and northing) coordinates.
+ *
+ *
+ * ERROR HANDLING
+ *
+ *    This component checks parameters for valid values.  If an 
+ *    invalid value is found the error code is combined with the 
+ *    current error code using the bitwise or.  This combining allows  
+ *    multiple error codes to be returned. The possible error codes 
+ *    are:
+ *
+ *         UPS_NO_ERROR           : No errors occurred in function
+ *         UPS_LAT_ERROR          : Latitude outside of valid range
+ *                                   (North Pole: 83.5 to 90,
+ *                                    South Pole: -79.5 to -90)
+ *         UPS_LON_ERROR          : Longitude outside of valid range
+ *                                   (-180 to 360 degrees)
+ *         UPS_HEMISPHERE_ERROR   : Invalid hemisphere ('N' or 'S')
+ *         UPS_EASTING_ERROR      : Easting outside of valid range,
+ *                                   (0 to 4,000,000m)
+ *         UPS_NORTHING_ERROR     : Northing outside of valid range,
+ *                                   (0 to 4,000,000m)
+ *         UPS_A_ERROR            : Semi-major axis less than or equal to zero
+ *         UPS_INV_F_ERROR        : Inverse flattening outside of valid range
+ *								  	               (250 to 350)
+ *
+ *
+ * REUSE NOTES
+ *
+ *    UPS is intended for reuse by any application that performs a Universal
+ *    Polar Stereographic (UPS) projection.
+ *
+ *
+ * REFERENCES
+ *
+ *    Further information on UPS can be found in the Reuse Manual.
+ *
+ *    UPS originated from :  U.S. Army Topographic Engineering Center
+ *                           Geospatial Information Division
+ *                           7701 Telegraph Road
+ *                           Alexandria, VA  22310-3864
+ *
+ *
+ * LICENSES
+ *
+ *    None apply to this component.
+ *
+ *
+ * RESTRICTIONS
+ *
+ *    UPS has no restrictions.
+ *
+ *
+ * ENVIRONMENT
+ *
+ *    UPS was tested and certified in the following environments:
+ *
+ *    1. Solaris 2.5 with GCC version 2.8.1
+ *    2. Windows 95 with MS Visual C++ version 6
+ *
+ *
+ * MODIFICATIONS
+ *
+ *    Date              Description
+ *    ----              -----------
+ *    06-11-95          Original Code
+ *    03-01-97          Original Code
+ *
+ *
+ */
+
+
+/************************************************************************/
+/*
+ *                               INCLUDES
+ */
+
+#include <math.h>
+#include "polarst.h"
+#include "ups.h"
+/*
+ *    math.h     - Is needed to call the math functions.
+ *    polar.h    - Is used to convert polar stereographic coordinates
+ *    ups.h      - Defines the function prototypes for the ups module.
+ */
+
+
+/************************************************************************/
+/*                               GLOBAL DECLARATIONS
+ *
+ */
+
+#define PI       3.14159265358979323e0  /* PI     */
+#define PI_OVER    (PI/2.0e0)           /* PI over 2 */
+#define MAX_LAT    ((PI * 90)/180.0)    /* 90 degrees in radians */
+#define MAX_ORIGIN_LAT ((81.114528 * PI) / 180.0)
+#define MIN_NORTH_LAT (83.5*PI/180.0)
+#define MIN_SOUTH_LAT (-79.5*PI/180.0)
+#define MIN_EAST_NORTH 0
+#define MAX_EAST_NORTH 4000000
+
+/* Ellipsoid Parameters, default to WGS 84  */
+static double UPS_a = 6378137.0;          /* Semi-major axis of ellipsoid in meters   */
+static double UPS_f = 1 / 298.257223563;  /* Flattening of ellipsoid  */
+const double UPS_False_Easting = 2000000;
+const double UPS_False_Northing = 2000000;
+static double UPS_Origin_Latitude = MAX_ORIGIN_LAT;  /*set default = North Hemisphere */
+static double UPS_Origin_Longitude = 0.0;
+
+
+/************************************************************************/
+/*                              FUNCTIONS
+ *
+ */
+
+
+long Set_UPS_Parameters( double a,
+                         double f)
+{
+/*
+ * The function SET_UPS_PARAMETERS receives the ellipsoid parameters and sets
+ * the corresponding state variables. If any errors occur, the error code(s)
+ * are returned by the function, otherwise UPS_NO_ERROR is returned.
+ *
+ *   a     : Semi-major axis of ellipsoid in meters (input)
+ *   f     : Flattening of ellipsoid					      (input)
+ */
+
+  double inv_f = 1 / f;
+  long Error_Code = UPS_NO_ERROR;
+
+  if (a <= 0.0)
+  { /* Semi-major axis must be greater than zero */
+    Error_Code |= UPS_A_ERROR;
+  }
+  if ((inv_f < 250) || (inv_f > 350))
+  { /* Inverse flattening must be between 250 and 350 */
+    Error_Code |= UPS_INV_F_ERROR;
+  }
+
+  if (!Error_Code)
+  { /* no errors */
+    UPS_a = a;
+    UPS_f = f;
+  }
+  return (Error_Code);
+}  /* END of Set_UPS_Parameters  */
+
+
+void Get_UPS_Parameters( double *a,
+                         double *f)
+{
+/*
+ * The function Get_UPS_Parameters returns the current ellipsoid parameters.
+ *
+ *  a      : Semi-major axis of ellipsoid, in meters (output)
+ *  f      : Flattening of ellipsoid					       (output)
+ */
+
+  *a = UPS_a;
+  *f = UPS_f;
+  return;
+} /* END OF Get_UPS_Parameters */
+
+
+long Convert_Geodetic_To_UPS ( double Latitude,
+                               double Longitude,
+                               char   *Hemisphere,
+                               double *Easting,
+                               double *Northing)
+{
+/*
+ *  The function Convert_Geodetic_To_UPS converts geodetic (latitude and
+ *  longitude) coordinates to UPS (hemisphere, easting, and northing)
+ *  coordinates, according to the current ellipsoid parameters. If any 
+ *  errors occur, the error code(s) are returned by the function, 
+ *  otherwide UPS_NO_ERROR is returned.
+ *
+ *    Latitude      : Latitude in radians                       (input)
+ *    Longitude     : Longitude in radians                      (input)
+ *    Hemisphere    : Hemisphere either 'N' or 'S'              (output)
+ *    Easting       : Easting/X in meters                       (output)
+ *    Northing      : Northing/Y in meters                      (output)
+ */
+
+  double tempEasting, tempNorthing;
+  long Error_Code = UPS_NO_ERROR;
+
+  if ((Latitude < -MAX_LAT) || (Latitude > MAX_LAT))
+  {   /* latitude out of range */
+    Error_Code |= UPS_LAT_ERROR;
+  }
+  if ((Latitude < 0) && (Latitude > MIN_SOUTH_LAT))
+    Error_Code |= UPS_LAT_ERROR;
+  if ((Latitude >= 0) && (Latitude < MIN_NORTH_LAT))
+    Error_Code |= UPS_LAT_ERROR;
+  if ((Longitude < -PI) || (Longitude > (2 * PI)))
+  {  /* slam out of range */
+    Error_Code |= UPS_LON_ERROR;
+  }
+
+  if (!Error_Code)
+  {  /* no errors */
+    if (Latitude < 0)
+    {
+      UPS_Origin_Latitude = -MAX_ORIGIN_LAT; 
+      *Hemisphere = 'S';
+    }
+    else
+    {
+      UPS_Origin_Latitude = MAX_ORIGIN_LAT; 
+      *Hemisphere = 'N';
+    }
+
+
+    Set_Polar_Stereographic_Parameters( UPS_a,
+                                        UPS_f,
+                                        UPS_Origin_Latitude,
+                                        UPS_Origin_Longitude,
+                                        UPS_False_Easting,
+                                        UPS_False_Northing);
+
+    Convert_Geodetic_To_Polar_Stereographic(Latitude,
+                                            Longitude,
+                                            &tempEasting,
+                                            &tempNorthing);
+
+    *Easting = tempEasting;
+    *Northing = tempNorthing;
+  }  /*  END of if(!Error_Code)   */
+
+  return Error_Code;
+}  /* END OF Convert_Geodetic_To_UPS  */
+
+
+long Convert_UPS_To_Geodetic(char   Hemisphere,
+                             double Easting,
+                             double Northing,
+                             double *Latitude,
+                             double *Longitude)
+{
+/*
+ *  The function Convert_UPS_To_Geodetic converts UPS (hemisphere, easting, 
+ *  and northing) coordinates to geodetic (latitude and longitude) coordinates
+ *  according to the current ellipsoid parameters.  If any errors occur, the 
+ *  error code(s) are returned by the function, otherwise UPS_NO_ERROR is 
+ *  returned.
+ *
+ *    Hemisphere    : Hemisphere either 'N' or 'S'              (input)
+ *    Easting       : Easting/X in meters                       (input)
+ *    Northing      : Northing/Y in meters                      (input)
+ *    Latitude      : Latitude in radians                       (output)
+ *    Longitude     : Longitude in radians                      (output)
+ */
+
+  long Error_Code = UPS_NO_ERROR;
+
+  if ((Hemisphere != 'N') && (Hemisphere != 'S'))
+    Error_Code |= UPS_HEMISPHERE_ERROR;
+  if ((Easting < MIN_EAST_NORTH) || (Easting > MAX_EAST_NORTH))
+    Error_Code |= UPS_EASTING_ERROR;
+  if ((Northing < MIN_EAST_NORTH) || (Northing > MAX_EAST_NORTH))
+    Error_Code |= UPS_NORTHING_ERROR;
+
+  if (Hemisphere =='N')
+  {UPS_Origin_Latitude = MAX_ORIGIN_LAT;}
+  if (Hemisphere =='S')
+  {UPS_Origin_Latitude = -MAX_ORIGIN_LAT;}
+
+  if (!Error_Code)
+  {   /*  no errors   */
+    Set_Polar_Stereographic_Parameters( UPS_a,
+                                        UPS_f,
+                                        UPS_Origin_Latitude,
+                                        UPS_Origin_Longitude,
+                                        UPS_False_Easting,
+                                        UPS_False_Northing);
+
+
+
+    Convert_Polar_Stereographic_To_Geodetic( Easting,
+                                             Northing,
+                                             Latitude,
+                                             Longitude); 
+
+
+    if ((*Latitude < 0) && (*Latitude > MIN_SOUTH_LAT))
+      Error_Code |= UPS_LAT_ERROR;
+    if ((*Latitude >= 0) && (*Latitude < MIN_NORTH_LAT))
+      Error_Code |= UPS_LAT_ERROR;
+  }  /*  END OF if(!Error_Code) */
+  return (Error_Code);
+}  /*  END OF Convert_UPS_To_Geodetic  */ 
+

src/ups.h

@@ -0,0 +1,175 @@
+#ifndef UPS_H
+  #define UPS_H
+/********************************************************************/
+/* RSC IDENTIFIER: UPS
+ *
+ *
+ * ABSTRACT
+ *
+ *    This component provides conversions between geodetic (latitude
+ *    and longitude) coordinates and Universal Polar Stereographic (UPS)
+ *    projection (hemisphere, easting, and northing) coordinates.
+ *
+ *
+ * ERROR HANDLING
+ *
+ *    This component checks parameters for valid values.  If an 
+ *    invalid value is found the error code is combined with the 
+ *    current error code using the bitwise or.  This combining allows  
+ *    multiple error codes to be returned. The possible error codes 
+ *    are:
+ *
+ *         UPS_NO_ERROR           : No errors occurred in function
+ *         UPS_LAT_ERROR          : Latitude outside of valid range
+ *                                   (North Pole: 83.5 to 90,
+ *                                    South Pole: -79.5 to -90)
+ *         UPS_LON_ERROR          : Longitude outside of valid range
+ *                                   (-180 to 360 degrees)
+ *         UPS_HEMISPHERE_ERROR   : Invalid hemisphere ('N' or 'S')
+ *         UPS_EASTING_ERROR      : Easting outside of valid range,
+ *                                   (0 to 4,000,000m)
+ *         UPS_NORTHING_ERROR     : Northing outside of valid range,
+ *                                   (0 to 4,000,000m)
+ *         UPS_A_ERROR            : Semi-major axis less than or equal to zero
+ *         UPS_INV_F_ERROR        : Inverse flattening outside of valid range
+ *								  	               (250 to 350)
+ *
+ *
+ * REUSE NOTES
+ *
+ *    UPS is intended for reuse by any application that performs a Universal
+ *    Polar Stereographic (UPS) projection.
+ *
+ *
+ * REFERENCES
+ *
+ *    Further information on UPS can be found in the Reuse Manual.
+ *
+ *    UPS originated from :  U.S. Army Topographic Engineering Center
+ *                           Geospatial Information Division
+ *                           7701 Telegraph Road
+ *                           Alexandria, VA  22310-3864
+ *
+ *
+ * LICENSES
+ *
+ *    None apply to this component.
+ *
+ *
+ * RESTRICTIONS
+ *
+ *    UPS has no restrictions.
+ *
+ *
+ * ENVIRONMENT
+ *
+ *    UPS was tested and certified in the following environments:
+ *
+ *    1. Solaris 2.5 with GCC version 2.8.1
+ *    2. Windows 95 with MS Visual C++ version 6
+ *
+ *
+ * MODIFICATIONS
+ *
+ *    Date              Description
+ *    ----              -----------
+ *    06-11-95          Original Code
+ *    03-01-97          Original Code
+ *
+ *
+ */
+
+
+/**********************************************************************/
+/*
+ *                        DEFINES
+ */
+
+  #define UPS_NO_ERROR                0x0000
+  #define UPS_LAT_ERROR               0x0001
+  #define UPS_LON_ERROR               0x0002
+  #define UPS_HEMISPHERE_ERROR        0x0004
+  #define UPS_EASTING_ERROR           0x0008
+  #define UPS_NORTHING_ERROR          0x0010
+  #define UPS_A_ERROR                 0x0020
+  #define UPS_INV_F_ERROR             0x0040
+
+
+/**********************************************************************/
+/*
+ *                        FUNCTION PROTOTYPES
+ *                          for UPS.C
+ */
+
+/* ensure proper linkage to c++ programs */
+  #ifdef __cplusplus
+extern "C" {
+  #endif
+
+  long Set_UPS_Parameters( double a,
+                           double f);
+/*
+ * The function SET_UPS_PARAMETERS receives the ellipsoid parameters and sets
+ * the corresponding state variables. If any errors occur, the error code(s)
+ * are returned by the function, otherwise UPS_NO_ERROR is returned.
+ *
+ *   a     : Semi-major axis of ellipsoid in meters (input)
+ *   f     : Flattening of ellipsoid                (input)
+ */
+
+
+  void Get_UPS_Parameters( double *a,
+                           double *f);
+/*
+ * The function Get_UPS_Parameters returns the current ellipsoid parameters.
+ *
+ *  a      : Semi-major axis of ellipsoid, in meters (output)
+ *  f      : Flattening of ellipsoid                 (output)
+ */
+
+
+  long Convert_Geodetic_To_UPS ( double Latitude,
+                                 double Longitude,
+                                 char   *Hemisphere,
+                                 double *Easting,
+                                 double *Northing);
+/*
+ *  The function Convert_Geodetic_To_UPS converts geodetic (latitude and
+ *  longitude) coordinates to UPS (hemisphere, easting, and northing)
+ *  coordinates, according to the current ellipsoid parameters. If any 
+ *  errors occur, the error code(s) are returned by the function, 
+ *  otherwide UPS_NO_ERROR is returned.
+ *
+ *    Latitude      : Latitude in radians                       (input)
+ *    Longitude     : Longitude in radians                      (input)
+ *    Hemisphere    : Hemisphere either 'N' or 'S'              (output)
+ *    Easting       : Easting/X in meters                       (output)
+ *    Northing      : Northing/Y in meters                      (output)
+ */
+
+
+  long Convert_UPS_To_Geodetic(char   Hemisphere,
+                               double Easting,
+                               double Northing,
+                               double *Latitude,
+                               double *Longitude);
+
+/*
+ *  The function Convert_UPS_To_Geodetic converts UPS (hemisphere, easting, 
+ *  and northing) coordinates to geodetic (latitude and longitude) coordinates
+ *  according to the current ellipsoid parameters.  If any errors occur, the 
+ *  error code(s) are returned by the function, otherwise UPS_NO_ERROR is 
+ *  returned.
+ *
+ *    Hemisphere    : Hemisphere either 'N' or 'S'              (input)
+ *    Easting       : Easting/X in meters                       (input)
+ *    Northing      : Northing/Y in meters                      (input)
+ *    Latitude      : Latitude in radians                       (output)
+ *    Longitude     : Longitude in radians                      (output)
+ */
+
+  #ifdef __cplusplus
+}
+  #endif
+
+#endif  /* UPS_H  */

src/utm.c

@@ -0,0 +1,354 @@
+/***************************************************************************/
+/* RSC IDENTIFIER: UTM
+ *
+ * ABSTRACT
+ *
+ *    This component provides conversions between geodetic coordinates 
+ *    (latitude and longitudes) and Universal Transverse Mercator (UTM)
+ *    projection (zone, hemisphere, easting, and northing) coordinates.
+ *
+ * ERROR HANDLING
+ *
+ *    This component checks parameters for valid values.  If an invalid value
+ *    is found, the error code is combined with the current error code using 
+ *    the bitwise or.  This combining allows multiple error codes to be
+ *    returned. The possible error codes are:
+ *
+ *          UTM_NO_ERROR           : No errors occurred in function
+ *          UTM_LAT_ERROR          : Latitude outside of valid range
+ *                                    (-80.5 to 84.5 degrees)
+ *          UTM_LON_ERROR          : Longitude outside of valid range
+ *                                    (-180 to 360 degrees)
+ *          UTM_EASTING_ERROR      : Easting outside of valid range
+ *                                    (100,000 to 900,000 meters)
+ *          UTM_NORTHING_ERROR     : Northing outside of valid range
+ *                                    (0 to 10,000,000 meters)
+ *          UTM_ZONE_ERROR         : Zone outside of valid range (1 to 60)
+ *          UTM_HEMISPHERE_ERROR   : Invalid hemisphere ('N' or 'S')
+ *          UTM_ZONE_OVERRIDE_ERROR: Zone outside of valid range
+ *                                    (1 to 60) and within 1 of 'natural' zone
+ *          UTM_A_ERROR            : Semi-major axis less than or equal to zero
+ *          UTM_INV_F_ERROR        : Inverse flattening outside of valid range
+ *								  	                (250 to 350)
+ *
+ * REUSE NOTES
+ *
+ *    UTM is intended for reuse by any application that performs a Universal
+ *    Transverse Mercator (UTM) projection or its inverse.
+ *    
+ * REFERENCES
+ *
+ *    Further information on UTM can be found in the Reuse Manual.
+ *
+ *    UTM originated from :  U.S. Army Topographic Engineering Center
+ *                           Geospatial Information Division
+ *                           7701 Telegraph Road
+ *                           Alexandria, VA  22310-3864
+ *
+ * LICENSES
+ *
+ *    None apply to this component.
+ *
+ * RESTRICTIONS
+ *
+ *    UTM has no restrictions.
+ *
+ * ENVIRONMENT
+ *
+ *    UTM was tested and certified in the following environments:
+ *
+ *    1. Solaris 2.5 with GCC, version 2.8.1
+ *    2. MSDOS with MS Visual C++, version 6
+ *
+ * MODIFICATIONS
+ *
+ *    Date              Description
+ *    ----              -----------
+ *    10-02-97          Original Code
+ *
+ */
+
+
+/***************************************************************************/
+/*
+ *                              INCLUDES
+ */
+#include "tranmerc.h"
+#include "utm.h"
+/*
+ *    tranmerc.h    - Is used to convert transverse mercator coordinates
+ *    utm.h         - Defines the function prototypes for the utm module.
+ */
+
+
+/***************************************************************************/
+/*
+ *                              DEFINES
+ */
+
+#define PI           3.14159265358979323e0    /* PI                        */
+#define MIN_LAT      ( (-80.5 * PI) / 180.0 ) /* -80.5 degrees in radians    */
+#define MAX_LAT      ( (84.5 * PI) / 180.0 )  /* 84.5 degrees in radians     */
+#define MIN_EASTING  100000
+#define MAX_EASTING  900000
+#define MIN_NORTHING 0
+#define MAX_NORTHING 10000000
+
+/***************************************************************************/
+/*
+ *                              GLOBAL DECLARATIONS
+ */
+
+static double UTM_a = 6378137.0;         /* Semi-major axis of ellipsoid in meters  */
+static double UTM_f = 1 / 298.257223563; /* Flattening of ellipsoid                 */
+static long   UTM_Override = 0;          /* Zone override flag                      */
+
+
+/***************************************************************************/
+/*
+ *                                FUNCTIONS
+ *
+ */
+
+long Set_UTM_Parameters(double a,      
+                        double f,
+                        long   override)
+{
+/*
+ * The function Set_UTM_Parameters receives the ellipsoid parameters and
+ * UTM zone override parameter as inputs, and sets the corresponding state
+ * variables.  If any errors occur, the error code(s) are returned by the 
+ * function, otherwise UTM_NO_ERROR is returned.
+ *
+ *    a                 : Semi-major axis of ellipsoid, in meters       (input)
+ *    f                 : Flattening of ellipsoid						            (input)
+ *    override          : UTM override zone, zero indicates no override (input)
+ */
+
+  double inv_f = 1 / f;
+  long Error_Code = UTM_NO_ERROR;
+
+  if (a <= 0.0)
+  { /* Semi-major axis must be greater than zero */
+    Error_Code |= UTM_A_ERROR;
+  }
+  if ((inv_f < 250) || (inv_f > 350))
+  { /* Inverse flattening must be between 250 and 350 */
+    Error_Code |= UTM_INV_F_ERROR;
+  }
+  if ((override < 0) || (override > 60))
+  {
+    Error_Code |= UTM_ZONE_OVERRIDE_ERROR;
+  }
+  if (!Error_Code)
+  { /* no errors */
+    UTM_a = a;
+    UTM_f = f;
+    UTM_Override = override;
+  }
+  return (Error_Code);
+} /* END OF Set_UTM_Parameters */
+
+
+void Get_UTM_Parameters(double *a,
+                        double *f,
+                        long   *override)
+{
+/*
+ * The function Get_UTM_Parameters returns the current ellipsoid
+ * parameters and UTM zone override parameter.
+ *
+ *    a                 : Semi-major axis of ellipsoid, in meters       (output)
+ *    f                 : Flattening of ellipsoid						            (output)
+ *    override          : UTM override zone, zero indicates no override (output)
+ */
+
+  *a = UTM_a;
+  *f = UTM_f;
+  *override = UTM_Override;
+} /* END OF Get_UTM_Parameters */
+
+
+long Convert_Geodetic_To_UTM (double Latitude,
+                              double Longitude,
+                              long   *Zone,
+                              char   *Hemisphere,
+                              double *Easting,
+                              double *Northing)
+{ 
+/*
+ * The function Convert_Geodetic_To_UTM converts geodetic (latitude and
+ * longitude) coordinates to UTM projection (zone, hemisphere, easting and
+ * northing) coordinates according to the current ellipsoid and UTM zone
+ * override parameters.  If any errors occur, the error code(s) are returned
+ * by the function, otherwise UTM_NO_ERROR is returned.
+ *
+ *    Latitude          : Latitude in radians                 (input)
+ *    Longitude         : Longitude in radians                (input)
+ *    Zone              : UTM zone                            (output)
+ *    Hemisphere        : North or South hemisphere           (output)
+ *    Easting           : Easting (X) in meters               (output)
+ *    Northing          : Northing (Y) in meters              (output)
+ */
+
+  long Lat_Degrees;
+  long Long_Degrees;
+  long temp_zone;
+  long Error_Code = UTM_NO_ERROR;
+  double Origin_Latitude = 0;
+  double Central_Meridian = 0;
+  double False_Easting = 500000;
+  double False_Northing = 0;
+  double Scale = 0.9996;
+
+  if ((Latitude < MIN_LAT) || (Latitude > MAX_LAT))
+  { /* Latitude out of range */
+    Error_Code |= UTM_LAT_ERROR;
+  }
+  if ((Longitude < -PI) || (Longitude > (2*PI)))
+  { /* Longitude out of range */
+    Error_Code |= UTM_LON_ERROR;
+  }
+  if (!Error_Code)
+  { /* no errors */
+    if((Latitude > -1.0e-9) && (Latitude < 0))
+      Latitude = 0.0;
+    if (Longitude < 0)
+      Longitude += (2*PI) + 1.0e-10;
+
+    Lat_Degrees = (long)(Latitude * 180.0 / PI);
+    Long_Degrees = (long)(Longitude * 180.0 / PI);
+
+    if (Longitude < PI)
+      temp_zone = (long)(31 + ((Longitude * 180.0 / PI) / 6.0));
+    else
+      temp_zone = (long)(((Longitude * 180.0 / PI) / 6.0) - 29);
+
+    if (temp_zone > 60)
+      temp_zone = 1;
+    /* UTM special cases */
+    if ((Lat_Degrees > 55) && (Lat_Degrees < 64) && (Long_Degrees > -1)
+        && (Long_Degrees < 3))
+      temp_zone = 31;
+    if ((Lat_Degrees > 55) && (Lat_Degrees < 64) && (Long_Degrees > 2)
+        && (Long_Degrees < 12))
+      temp_zone = 32;
+    if ((Lat_Degrees > 71) && (Long_Degrees > -1) && (Long_Degrees < 9))
+      temp_zone = 31;
+    if ((Lat_Degrees > 71) && (Long_Degrees > 8) && (Long_Degrees < 21))
+      temp_zone = 33;
+    if ((Lat_Degrees > 71) && (Long_Degrees > 20) && (Long_Degrees < 33))
+      temp_zone = 35;
+    if ((Lat_Degrees > 71) && (Long_Degrees > 32) && (Long_Degrees < 42))
+      temp_zone = 37;
+
+    if (UTM_Override)
+    {
+      if ((temp_zone == 1) && (UTM_Override == 60))
+        temp_zone = UTM_Override;
+      else if ((temp_zone == 60) && (UTM_Override == 1))
+        temp_zone = UTM_Override;
+      else if ((Lat_Degrees > 71) && (Long_Degrees > -1) && (Long_Degrees < 42))
+      {
+        if (((temp_zone-2) <= UTM_Override) && (UTM_Override <= (temp_zone+2)))
+          temp_zone = UTM_Override;
+        else
+          Error_Code = UTM_ZONE_OVERRIDE_ERROR;
+      }
+      else if (((temp_zone-1) <= UTM_Override) && (UTM_Override <= (temp_zone+1)))
+        temp_zone = UTM_Override;
+      else
+        Error_Code = UTM_ZONE_OVERRIDE_ERROR;
+    }
+    if (!Error_Code)
+    {
+      if (temp_zone >= 31)
+        Central_Meridian = (6 * temp_zone - 183) * PI / 180.0;
+      else
+        Central_Meridian = (6 * temp_zone + 177) * PI / 180.0;
+      *Zone = temp_zone;
+      if (Latitude < 0)
+      {
+        False_Northing = 10000000;
+        *Hemisphere = 'S';
+      }
+      else
+        *Hemisphere = 'N';
+      Set_Transverse_Mercator_Parameters(UTM_a, UTM_f, Origin_Latitude,
+                                         Central_Meridian, False_Easting, False_Northing, Scale);
+      Convert_Geodetic_To_Transverse_Mercator(Latitude, Longitude, Easting,
+                                              Northing);
+      if ((*Easting < MIN_EASTING) || (*Easting > MAX_EASTING))
+        Error_Code = UTM_EASTING_ERROR;
+      if ((*Northing < MIN_NORTHING) || (*Northing > MAX_NORTHING))
+        Error_Code |= UTM_NORTHING_ERROR;
+    }
+  } /* END OF if (!Error_Code) */
+  return (Error_Code);
+} /* END OF Convert_Geodetic_To_UTM */
+
+
+long Convert_UTM_To_Geodetic(long   Zone,
+                             char   Hemisphere,
+                             double Easting,
+                             double Northing,
+                             double *Latitude,
+                             double *Longitude)
+{
+/*
+ * The function Convert_UTM_To_Geodetic converts UTM projection (zone, 
+ * hemisphere, easting and northing) coordinates to geodetic(latitude
+ * and  longitude) coordinates, according to the current ellipsoid
+ * parameters.  If any errors occur, the error code(s) are returned
+ * by the function, otherwise UTM_NO_ERROR is returned.
+ *
+ *    Zone              : UTM zone                               (input)
+ *    Hemisphere        : North or South hemisphere              (input)
+ *    Easting           : Easting (X) in meters                  (input)
+ *    Northing          : Northing (Y) in meters                 (input)
+ *    Latitude          : Latitude in radians                    (output)
+ *    Longitude         : Longitude in radians                   (output)
+ */
+  long Error_Code = UTM_NO_ERROR;
+  long tm_error_code = UTM_NO_ERROR;
+  double Origin_Latitude = 0;
+  double Central_Meridian = 0;
+  double False_Easting = 500000;
+  double False_Northing = 0;
+  double Scale = 0.9996;
+
+  if ((Zone < 1) || (Zone > 60))
+    Error_Code |= UTM_ZONE_ERROR;
+  if ((Hemisphere != 'S') && (Hemisphere != 'N'))
+    Error_Code |= UTM_HEMISPHERE_ERROR;
+  if ((Easting < MIN_EASTING) || (Easting > MAX_EASTING))
+    Error_Code |= UTM_EASTING_ERROR;
+  if ((Northing < MIN_NORTHING) || (Northing > MAX_NORTHING))
+    Error_Code |= UTM_NORTHING_ERROR;
+  if (!Error_Code)
+  { /* no errors */
+    if (Zone >= 31)
+      Central_Meridian = ((6 * Zone - 183) * PI / 180.0 /*+ 0.00000005*/);
+    else
+      Central_Meridian = ((6 * Zone + 177) * PI / 180.0 /*+ 0.00000005*/);
+    if (Hemisphere == 'S')
+      False_Northing = 10000000;
+    Set_Transverse_Mercator_Parameters(UTM_a, UTM_f, Origin_Latitude,
+                                       Central_Meridian, False_Easting, False_Northing, Scale);
+
+    tm_error_code = Convert_Transverse_Mercator_To_Geodetic(Easting, Northing, Latitude, Longitude);
+    if(tm_error_code)
+    {
+      if(tm_error_code & TRANMERC_EASTING_ERROR)
+        Error_Code |= UTM_EASTING_ERROR;
+      if(tm_error_code & TRANMERC_NORTHING_ERROR)
+        Error_Code |= UTM_NORTHING_ERROR;
+    }
+
+    if ((*Latitude < MIN_LAT) || (*Latitude > MAX_LAT))
+    { /* Latitude out of range */
+      Error_Code |= UTM_NORTHING_ERROR;
+    }
+  }
+  return (Error_Code);
+} /* END OF Convert_UTM_To_Geodetic */

src/utm.h

@@ -0,0 +1,178 @@
+#ifndef UTM_H
+  #define UTM_H
+
+/***************************************************************************/
+/* RSC IDENTIFIER: UTM
+ *
+ * ABSTRACT
+ *
+ *    This component provides conversions between geodetic coordinates 
+ *    (latitude and longitudes) and Universal Transverse Mercator (UTM)
+ *    projection (zone, hemisphere, easting, and northing) coordinates.
+ *
+ * ERROR HANDLING
+ *
+ *    This component checks parameters for valid values.  If an invalid value
+ *    is found, the error code is combined with the current error code using 
+ *    the bitwise or.  This combining allows multiple error codes to be
+ *    returned. The possible error codes are:
+ *
+ *          UTM_NO_ERROR           : No errors occurred in function
+ *          UTM_LAT_ERROR          : Latitude outside of valid range
+ *                                    (-80.5 to 84.5 degrees)
+ *          UTM_LON_ERROR          : Longitude outside of valid range
+ *                                    (-180 to 360 degrees)
+ *          UTM_EASTING_ERROR      : Easting outside of valid range
+ *                                    (100,000 to 900,000 meters)
+ *          UTM_NORTHING_ERROR     : Northing outside of valid range
+ *                                    (0 to 10,000,000 meters)
+ *          UTM_ZONE_ERROR         : Zone outside of valid range (1 to 60)
+ *          UTM_HEMISPHERE_ERROR   : Invalid hemisphere ('N' or 'S')
+ *          UTM_ZONE_OVERRIDE_ERROR: Zone outside of valid range
+ *                                    (1 to 60) and within 1 of 'natural' zone
+ *          UTM_A_ERROR            : Semi-major axis less than or equal to zero
+ *          UTM_INV_F_ERROR        : Inverse flattening outside of valid range
+ *								  	                (250 to 350)
+ *
+ * REUSE NOTES
+ *
+ *    UTM is intended for reuse by any application that performs a Universal
+ *    Transverse Mercator (UTM) projection or its inverse.
+ *    
+ * REFERENCES
+ *
+ *    Further information on UTM can be found in the Reuse Manual.
+ *
+ *    UTM originated from :  U.S. Army Topographic Engineering Center
+ *                           Geospatial Information Division
+ *                           7701 Telegraph Road
+ *                           Alexandria, VA  22310-3864
+ *
+ * LICENSES
+ *
+ *    None apply to this component.
+ *
+ * RESTRICTIONS
+ *
+ *    UTM has no restrictions.
+ *
+ * ENVIRONMENT
+ *
+ *    UTM was tested and certified in the following environments:
+ *
+ *    1. Solaris 2.5 with GCC, version 2.8.1
+ *    2. MSDOS with MS Visual C++, version 6
+ *
+ * MODIFICATIONS
+ *
+ *    Date              Description
+ *    ----              -----------
+ *    10-02-97          Original Code
+ *
+ */
+
+
+/***************************************************************************/
+/*
+ *                              DEFINES
+ */
+
+  #define UTM_NO_ERROR            0x0000
+  #define UTM_LAT_ERROR           0x0001
+  #define UTM_LON_ERROR           0x0002
+  #define UTM_EASTING_ERROR       0x0004
+  #define UTM_NORTHING_ERROR      0x0008
+  #define UTM_ZONE_ERROR          0x0010
+  #define UTM_HEMISPHERE_ERROR    0x0020
+  #define UTM_ZONE_OVERRIDE_ERROR 0x0040
+  #define UTM_A_ERROR             0x0080
+  #define UTM_INV_F_ERROR         0x0100
+
+
+/***************************************************************************/
+/*
+ *                              FUNCTION PROTOTYPES
+ *                                for UTM.C
+ */
+
+/* ensure proper linkage to c++ programs */
+  #ifdef __cplusplus
+extern "C" {
+  #endif
+
+  long Set_UTM_Parameters(double a,      
+                          double f,
+                          long   override);
+/*
+ * The function Set_UTM_Parameters receives the ellipsoid parameters and
+ * UTM zone override parameter as inputs, and sets the corresponding state
+ * variables.  If any errors occur, the error code(s) are returned by the 
+ * function, otherwise UTM_NO_ERROR is returned.
+ *
+ *    a                 : Semi-major axis of ellipsoid, in meters       (input)
+ *    f                 : Flattening of ellipsoid                       (input)
+ *    override          : UTM override zone, zero indicates no override (input)
+ */
+
+
+  void Get_UTM_Parameters(double *a,
+                          double *f,
+                          long   *override);
+/*
+ * The function Get_UTM_Parameters returns the current ellipsoid
+ * parameters and UTM zone override parameter.
+ *
+ *    a                 : Semi-major axis of ellipsoid, in meters       (output)
+ *    f                 : Flattening of ellipsoid                       (output)
+ *    override          : UTM override zone, zero indicates no override (output)
+ */
+
+
+  long Convert_Geodetic_To_UTM (double Latitude,
+                                double Longitude,
+                                long   *Zone,
+                                char   *Hemisphere,
+                                double *Easting,
+                                double *Northing); 
+/*
+ * The function Convert_Geodetic_To_UTM converts geodetic (latitude and
+ * longitude) coordinates to UTM projection (zone, hemisphere, easting and
+ * northing) coordinates according to the current ellipsoid and UTM zone
+ * override parameters.  If any errors occur, the error code(s) are returned
+ * by the function, otherwise UTM_NO_ERROR is returned.
+ *
+ *    Latitude          : Latitude in radians                 (input)
+ *    Longitude         : Longitude in radians                (input)
+ *    Zone              : UTM zone                            (output)
+ *    Hemisphere        : North or South hemisphere           (output)
+ *    Easting           : Easting (X) in meters               (output)
+ *    Northing          : Northing (Y) in meters              (output)
+ */
+
+
+  long Convert_UTM_To_Geodetic(long   Zone,
+                               char   Hemisphere,
+                               double Easting,
+                               double Northing,
+                               double *Latitude,
+                               double *Longitude);
+/*
+ * The function Convert_UTM_To_Geodetic converts UTM projection (zone, 
+ * hemisphere, easting and northing) coordinates to geodetic(latitude
+ * and  longitude) coordinates, according to the current ellipsoid
+ * parameters.  If any errors occur, the error code(s) are returned
+ * by the function, otherwise UTM_NO_ERROR is returned.
+ *
+ *    Zone              : UTM zone                               (input)
+ *    Hemisphere        : North or South hemisphere              (input)
+ *    Easting           : Easting (X) in meters                  (input)
+ *    Northing          : Northing (Y) in meters                 (input)
+ *    Latitude          : Latitude in radians                    (output)
+ *    Longitude         : Longitude in radians                   (output)
+ */
+
+  #ifdef __cplusplus
+}
+  #endif
+
+#endif /* UTM_H */

tests/test-all.R

@@ -0,0 +1,2 @@
+library(testthat)
+test_check("mgrs")

tests/testthat/test-mgrs.R

@@ -0,0 +1,6 @@
+context("basic functionality")
+test_that("we can do something", {
+
+  #expect_that(some_function(), is_a("data.frame"))
+
+})