mirror of https://github.com/wb2osz/direwolf.git
619 lines
24 KiB
C
619 lines
24 KiB
C
/***************************************************************************/
|
|
/* 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 */
|