using System;
namespace NZTM_NZGD_Converter
{
///
/// Conversion of the c project by Land Information New Zealand
///
public class Converter
{
struct TransverseMercatorProjection
{
public double Meridian { get; set; }
public double ScaleFactor { get; set; }
public double OriginLatitude { get; set; }
public double FalseEasting { get; set; }
public double FalseNorthing { get; set; }
public double UnitsToMeters { get; set; }
public double a, rf, f, e2, ep2;
public double om;
///
/// Constructs a new TM Projected point
///
///
///
///
///
///
///
///
///
public TransverseMercatorProjection(double Meridian, double ScaleFactor, double OriginLatitude, double FalseEasting, double FalseNorthing, double UnitsToMeters, double a, double rf)
{
this.Meridian = Meridian;
this.ScaleFactor = ScaleFactor;
this.OriginLatitude = OriginLatitude;
this.FalseEasting = FalseNorthing;
this.FalseNorthing = FalseNorthing;
this.UnitsToMeters = UnitsToMeters;
if (rf != 0.0) f = 1.0 / rf;
else f = 0.0;
this.a = a;
this.rf = rf;
e2 = 2.0 * f - f * f;
ep2 = e2 / (1.0 - e2);
om = GetMeridianArc(e2, a, OriginLatitude);
}
}
struct GeodeticProjection
{
public double Latitude { get; set; }
public double Longitude { get; set; }
}
///
/// Returns the length of the meridional arc (Helmert formula)
/// Method based on Redfearn's formulation as expressed in GDA technical
/// manual at http://www.anzlic.org.au/icsm/gdatm/index.html
///
/// Source e2
/// Source a
/// Source lt
/// Arc-length in meters
static double GetMeridianArc(double e2, double a, double lt)
{
double e4 = e2 * e2;
double e6 = e4 * e2;
double A0 = 1 - (e2 / 4.0) - (3.0 * e4 / 64.0) - (5.0 * e6 / 256.0);
double A2 = (3.0 / 8.0) * (e2 + e4 / 4.0 + 15.0 * e6 / 128.0);
double A4 = (15.0 / 256.0) * (e4 + 3.0 * e6 / 4.0);
double A6 = 35.0 * e6 / 3072.0;
return a * (A0 * lt - A2 * Math.Sin(2 * lt) + A4 * Math.Sin(4 * lt) - A6 * Math.Sin(6 * lt));
}
/*************************************************************************/
/* */
/* foot_point_lat */
/* */
/* Calculates the foot point latitude from the meridional arc */
/* Method based on Redfearn's formulation as expressed in GDA technical*/
/* manual at http://www.anzlic.org.au/icsm/gdatm/index.html */
/* */
/* Takes parameters */
/* tm definition (for scale factor) */
/* meridional arc (metres) */
/* */
/* Returns the foot point latitude (radians) */
/*************************************************************************/
static double GetFootPointLatitude(double f, double a, double m)
{
double n = f / (2.0 - f);
double n2 = n * n;
double n3 = n2 * n;
double n4 = n2 * n2;
double g = a * (1.0 - n) * (1.0 - n2) * (1 + 9.0 * n2 / 4.0 + 225.0 * n4 / 64.0);
double sig = m / g;
return sig + (3.0 * n / 2.0 - 27.0 * n3 / 32.0) * Math.Sin(2.0 * sig)
+ (21.0 * n2 / 16.0 - 55.0 * n4 / 32.0) * Math.Sin(4.0 * sig)
+ (151.0 * n3 / 96.0) * Math.Sin(6.0 * sig)
+ (1097.0 * n4 / 512.0) * Math.Sin(8.0 * sig);
}
/***************************************************************************/
/* */
/* tmgeod */
/* */
/* Routine to convert from Tranverse Mercator to latitude and longitude. */
/* Method based on Redfearn's formulation as expressed in GDA technical */
/* manual at http://www.anzlic.org.au/icsm/gdatm/index.html */
/* */
/* Takes parameters */
/* input easting (metres) */
/* input northing (metres) */
/* output latitude (radians) */
/* output longitude (radians) */
/* */
/***************************************************************************/
static GeodeticProjection TransversMercatorToGeodetic(TransverseMercatorProjection projection, double ce, double cn)
{
double cn1 = (cn - projection.FalseNorthing) * projection.UnitsToMeters / projection.ScaleFactor + projection.om;
double fphi = GetFootPointLatitude(projection.f, projection.a, cn1);
double slt = Math.Sin(fphi);
double clt = Math.Cos(fphi);
double eslt = (1.0 - projection.e2 * slt * slt);
double eta = projection.a / Math.Sqrt(eslt);
double rho = eta * (1.0 - projection.e2) / eslt;
double psi = eta / rho;
double E = (ce - projection.FalseEasting) * projection.UnitsToMeters;
double x = E / (eta * projection.ScaleFactor);
double x2 = x * x;
double t = slt / clt;
double t2 = t * t;
double t4 = t2 * t2;
double trm1 = 1.0 / 2.0;
double trm2 = ((-4.0 * psi
+ 9.0 * (1 - t2)) * psi
+ 12.0 * t2) / 24.0;
double trm3 = ((((8.0 * (11.0 - 24.0 * t2) * psi
- 12.0 * (21.0 - 71.0 * t2)) * psi
+ 15.0 * ((15.0 * t2 - 98.0) * t2 + 15)) * psi
+ 180.0 * ((-3.0 * t2 + 5.0) * t2)) * psi + 360.0 * t4) / 720.0;
double trm4 = (((1575.0 * t2 + 4095.0) * t2 + 3633.0) * t2 + 1385.0) / 40320.0;
GeodeticProjection geodeticProjection = new ();
geodeticProjection.Latitude = fphi + (t * x * E / (projection.ScaleFactor * rho)) * (((trm4 * x2 - trm3) * x2 + trm2) * x2 - trm1);
trm1 = 1.0;
trm2 = (psi + 2.0 * t2) / 6.0;
trm3 = (((-4.0 * (1.0 - 6.0 * t2) * psi + (9.0 - 68.0 * t2)) * psi + 72.0 * t2) * psi + 24.0 * t4) / 120.0;
trm4 = (((720.0 * t2 + 1320.0) * t2 + 662.0) * t2 + 61.0) / 5040.0;
geodeticProjection.Longitude = projection.Meridian - (x / clt) * (((trm4 * x2 - trm3) * x2 + trm2) * x2 - trm1);
return geodeticProjection;
}
static TransverseMercatorProjection GeodeticToTransverseMercator(GeodeticProjection geodeticProjection)
{
TransverseMercatorProjection projection = new();
double dlon = geodeticProjection.Longitude -
}
}
}