高斯投影正反算

public static void xy_BL(double x, double y, out double B, out double L, double a,double f, int beltWidth)
{
//如果为假定坐标，转换为自然坐标
int beltNum = 0;
if (y > 1000000)
}
/// <summary>
///从大地坐标到平面坐标的高斯正算
///
///默认的是六度带投影
/// </summary>
/// <param name="B">大地纬度</param>
/// <param name="L">大地经度</param>
/// <param name="x">平面纵轴</param>
{
double ee = (2 * f - 1) / f / f; //第一偏心率的平方
double ee2 = ee / (1 - ee); //第二偏心率的平方
double rB, tB, m;
rB = B * Math.PI / 180;
tB = Math.Tan(rB);
m = Math.Cos(rB) * dL * Math.PI / 180;
public static void BL_xy(double B, double L, out double x, out double y, double a, double f, int beltWidth, bool assumedCoord)
{
int beltNum; //投影分带的带号
beltNum = (int)Math.Ceiling((L - (beltWidth == 3 ? 1.5 : 0)) / beltWidth);
/// <param name="B">大地纬度</param>
/// <param name="L">大地经度</param>
/// <param name="CenterL">中央子午线</param>
/// <param name="x">平面纵轴</param>
/// <param name="y">平面横轴</param>
//换算成假定坐标，平移500km，前面加带号
if (assumedCoord) y += 500000 + beltNum * 1000000;
}
/// <summary>
///从大地坐标到平面坐标的高斯正算
///
///指定中央子午线，用于进行邻带换算，此时必不使用假定坐标
/// </summary>
if (beltWidth == 3 && beltNum * 3 == L - 1.5) beltNum += 1;
L -= beltNum * beltWidth - (beltWidth == 6 ? 3 : 0);
Bl_xy(B, L, out x, out y, a, f, beltWidth);
x = MeridianLength(B, a, f) + N * tB * x;
y = N * (m + (1 - tB * tB + it2) * Math.Pow(m, 3) / 6 + (5 - 18 * tB * tB + Math.Pow(tB, 4) + 14 * it2 - 58 * tB * tB * it2) * Math.Pow(m, 5) / 120);
class Gauss
{
#region高斯投影正反算
/// <summary>
///从大地坐标到平面坐标的高斯正算
///
///默认的是使用假定坐标的六度带投影
/// </summary>
/// <param name="B">大地纬度</param>
/// <param name="L">大地经度</param>
cA = 1 + 3 * ee / 4 + 45 * ee * ee / 64 + 175 * Math.Pow(ee, 3) / 256 + 11025 * Math.Pow(ee, 4) / 16384;
cB = 3 * ee / 4 + 15 * ee * ee / 16 + 525 * Math.Pow(ee, 3) / 512 + 2205 * Math.Pow(ee, 4) / 2048;
{
BL_xy(B, L, out x, out y, a, f, 6, true);
}
/// <summary>
///从大地坐标到平面坐标的高斯正算
///
///默认的是使用假定坐标
/// </summary>
/// <param name="B">大地纬度</param>
/// <param name="L">大地经度</param>
/// <param name="L">大地经度</param>
/// <param name="a">参考椭球长半轴</param>
/// <param name="f">参考椭球扁率倒数</param>
public static void xy_BL(double x, double y, out double B, out double L, double a, double f)
}
/// <summary>
///平面坐标（自然坐标）到大地坐标的高斯反算
/// </summary>
/// <param name="x">平面纵轴</param>
/// <param name="y">平面横轴</param>
/// <param name="B">大地纬度</param>
/// <param name="l">经度差</param>
/// <param name="a">参考椭球长半轴</param>
/// <param name="f">参考椭球扁率倒数</param>
/// <param name="beltWidth">投影分带的带宽</param>
private static void xy_Bl(double x, double y, out double B, out double l, double a, double f, int beltWidth)
/// <param name="B">大地纬度</param>
/// <param name="L">大地经度</param>
/// <param name="a">参考椭球长半轴</param>
/// <param name="f">参考椭球扁率倒数</param>
/// <param name="beltWidth">投影分带的带宽</param>
}
/// <summary>
///从大地坐标到平面坐标的高斯正算
/// </summary>
/// <param name="B">大地纬度</param>
/// <param name="L">大地经度</param>
/// <param name="x">平面纵轴</param>
/// <param name="y">平面横轴</param>
double N = a / Math.Sqrt(1 - ee * Math.Sin(rB) * Math.Sin(rB));
double it2 = ee2 * Math.Pow(Math.Cos(rB), 2);
x = m * m / 2 + (5 - tB * tB + 9 * it2 + 4 * it2 * it2) * Math.Pow(m, 4) / 24 + (61 - 58 * tB * tB + Math.Pow(tB, 4)) * Math.Pow(m, 6) / 720;
/// <param name="a">参考椭球长半轴</param>
/// <param name="f">参考椭球扁率倒数</param>
/// <param name="beltWidth">投影分带的带宽</param>
public static void Bl_xy(double B, double dL, out double x, out double y, double a, double f, int beltWidth)
/// Hale Waihona Puke Baiduparam name="x">平面纵轴</param>
/// <param name="y">平面横轴</param>
/// <param name="a">参考椭球长半轴</param>
/// <param name="f">参考椭球扁率倒数</param>
/// <param name="beltWidth">投影分带的带宽</param>
}
/// <summary>
///平面坐标（自然坐标或假定坐标）到大地坐标的高斯反算
///
///默认使用六度带
/// </summary>
/// <param name="x">平面纵轴</param>
/// <param name="y">平面横轴</param>
/// <param name="B">大地纬度</param>
/// <param name="y">平面横轴</param>
/// <param name="a">参考椭球长半轴</param>
/// <param name="f">参考椭球扁率倒数</param>
/// <param name="assumedCoord">是否使用假定坐标</param>
/// <param name="a">参考椭球长半轴</param>
/// <param name="f">参考椭球扁率倒数</param>
/// <param name="beltWidth">投影分带的带宽</param>
/// <param name="assumedCoord">是否使用假定坐标</param>
{
if (y > 1000000)
{
throw new Exception("坐标类型错误，应使用自然坐标");
}
double ee = (2 * f - 1) / f / f; //第一偏心率的平方
double ee2 = ee / (1 - ee); //第二偏心率的平方
double cA, cB, cC, cD, cE;
{
xy_BL(x, y, out B, out L, a, f, 6);
}
/// <summary>
///平面坐标（自然坐标或假定坐标）到大地坐标的高斯反算
/// </summary>
/// <param name="x">平面纵轴</param>
/// <param name="y">平面横轴</param>
public static void BL_xy(double B, double L, out double x, out double y, double a, double f, bool assumedCoord)
{
BL_xy(B, L, out x, out y, a, f, 6, assumedCoord);
{
beltNum = (int)Math.Ceiling(y / 1000000) - 1;
y -= 1000000 * beltNum + 500000;
}
//求解纬度与经差
xy_Bl(x, y, out B, out L, a, f, beltWidth);
//求解经度
L += beltWidth * beltNum - ((beltWidth == 6) ? 3 : 0);
/// <param name="x">平面纵轴</param>
/// <param name="y">平面横轴</param>
/// <param name="a">参考椭球长半轴</param>
/// <param name="f">参考椭球扁率倒数</param>
public static void BL_xy(double B, double L, out double x, out double y, double a, double f)
public static void BL_xy(double B, double L, out double x, out double y, double a, double f, int beltWidth)
{
BL_xy(B, L, out x, out y, a, f, beltWidth, true);