高斯正算
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace ArcGISPublic.GaussHelper
{
class JWtoGauss
{
private double E = -9999;
private double N = -9999;
private double a3 = -9999;
public double _a3
{
set
{
if (a3 == -9999)
{
GEOToGuass(_a, _b, _FE, _FN, _E, _N, _CM, _Nψ0, k);
}
}
get { return a3; }
}
private double a2 = -9999;
public double _a2
{
set
{
if (a2 == -9999)
{
GEOToGuass(_a, _b, _FE, _FN, _E, _N, _CM, _Nψ0, k);
}
}
get { return a2; }
}
private double _a = 6378245.0;
public double longR
{
set { _a = value; }
get { return _a; }
}
private double _b = 6356863.01877304730;
public double shortR
{
set { _b = value; }
get { return _b; }
}
private double _FE = 500000.0;
private double _FN = 0.0;
private double _CM = 105 / (180 / 3.1415926);//注意順序括號在後！
public double jisuancenter
{
set { _CM = CenterLonNum6(value); }
}
private double _Nψ0 = 0.0;
private double k = 1.0;
public double X
{
get
{
if (E == -9999)
{
GEOToGuass(_a, _b, _FE, _FN, _E, _N, _CM, _Nψ0, k);
}
return E;
}
}
public double Y
{
get
{
if (N == -9999)
{
GEOToGuass(_a, _b, _FE, _FN, _E, _N, _CM, _Nψ0, k);
}
return N;
}
}
private double _E;
public double _EX
{
set { _E = value / 57.29578; }
}
private double _N;
public double _NY
{
set { _N = value / 57.29578; }
}
private double atanh(double z)
{
return 0.5 * Math.Log((1 + z) / (1 - z));
}
private double asinh(double z)
{
return Math.Log(z + Math.Sqrt(z * z + 1));
}
private double CenterLonNum6(double Eλ)
{
int a = (int)((Eλ * 180 / 3.1415926) / 6);
double center = ((a + 1) * 6 - 3);
return center / (180 / 3.1415926);
}
private double CenterLonNum3(double Eλ)
{
//if (Eλ < 1.5)
//{
// //jc=jd;
// center = Eλ;
//}
double n = Eλ - 1.5;
int center = ((int)(n / 3) + 1) * 3;
//int center=n*3;
return center * 1.0;
}
private void GEOToGuass(double a, double b, double FE, double FN, double Eλ, double Nψ, double Eλ0, double Nψ0, double k0)
{
double f = (a - b) / a;
double e = Math.Sqrt((a * a - b * b)) / a;
double n = f / (2 - f);
double B = (a / (1 + n)) * (1 + n * n / 4 + n * n * n * n / 64);
double h1 = n / 2 - ((2.0 / 3) * n * n) + ((5.0 / 16) * n * n * n) + ((41.0 / 180) * n * n * n * n);
double h2 = ((13.0 / 48) * n * n) - ((3.0 / 5) * n * n * n) + ((557.0 / 1440) * n * n * n * n);
double h3 = ((61.0 / 240) * n * n * n) - ((103.0 / 140) * n * n * n * n);
double h4 = (49561.0 / 161280) * n * n * n * n;
//double ψ=0.0;
//double Q0 = asinh(Math.Tan(ψ))-(e*atanh(e*Math.Sin(ψ)));
//double β0 = atanh(Math.Sinh(Q0));
//double ξq0 = Math.Asin(Math.Sin(β0));
double Q = asinh(Math.Tan(Nψ)) - (e * atanh(e * Math.Sin(Nψ)));
double β = Math.Atan(Math.Sinh(Q));
double η0 = atanh(Math.Cos(β) * Math.Sin(Eλ - Eλ0));
//double η = (E - FE) / (B * k0);
//double ζ = (N - FN) / (B * k0);
double ζ0 = Math.Asin(Math.Sin(β) * Math.Cosh(η0));
double ζ1 = h1 * Math.Sin(2 * ζ0) * Math.Cosh(2 * η0);
double ζ2 = h2 * Math.Sin(4 * ζ0) * Math.Cosh(4 * η0);
double ζ3 = h3 * Math.Sin(6 * ζ0) * Math.Cosh(6 * η0);
double ζ4 = h4 * Math.Sin(8 * ζ0) * Math.Cosh(8 * η0);
double ζ = ζ0 + ζ1 + ζ2 + ζ3 + ζ4;
double η1 = h1 * Math.Cos(2 * ζ0) * Math.Sinh(2 * η0);
double η2 = h2 * Math.Cos(4 * ζ0) * Math.Sinh(4 * η0);
double η3 = h3 * Math.Cos(6 * ζ0) * Math.Sinh(6 * η0);
double η4 = h4 * Math.Cos(8 * ζ0) * Math.Sinh(8 * η0);
double η = η0 + η1 + η2 + η3 + η4;
E = FE + k0 * B * η;
N = FN + k0 * B * ζ;
a3 = 4495668.7826 - N;
a2 = 542576.2338 - E;
}
}
}