高斯投影正反算

{
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;
/// <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)
/// <param name="a">参考椭球长半轴</param>
/// <param name="f">参考椭球扁率倒数</param>
/// <param name="beltWidth">投影分带的带宽</param>
/// <param name="assumedCoord">是否使用假定坐标</param>
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);
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);
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);
//换算成假定坐标，平移500km，前面加带号
if (assumedCoord) y += 500000 + beltNum * 1000000;
}
/// <summary>
///从大地坐标到平面坐标的高斯正算
///
///指定中央子午线，用于进行邻带换算，此时必不使用假定坐标
/// </summary>
}
/// <summary>
///从大地坐标到平面坐标的高斯正算
/// </summary>
/// <param name="B">大地纬度</param>
/// <param name="L">大地经度</param>
/// <param name="x">平面纵轴</param>
/// <param name="y">平面横轴</param>
}
/// <summary>
///平面坐标（自然坐标）到大地坐标的高斯反算
/// </summary>
/// <param name="x">平面纵轴</param>
/// <param name="y">平面横轴</param>
/// <param name="B">大地纬度</param>
/// <param name="l">经度差</param>
/// <param name="B">大地纬度</param>
/// <param name="L">大地经度</param>
/// <param name="CenterL">中央子午线</param>
/// <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);
}
/// <summary>
///平面坐标（自然坐标或假定坐标）到大地坐标的高斯反算
///
///默认使用六度带
/// </summary>
/// <param name="x">平面纵轴</param>
/// <param name="y">平面横轴</param>
/// <param name="B">大地纬度</param>
/// <param name="x">平面纵轴</param>
/// <param name="y">平面横轴</param>
/// <param name="a">参考椭球长半轴</param>
/// <param name="f">参考椭球扁率倒数</param>
/// <param name="beltWidth">投影分带的带宽</param>
/// <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)
/// <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="y">平面横轴</param>
/// <param name="a">参考椭球长半轴</param>
/// <param name="f">参考椭球扁率倒数</param>
/// <param name="assumedCoord">是否使用假定坐标</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;
{
beltNum = (int)Math.Ceiling(y / 1000000) - 1;
y -= 1000000 * beltNum + 500000;