Geo-Coordinates-MGRS-XS
view release on metacpan or search on metacpan
* 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;
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)
( run in 4.013 seconds using v1.01-cache-2.11-cpan-a49fcb8fa48 )