c语言坐标转换程序,高斯投影转换的C语言程序

高斯投影转换的C语言程序

(2008-11-17 09:47:55)

标签：

高斯投影

杂谈

#include

#include

const double pi = 3.1415926535;

void main()

{

// cout << "select a mode:"

<< endl;

cout << "选择大地基准面:"

<< endl

<

<< endl

<< "2 -西安1980"

<< endl

<< "3 -WGS84"

<< endl;

char i;

cin >> i;

double a,b;

switch(i) {

case '1':

a = 6378245; b = 6356863.0188;

break;

case '2':

a = 6378140; b = 6356755.2882;

break;

default:

a = 6378137; b = 6356752.3142;

}

double f = (a - b) / a;

double e = sqrt(2*f - f*f);

double ei = sqrt(e*e / (1-e*e));

double FE = 500000;

double FN = 0;

double k0 = 1;

double latitude;

double longtitude;

double longtitude0;

double latitude0 = 0;

while (1){

cout << "form 经纬度 to 高斯坐标"

<< endl;

cout << "请输入中央经度： " ;

cin >> longtitude0;

double field = (longtitude0 + 3) / 6;

cout <

;

cin >> latitude;

cout << "请输入经度坐标：" ;

cin >> longtitude;

cout << endl;

latitude = latitude * pi / 180;

longtitude = longtitude * pi / 180;

longtitude0 = longtitude0 * pi / 180;

double T = tan(latitude) * tan(latitude);

double C = e*e * cos(latitude)*cos(latitude) / (1 - e*e);

double A = (longtitude - longtitude0) * cos(latitude);

double v = a / sqrt((1 - e*e * sin(latitude)*sin(latitude)));

double M = a * ((1 - e*e/4 - 3*e*e*e*e/64 -

5*e*e*e*e*e*e/256)*latitude

- (3*e*e/8 + 3*e*e*e*e/32 +

45*e*e*e*e*e*e/1024)*sin(2*latitude)

+ (15*e*e*e*e/256 + 45*e*e*e*e*e*e/1024)*sin(4*latitude)

- 35*e*e*e*e*e*e/3072 * sin(6*latitude));

double M0 = a * ((1 - e*e/4 - 3*e*e*e*e/64 -

5*e*e*e*e*e*e/256)*latitude0

- (3*e*e/8 + 3*e*e*e*e/32 +

45*e*e*e*e*e*e/1024)*sin(2*latitude0)

+ (15*e*e*e*e/256 + 45*e*e*e*e*e*e/1024)*sin(4*latitude0)

- 35*e*e*e*e*e*e/3072 * sin(6*latitude0));

double Northing = FN + k0*(M - M0 + v*tan(latitude)*(A*A/2 +

(5-T+9*C+4*C*C)*A*A*A*A/24

+ (61-58*T+T*T+600*C-330*ei*ei)*A*A*A*A*A*A/720));

double Easting = FE + k0*v*(A + (1-T+C)*A*A*A/6 +

(5-18*T+T*T+72*C-58*ei*ei)*A*A*A*A*A/120) + field*1000000;

cout << "Northing = "

<< Northing

<< endl;

cout << "Easting = "

<< Easting

<< endl;

char key;

cout << "Quit? y/n";

cin >> key;

if (key == 'y')

break;

}

while(1){

cout << "from 高斯坐标 to 经纬度"

<< endl;

cout << "请输入中央经度： " ;

cin >> longtitude0;

double field = (longtitude0 + 3) / 6;

double north ;

double east ;

cout << "请输入north坐标： " ;

cin >> north;

cout << "请输入east坐标: " ;

cin >> east;

east = east - field * 1000000;

longtitude0 = longtitude0 * pi / 180;

double e1 = (1 - sqrt(1-e*e)) / (1 + sqrt(1 - e*e));

double M0 = a * ((1 - e*e/4 - 3*e*e*e*e/64 -

5*e*e*e*e*e*e/256)*latitude0

- (3*e*e/8 + 3*e*e*e*e/32 +

45*e*e*e*e*e*e/1024)*sin(2*latitude0)

+ (15*e*e*e*e/256 + 45*e*e*e*e*e*e/1024)*sin(4*latitude0)

- 35*e*e*e*e*e*e/3072 * sin(6*latitude0));

double M1 = M0 + (north - FN) * k0;

double u1 = M1 /

(a*(1-e*e/4-3*e*e*e*e/64-5*e*e*e*e*e*e/256));

double latitude1 = u1 + (3*e1/2 - 27*e1*e1*e1/32)*sin(2*u1) +

(21*e1*e1/16

- 55*e1*e1*e1*e1/32)*sin(4*u1) + 151*e1*e1*e1/96*sin(6*u1)

+ 1097*e1*e1*e1*e1/512*sin(8*u1);

double T1 = tan(latitude1) * tan(latitude1);

double C1 = ei*ei*cos(latitude1)*cos(latitude1);

double v1 = a / sqrt(1 - e*e*sin(latitude1)*sin(latitude1));

double D = (east - FE) / (v1*k0);

double p1 = a*(1-e*e) /

((1-e*e*sin(latitude1)*sin(latitude1))*sqrt(1-e*e*sin(latitude1)*sin(latitude1)));

double latitude_out = latitude1 - (v1*tan(latitude1)/p1) * (D*D/2 -

(5+3*T1+10*C1-4*C1*C1-9*ei*ei)*D*D*D*D/24

+

(61+90*T1+298*C1+45*T1*T1-252*ei*ei-3*C1*C1)*D*D*D*D*D*D/720);

double longtitude_out = longtitude0 + (D - (1+2*T1+C1)*D*D*D/6 +

(5-2*C1+28*T1-3*C1+8*ei*ei

+ 24*T1*T1)*D*D*D*D*D/120)/cos(latitude1);

cout << "latitude = "

<< latitude_out*180/pi

<< endl;

cout << "longtitude = "

<< longtitude_out*180/pi

<< endl;

char key;

cout << "Quit? y/n";

cin >> key;

if (key == 'y')

break;

}

分享：

喜欢

0

赠金笔

加载中，请稍候......

评论加载中，请稍候...

发评论

登录名： 密码： 找回密码 注册记住登录状态

昵   称：

评论并转载此博文

发评论

以上网友发言只代表其个人观点，不代表新浪网的观点或立场。