高斯反算
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace ArcGISPublic.GaussHelper
{
class Gausstojw
{
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 ψ = -9999;
private double lambda = -9999;
public double Longitude
{
get
{
if (lambda == -9999)
{
GuassToGeo(_a, _b, 1 / _invf, _FE, _FN, _X, _Y, _CM, _k);
}
return lambda;
}
}
private double _a = 6378245.0;
private double _b = (298.300000000000010000 - 1) / 298.300000000000010000 * 6378245.0;
private double _invf = 298.300000000000010000;
private double _FE = 500000;
private double _FN = 0;
private double _CM = 105;
private double _k = 1;
public int Zone
{
set
{
_CM = value * 6 - 3;
}
}
public double SemiMajor
{
set
{
_a = value;
}
get
{
return _a;
}
}
public double Lantitude
{
get
{
if (lambda == -9999)
{
GuassToGeo(_a, _b, 1 / _invf, _FE, _FN, _X, _Y, _CM, _k);
}
return ψ;
}
}
private double _X;
private double _Y;
public double X
{
set
{
_X = value;
}
}
public double Y
{
set
{
_Y = value;
}
}
private void GuassToGeo(double a, double b, double f, double FE, double FN, double E, double N, double λ0, double k0)
{
//double f = (a - b) / a;
double n = f / (2 - f);
double e = Math.Sqrt((a * a - b * b)) / a;
double B = (a / (1 + n)) * (1 + n * n / 4 + n * n * n * n / 64);
double h1 = (n / 2) - (2.0 / 3) * n * n + (37.0 / 96) * n * n * n - (1.0 / 360) * n * n * n * n;
double h2 = (1.0 / 48) * n * n + (1.0 / 15) * n * n * n - (437.0 / 1440) * n * n * n * n;
double h3 = (17.0 / 480) * n * n * n - (37.0 / 840) * n * n * n * n;
double h4 = (4397.0 / 161280) * n * n * n * n;
double η = (E - FE) / (B * k0);
double ζ = (N - FN) / (B * k0);
double ζ1 = h1 * Math.Sin(2 * ζ) * Math.Cosh(2 * η);
double ζ2 = h2 * Math.Sin(4 * ζ) * Math.Cosh(4 * η);
double ζ3 = h3 * Math.Sin(6 * ζ) * Math.Cosh(6 * η);
double ζ4 = h4 * Math.Sin(8 * ζ) * Math.Cosh(8 * η);
double ζ0 = ζ - (ζ1 + ζ2 + ζ3 + ζ4);
double η1 = h1 * Math.Cos(2 * ζ) * Math.Sinh(2 * η);
double η2 = h2 * Math.Cos(4 * ζ) * Math.Sinh(4 * η);
double η3 = h3 * Math.Cos(6 * ζ) * Math.Sinh(6 * η);
double η4 = h4 * Math.Cos(8 * ζ) * Math.Sinh(8 * η);
double η0 = η - (η1 + η2 + η3 + η4);
double Β = Math.Asin(Math.Sin(ζ0) / Math.Cosh(η0));
double Q = asinh(Math.Tan(Β));
double Q1 = Q + (e * atanh(e * Math.Tanh(Q)));
for (int i = 1; i < 400; i++)
{
Q1 = Q + (e * atanh(e * Math.Tanh(Q1)));
}
ψ = Math.Atan(Math.Sinh(Q1)) * 180 / 3.1415926;
lambda = λ0 + Math.Asin(Math.Tanh(η0) / Math.Cos(Β)) * 180 / 3.1415926;
//System.Console.WriteLine("wd:{0},jd:{1}", ψ, λ);
}
}
}