private void GaussProjCal(double longitude, double latitude, ref double xval, ref double yval)
{
int ProjNo = 0; int ZoneWide; 带宽
double longitude1, latitude1, longitude0, latitude0, X0, Y0;
double a, f, e2, ee, NN, T, C, A, M, iPI;
iPI = 0.0174532925199433; 3.1415926535898/180.0;
ZoneWide = 6; 6度带宽
a = 6378137;
f = 1 / 298.257223563;
// a=6378245.0; f=1.0/298.3; //54年北京坐标系参数
a=6378140.0; f=1/298.257; //80年西安坐标系参数
ProjNo = (int)(longitude / ZoneWide);
longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
longitude0 = longitude0 * iPI; //中央经线转换为弧度
latitude0 = 0;
longitude1 = longitude * iPI; //经度转换为弧度
latitude1 = latitude * iPI; //纬度转换为弧度
e2 = 2 * f - f * f;
ee = e2 * (1.0 - e2);
NN = a / Math.Sqrt(1.0 - e2 * Math.Sin(latitude1) * Math.Sin(latitude1));// NN卯酉圈曲率半径,测量学里面用N表示
T = Math.Tan(latitude1) * Math.Tan(latitude1);
C = ee * Math.Cos(latitude1) * Math.Cos(latitude1);
A = (longitude1 - longitude0) * Math.Cos(latitude1);
M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latitude1 - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2
* e2 / 1024) * Math.Sin(2 * latitude1)
+ (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.Sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072) * Math.Sin(6 * latitude1));// M为子午线弧长,测量学里用大X表示
xval = (NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120));
yval = (M + NN * Math.Tan(latitude1) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720));
// X0 = 1000000*(ProjNo+1)+500000;
X0 = 500000;
Y0 = 0;
xval = xval + X0 ;
yval = yval + Y0;
}
{
int ProjNo = 0; int ZoneWide; 带宽
double longitude1, latitude1, longitude0, latitude0, X0, Y0;
double a, f, e2, ee, NN, T, C, A, M, iPI;
iPI = 0.0174532925199433; 3.1415926535898/180.0;
ZoneWide = 6; 6度带宽
a = 6378137;
f = 1 / 298.257223563;
// a=6378245.0; f=1.0/298.3; //54年北京坐标系参数
a=6378140.0; f=1/298.257; //80年西安坐标系参数
ProjNo = (int)(longitude / ZoneWide);
longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
longitude0 = longitude0 * iPI; //中央经线转换为弧度
latitude0 = 0;
longitude1 = longitude * iPI; //经度转换为弧度
latitude1 = latitude * iPI; //纬度转换为弧度
e2 = 2 * f - f * f;
ee = e2 * (1.0 - e2);
NN = a / Math.Sqrt(1.0 - e2 * Math.Sin(latitude1) * Math.Sin(latitude1));// NN卯酉圈曲率半径,测量学里面用N表示
T = Math.Tan(latitude1) * Math.Tan(latitude1);
C = ee * Math.Cos(latitude1) * Math.Cos(latitude1);
A = (longitude1 - longitude0) * Math.Cos(latitude1);
M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latitude1 - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2
* e2 / 1024) * Math.Sin(2 * latitude1)
+ (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.Sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072) * Math.Sin(6 * latitude1));// M为子午线弧长,测量学里用大X表示
xval = (NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120));
yval = (M + NN * Math.Tan(latitude1) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720));
// X0 = 1000000*(ProjNo+1)+500000;
X0 = 500000;
Y0 = 0;
xval = xval + X0 ;
yval = yval + Y0;
}