白塞尔_白塞尔大地主题解算正反解_

// 白塞尔主题解算正反解.cpp: 定义控制台应用程序的入口点。 // #include"stdafx.h" #include<iostream> #include<cmath> using namespace std; #define e2 0.006693421622966//克拉索夫斯基椭球元素 #define PI 3.141592653 double hudu(double, double, double); void zhengsuan(double, double, double, double); void fansuan(double, double, double, double); void main(void) { int k; cout << "该主题解算用的是克拉索夫斯基椭球元素" << endl; cout << "请选择正算还是反算,正算输入1,反算输入2:" << endl; cin >> k; if (k == 1) { double S; double m1, m2, m3, d1, f1, d2, f2, d3, f3; double B1, L1, A1; cout << "请按次序输入A1点的大地纬度(B1)的度分秒，中间用空格隔开:" << endl; cin >> d1 >> f1 >> m1; cout << "请按次序输入A1点的大地经度(L1)的度分秒，中间用空格隔开:" << endl; cin >> d2 >> f2 >> m2; cout << "请按次序输入A1点的大地方位角(A1)的度分秒，中间用空格隔开:" << endl; cin >> d3 >> f3 >> m3; cout << "请输入大地线长度S(m):" << endl; cin >> S; B1 = hudu(d1, f1, m1); L1 = hudu(d2, f2, m2); A1 = hudu(d3, f3, m3); zhengsuan(B1, L1, A1, S); } else if (k == 2) { double B1, B2, L1, L2; double m1, m2, m3, m4, d1, d2, d3, d4, f1, f2, f3, f4; cout << "请按次序输入起点A1点的大地纬度(B1)的度分秒，中间用空格隔开:" << endl; cin >> d1 >> f1 >> m1; cout << "请按次序输入起点A1点的大地经度(L1)的度分秒，中间用空格隔开:" << endl; cin >> d2 >> f2 >> m2; cout << "请按次序输入终点A2点的大地纬度(B2)的度分秒，中间用空格隔开:" << endl; cin >> d3 >> f3 >> m3; cout << "请按次序输入终点A2点的大地经度(L2)的度分秒，中间用空格隔开:" << endl; cin >> d4 >> f4 >> m4; B1 = hudu(d1, f1, m1); L1 = hudu(d2, f2, m2); B2 = hudu(d3, f3, m3); L2 = hudu(d4, f4, m4); fansuan(B1, L1, B2, L2); } else printf("输入错误！\n"); } //度分秒转换函数 double hudu(double d, double f, double m) { double D; if (d >= 0) D = (d + f / 60 + m / 3600)*PI / 180; else D = -(fabs(d) + fabs(f / 60) + fabs(m / 3600))*PI / 180; return D; } //正算函数 void zhengsuan(double B1, double L1, double A1, double S) { //σ=Q,μ=u,α=a,β=b,δ=r,λ=v double W1, B2, L2, u1, u2, A0, A2, Q1, Q0, Q, a, A, B, C, b, r, v; int B2D, B2F, A2D, A2F, L2D, L2F; double B2M, A2M, L2M; W1 = sqrt(1 - e2 * sin(B1)*sin(B1)); u1 = acos(cos(B1) / W1); A0 = asin(cos(u1)*sin(A1)); Q1 = atan(sin(u1) / (cos(u1)*cos(A1))); A = 6356863.02 + (10708.949 - 13.474*cos(A0)*cos(A0))*cos(A0)*cos(A0); B = (5354.469 - 8.978*cos(A0)*cos(A0))*cos(A0)*cos(A0); C = (2.238*cos(A0)*cos(A0))*cos(A0)*cos(A0) + 0.006; a = 691.46768 - (0.58143 - 0.00144*cos(A0)*cos(A0))*cos(A0)*cos(A0); b = (0.2907 - 0.001*cos(A0)*cos(A0))*cos(A0)*cos(A0); Q0 = (S - (B + C * cos(2.0*Q1))*sin(2.0*Q1)) / A; Q = Q0 + (B + 5.0*C*cos(2.0*(Q1 + Q0)))*sin(2.0*(Q1 + Q0)) / A; r = (a*Q + b * (sin(2 * (Q1 + Q0)) - sin(2 * Q1)))*sin(A0); r = r / 3600; //计算终点大地坐标及大地方位角 u2 = asin(sin(u1)*cos(Q) + cos(u1)*cos(A1)*sin(Q)); B2 = atan(sin(u2) / (sqrt(1 - e2)*sqrt(1 - sin(u2)*sin(u2)))); v = atan(sin(A1)*sin(Q) / (cos(u1)*cos(Q) - sin(u1)*sin(Q)*cos(A1))); if (sin(A1)>0 && tan(v)>0) v = fabs(v); else if (sin(A1)>0 && tan(v)<0) v = PI - fabs(v); else if (sin(A1)<0 && tan(v)<0) v = 0 - fabs(v); else v = fabs(v) - PI; L2 = L1 + v - r * PI / 180; A2 = atan((cos(u1)*sin(A1)) / (cos(u1)*cos(Q)*cos(A1) - sin(u1)*sin(Q))); if (sin(A1)<0 && tan(A2)>0) A2 = fabs(A2); else if (sin(A1)<0 && tan(A2)<0) A2 = PI - fabs(A2); else if (sin(A1)>0 && tan(A2)>0) A2 = PI + fabs(A2); else A2 = 2 * PI - fabs(A2); B2D = (int)(B2 * 180 / PI); B2F = (int)((fabs(B2 * 180 / PI) - (double)fabs(B2D)) * 60); B2M = ((fabs(B2 * 180 / PI) - (double)fabs(B2D)) * 60 - (double)B2F) * 60; cout << "大地线终点的纬度B2为:" << B2D << "゜" << B2F << "′" << B2M << "″" << endl; L2D = (int)(L2 * 180 / PI); L2F = (int)((fabs(L2 * 180 / PI) - (double)fabs(L2D)) * 60); L2M = ((fabs(L2 * 180 / PI) - (double)fabs(L2D)) * 60 - (double)L2F) * 60; cout << "大地线终点的经度L2为:" << L2D << "゜" << L2F << "′" << L2M << "″" << endl; A2D = (int)(A2 * 180 / PI); A2F = (int)((fabs(A2 * 180 / PI) - (double)fabs(A2D)) * 60); A2M = ((fabs(A2 * 180 / PI) - (double)fabs(A2D)) * 60 - (double)A2F) * 60; cout << "大地方�