反算计算公式不会写 ！有没有测绘方面的编程厉害选手 出手帮帮我 ꒦ິ^꒦ິ已经有一堆数据并且是读取到表格控件之后再进行计算
https://blog.csdn.net/weixin_38154901/article/details/111768365
是指以改正数 Vi 代替真误差 △i 得到按观测值的改正值计算观测值的中误差公式，这个么？
这这这这
可以试试下面的代码参考
using System;
class Program {
static void Main() {
double[] r0 = { 7000, 0, 0 }; // 初始位置矢量
double[] v0 = { 0, 7.5, 7.5 }; // 初始速度矢量
double mu = 398600; // 地球引力常数
double[] r = new double[3];
double[] v = new double[3];
double a, e, i, Omega, omega, M;
double[] h = new double[3];
double[] N = new double[3];
double[] eVec = new double[3];
double eMag, rMag, vMag, vr, vp, alpha, beta;
// 计算角动量矢量
h[0] = r0[1] * v0[2] - r0[2] * v0[1];
h[1] = r0[2] * v0[0] - r0[0] * v0[2];
h[2] = r0[0] * v0[1] - r0[1] * v0[0];
// 计算升交点赤经矢量
N[0] = -h[1];
N[1] = h[0];
N[2] = 0;
// 计算轨道参数
rMag = Math.Sqrt(r0[0] * r0[0] + r0[1] * r0[1] + r0[2] * r0[2]);
vMag = Math.Sqrt(v0[0] * v0[0] + v0[1] * v0[1] + v0[2] * v0[2]);
vr = (r0[0] * v0[0] + r0[1] * v0[1] + r0[2] * v0[2]) / rMag;
vp = Math.Sqrt(vMag * vMag - vr * vr);
eVec[0] = (v0[1] * h[2] - v0[2] * h[1]) / mu - r0[0] / rMag;
eVec[1] = (v0[2] * h[0] - v0[0] * h[2]) / mu - r0[1] / rMag;
eVec[2] = (v0[0] * h[1] - v0[1] * h[0]) / mu - r0[2] / rMag;
eMag = Math.Sqrt(eVec[0] * eVec[0] + eVec[1] * eVec[1] + eVec[2] * eVec[2]);
a = 1 / (2 / rMag - vMag * vMag / mu);
e = eMag;
i = Math.Acos(h[2] / Math.Sqrt(h[0] * h[0] + h[1] * h[1] + h[2] * h[2]));
Omega = Math.Acos(N[0] / Math.Sqrt(N[0] * N[0] + N[1] * N[1]));
if (N[1] < 0) {
Omega = 2 * Math.PI - Omega;
}
omega = Math.Acos((N[0] * eVec[0] + N[1] * eVec[1] + N[2] * eVec[2]) / (N[0] * N[0] + N[1] * N[1]));
if (eVec[2] < 0) {
omega = 2 * Math.PI - omega;
}
M = Math.Acos((eVec[0] * r0[0] + eVec[1] * r0[1] + eVec[2] * r0[2]) / (eMag * rMag));
if (vr < 0) {
M = 2 * Math.PI - M;
}
// 输出轨道参数
Console.WriteLine("Semi-major axis (a): {0} km", a);
Console.WriteLine("Eccentricity (e): {0}", e);
Console.WriteLine("Inclination (i): {0} deg", i * 180 / Math.PI);
Console.WriteLine("Right ascension of the ascending node (Omega): {0} deg", Omega * 180 / Math.PI);
Console.WriteLine("Argument of perigee (omega): {0} deg", omega * 180 / Math.PI);
Console.WriteLine("Mean anomaly at epoch (M): {0} deg", M * 180 / Math.PI);
}
}
引用网络资料指引作答：
算法原理：
已知两端点经纬度，反求大地线长度及正反大地方位角。
算法流程：
首先输入两端点经纬度，将角度转化为弧度参与后续计算。然后输入椭球参数（本次使用克氏椭球），计算辅助值。再用逐次趋近法同时计算起点大地方位角、球面长度及经差（第一次趋近时，取δ=0），直到最后两次δ小于限值（此处考虑到计算机省略小数点后几位的特性，直接令两者相等）。然后计算辅助系数，之后计算大底线长度，计算反方位角。
#include <iostream>
#include<cmath>
using namespace std;
int main()
{
double angle_to_radian(double degree, double min, double second);
void radian_to_angle(double rad);
double B1,L1,B2,L2,B1d, B1f, B1m, L1d, L1f, L1m, B2d, B2f, B2m, L2d, L2f, L2m;
cout << "请输入大地线起点坐标B1=" << "°" << "′" << "″" << endl;
cin >> B1d >> B1f >> B1m;
cout << "经度L1=" << "°" << "′" << "″" << endl;
cin >> L1d >> L1f >> L1m;
cout << "请输入大地线终点坐标B2=" << "°" << "′" << "″" << endl;
cin >> B2d >> B2f >> B2m;
cout << "经度L2=" << "°" << "′" << "″" << endl;
cin >> L2d >> L2f >> L2m;
B1 = angle_to_radian(B1d, B1f, B1m);
L1 = angle_to_radian(L1d, L1f, L1m);
B2 = angle_to_radian(B2d, B2f, B2m);
L2 = angle_to_radian(L2d, L2f, L2m);
//辅助计算
double W1, W2, sinu1, sinu2, cosu1, cosu2, L, a1, a2, b1, b2;
double e2= 0.006693421622966;//克拉索夫斯基椭球
W1 = sqrt(1 - e2 * sin(B1) * sin(B1));
W2 = sqrt(1 - e2 * sin(B2) * sin(B2));
sinu1 = sin(B1) * sqrt(1 - e2) / W1;
cosu1=cos(B1) / W1;
sinu2 = sin(B2) * sqrt(1 - e2) / W2;
cosu2 = cos(B2) / W2;
L = L2 - L1;
a1 = sinu1 * sinu2;
a2 = cosu1 * cosu2;
b1 = cosu1 * sinu2;
b2 = sinu1 * cosu2;
//逐次趋近法求起点大地方位角、球面长度及经差
double l, d = 0, p, q, A1, c,sinc,cosc,t,i=0;
double sinA0, x, a, b;
double pi = 4 * atan(1);
do{
t = d;
l = L + d;
p = cosu2 * sin(l);
q = b1 - b2 * cos(l);
A1 = atan(p / q);
if (p > 0)
{
if (q > 0)
A1 = fabs(A1);
else
A1 = pi - fabs(A1);
}
else
{
if (q < 0)
A1 = pi + fabs(A1);
else
A1 = 2 * pi - fabs(A1);
}
sinc = p * sin(A1) + q * cos(A1);
cosc = a1 + a2 * cos(l);
c = atan(sinc / cosc);
if (cosc > 0)
c = fabs(c);
else
c = pi - fabs(c);
sinA0 = cosu1 * sin(A1);
x = 2 * a1 - (1 - sinA0 * sinA0) * cosc;
a = (33523299 - (28189 - 70 * (1 - sinA0 * sinA0)) * (1 - sinA0 * sinA0)) * 1e-10;
b = (28189 - 94 * (1 - sinA0 * sinA0)) * 1e-10;
d = (a * c - b * x * sinc) * sinA0;
i = i + 1;
/* cout << "d"<<i<<" = " << d<<endl<<"l"<<i<<"="<<l << endl << "A1" << i << "=" << A1 << endl << "q" << i << "=" << q<<endl;*/
} while (t != d);
//计算A,B''，C'',S,A2
double A,B_,C_,y,S,A2;
double cos2A0 = 1 - sinA0 * sinA0;
A = 6356863.020 + (10708.949 - 13.474 * cos2A0) * cos2A0;
B_ = 10708.938 - 17.956 * cos2A0;
C_ = 4.487;
y = (cos2A0 * cos2A0 - 2 * x * x) * cosc;
S = A * c + (B_ * x + C_ * y) * sinc;
A2 = atan(cosu1 * sin(l) / (b1 * cos(l) - b2));
if (p < 0)
{
if (q< 0)
A2 = fabs(A2);
else
A2 = pi - fabs(A2);
}
else
{
if (q> 0)
A2 = pi + fabs(A2);
else
A2 = 2 * pi - fabs(A2);
}//A2符号判定有问题，待修改
cout << "S=" << S << endl;
cout << "A1=";
radian_to_angle(A1);
cout<< "A2=" ;
radian_to_angle(A2);
}
double angle_to_radian(double degree, double min, double second)
{
double flag = (degree < 0) ? -1.0 : 1.0; //判断正负
double pi = 4 * atan(1);
if (degree < 0)
{
degree = degree * (-1.0);
}
double angle = degree + min / 60 + second / 3600;
double result = flag * (angle * pi) / 180;
return result;
//cout<<result<<endl;eg
}
void radian_to_angle(double rad)
{
double flag = (rad < 0) ? -1.0 : 1.0;
double pi = 4 * atan(1);
if (rad < 0)
{
rad = rad * (-1.0);
}
double result = (rad * 180) / pi;
double degree = int(result);
double min = (result - degree) * 60;
double second = (min - int(min)) * 60;
cout << flag * degree << "°" << int(min) << "′" << second << "″"<<endl;
}
我写了一个，你可以先参考着 看看先：
using System;
class Program
{
static void Main()
{
double a = 6378137; // 赤道半径
double b = 6356752.314245; // 极半径
double f = 1 / 298.257223563; //扁率
double lon1 = 118.5; //点1的经度(单位：度)
double lat1 = 32.5; //点1的纬度(单位：度)
double alpha1 = 30; //点1的方位角(单位：度)
double s12 = 10000; //点1到点2的距离(单位：米)
double a1 = alpha1 * Math.PI / 180; //将方位角转换为弧度
double f1 = (1 - f) * Math.Tan(lat1 * Math.PI / 180); //计算一些中间变量
double cos2a1 = Math.Cos(a1) * Math.Cos(a1);
double sin2a1 = Math.Sin(a1) * Math.Sin(a1);
double u1 = Math.Atan((1 - f) * Math.Tan(lat1 * Math.PI / 180));
double sinu1 = Math.Sin(u1);
double cosu1 = Math.Cos(u1);
double sigma1 = Math.Atan2(Math.Tan(u1), Math.Cos(a1)); //计算初始子午圈弧长度
double alpha = 0;
double sigma = s12 / (b * (1 - f) * cosu1); //计算迭代的第一个步长
double sinSigma = 0;
double cosSigma = 0;
double deltaSigma = 0;
while (Math.Abs(sigma - deltaSigma) > 1e-12) //迭代计算
{
sinSigma = Math.Sin(sigma);
cosSigma = Math.Cos(sigma);
double sin2Sigma = 2 * sinSigma * cosSigma;
double cos2Sigma = cosSigma * cosSigma - sinSigma * sinSigma;
double sinSqSigma = sinSigma * sinSigma;
double cosSqSigma = cosSigma * cosSigma;
double cos2Sigmam = cos2Sigma * cos2a1 - 2 * sinSigma * sin2a1;
double C = f / 16 * cosSqSigma * (4 + f * (4 - 3 * cosSqSigma));
deltaSigma = sigma;
sigma = s12 / (b * (1 - f) * cosu1) + C * sinSigma * (cos2Sigmam + C / 4 * (cosSigma * (-1 + 2 * cos2Sigmam * cos2Sigmam) - C / 6 * cos2Sigmam * (-3 + 4 * sinSqSigma) * (-3 + 4 * cos2Sigmam * cos2Sigmam)));
}
double tmp = sinu1 * sinSigma - cosu1 * cosSigma * Math.Cos(alpha1 * Math.PI / 180);
double lat2 = Math.Atan2(sinu1 * cosSigma + cosu1 * sinSigma * Math.Cos(alpha1 * Math.PI / 180), (1 - f) * Math.Sqrt(sin2a1 + tmp * tmp));
double lambda = Math.Atan2(sinSigma * Math.Sin(alpha1 * Math.PI / 180), cosu1 * cosSigma - sinu1 * sinSigma * Math.Cos(alpha1 * Math.PI / 180));
double C2 = f / 16 * cosSqSigma * (4 + f * (4 - 3 * cosSqSigma));
double L = lambda - (1 - C2) * f * sinSigma * (sigma + C2 * sinSigma * (cos2Sigmam + C2 * cosSigma * (-1 + 2 * cos2Sigmam * cos2Sigmam)));
double lon2 = lon1 + L * 180 / Math.PI;
Console.WriteLine("点2的纬度为：" + lat2 * 180 / Math.PI);
Console.WriteLine("点2的经度为：" + lon2);
}
}
#include <iostream>
#include<cmath>
using namespace std;
int main()
{
double angle_to_radian(double degree, double min, double second);
void radian_to_angle(double rad);
double B1,L1,B2,L2,B1d, B1f, B1m, L1d, L1f, L1m, B2d, B2f, B2m, L2d, L2f, L2m;
cout << "请输入大地线起点坐标B1=" << "°" << "′" << "″" << endl;
cin >> B1d >> B1f >> B1m;
cout << "经度L1=" << "°" << "′" << "″" << endl;
cin >> L1d >> L1f >> L1m;
cout << "请输入大地线终点坐标B2=" << "°" << "′" << "″" << endl;
cin >> B2d >> B2f >> B2m;
cout << "经度L2=" << "°" << "′" << "″" << endl;
cin >> L2d >> L2f >> L2m;
B1 = angle_to_radian(B1d, B1f, B1m);
L1 = angle_to_radian(L1d, L1f, L1m);
B2 = angle_to_radian(B2d, B2f, B2m);
L2 = angle_to_radian(L2d, L2f, L2m);
//辅助计算
double W1, W2, sinu1, sinu2, cosu1, cosu2, L, a1, a2, b1, b2;
double e2= 0.006693421622966;//克拉索夫斯基椭球
W1 = sqrt(1 - e2 * sin(B1) * sin(B1));
W2 = sqrt(1 - e2 * sin(B2) * sin(B2));
sinu1 = sin(B1) * sqrt(1 - e2) / W1;
cosu1=cos(B1) / W1;
sinu2 = sin(B2) * sqrt(1 - e2) / W2;
cosu2 = cos(B2) / W2;
L = L2 - L1;
a1 = sinu1 * sinu2;
a2 = cosu1 * cosu2;
b1 = cosu1 * sinu2;
b2 = sinu1 * cosu2;
//逐次趋近法求起点大地方位角、球面长度及经差
double l, d = 0, p, q, A1, c,sinc,cosc,t,i=0;
double sinA0, x, a, b;
double pi = 4 * atan(1);
do{
t = d;
l = L + d;
p = cosu2 * sin(l);
q = b1 - b2 * cos(l);
A1 = atan(p / q);
if (p > 0)
{
if (q > 0)
A1 = fabs(A1);
else
A1 = pi - fabs(A1);
}
else
{
if (q < 0)
A1 = pi + fabs(A1);
else
A1 = 2 * pi - fabs(A1);
}
sinc = p * sin(A1) + q * cos(A1);
cosc = a1 + a2 * cos(l);
c = atan(sinc / cosc);
if (cosc > 0)
c = fabs(c);
else
c = pi - fabs(c);
sinA0 = cosu1 * sin(A1);
x = 2 * a1 - (1 - sinA0 * sinA0) * cosc;
a = (33523299 - (28189 - 70 * (1 - sinA0 * sinA0)) * (1 - sinA0 * sinA0)) * 1e-10;
b = (28189 - 94 * (1 - sinA0 * sinA0)) * 1e-10;
d = (a * c - b * x * sinc) * sinA0;
i = i + 1;
/* cout << "d"<<i<<" = " << d<<endl<<"l"<<i<<"="<<l << endl << "A1" << i << "=" << A1 << endl << "q" << i << "=" << q<<endl;*/
} while (t != d);
//计算A,B''，C'',S,A2
double A,B_,C_,y,S,A2;
double cos2A0 = 1 - sinA0 * sinA0;
A = 6356863.020 + (10708.949 - 13.474 * cos2A0) * cos2A0;
B_ = 10708.938 - 17.956 * cos2A0;
C_ = 4.487;
y = (cos2A0 * cos2A0 - 2 * x * x) * cosc;
S = A * c + (B_ * x + C_ * y) * sinc;
A2 = atan(cosu1 * sin(l) / (b1 * cos(l) - b2));
if (p < 0)
{
if (q< 0)
A2 = fabs(A2);
else
A2 = pi - fabs(A2);
}
else
{
if (q> 0)
A2 = pi + fabs(A2);
else
A2 = 2 * pi - fabs(A2);
}//A2符号判定有问题，待修改
cout << "S=" << S << endl;
cout << "A1=";
radian_to_angle(A1);
cout<< "A2=" ;
radian_to_angle(A2);
}
double angle_to_radian(double degree, double min, double second)
{
double flag = (degree < 0) ? -1.0 : 1.0; //判断正负
double pi = 4 * atan(1);
if (degree < 0)
{
degree = degree * (-1.0);
}
double angle = degree + min / 60 + second / 3600;
double result = flag * (angle * pi) / 180;
return result;
//cout<<result<<endl;eg
}
void radian_to_angle(double rad)
{
double flag = (rad < 0) ? -1.0 : 1.0;
double pi = 4 * atan(1);
if (rad < 0)
{
rad = rad * (-1.0);
}
double result = (rad * 180) / pi;
double degree = int(result);
double min = (result - degree) * 60;
double second = (min - int(min)) * 60;
cout << flag * degree << "°" << int(min) << "′" << second << "″"<<endl;
}
白塞尔公式是一种计算椭圆轨道位置的数学公式，它可以用来计算行星、卫星等在轨道上的位置。下面是 C# 语言的白塞尔公式反算程序示例：
using System;
class Program
{
static void Main(string[] args)
{
double a = 6378137.0; // 椭球长轴
double b = 6356752.314245; // 椭球短轴
double e = Math.Sqrt(1 - (b / a) * (b / a)); // 第一偏心率
double e2 = Math.Sqrt((a / b) * (a / b) - 1); // 第二偏心率
double f = 1 / 298.257223563; // 椭球扁率
double longitude = 116.397458591; // 经度
double latitude = 39.908714904; // 纬度
double height = 0.0; // 高度
double N = a / Math.Sqrt(1 - e * e * Math.Sin(latitude * Math.PI / 180) * Math.Sin(latitude * Math.PI / 180));
double X = (N + height) * Math.Cos(latitude * Math.PI / 180) * Math.Cos(longitude * Math.PI / 180);
double Y = (N + height) * Math.Cos(latitude * Math.PI / 180) * Math.Sin(longitude * Math.PI / 180);
double Z = ((1 - e * e) * N + height) * Math.Sin(latitude * Math.PI / 180);
double p = Math.Sqrt(X * X + Y * Y);
double theta = Math.Atan(Z * a / (p * b));
double longitude0 = Math.Atan(Y / X);
double latitude0 = Math.Atan((Z + e2 * e2 * b * Math.Pow(Math.Sin(theta), 3)) / (p - e * e * a * Math.Pow(Math.Cos(theta), 3)));
double N0 = a / Math.Sqrt(1 - e * e * Math.Sin(latitude0) * Math.Sin(latitude0));
double height0 = p / Math.Cos(latitude0) - N0;
Console.WriteLine("经度：{0}，纬度：{1}，高度：{2}", longitude0 * 180 / Math.PI, latitude0 * 180 / Math.PI, height0);
}
}
这段代码可以计算出给定经纬度的点在椭球上的位置，具体的计算过程请参考代码注释。注意，这段代码中的参数 a、b、e、f 等是根据 WGS84 椭球参数进行计算的，如果使用其他椭球参数，需要相应地修改这些参数。
csharp
using System;
class BezierCurve
{
// 计算阶乘
static int Factorial(int n)
{
if (n == 0)
{
return 1;
}
else
{
return n * Factorial(n - 1);
}
}
// 计算组合数
static int Combination(int n, int i)
{
return Factorial(n) / (Factorial(i) * Factorial(n - i));
}
// 计算贝塞尔曲线上的点
static void CalculateBezierPoint(double t, double[] x, double[] y, out double px, out double py)
{
int n = x.Length - 1;
px = 0.0;
py = 0.0;
for (int i = 0; i <= n; i++)
{
double b = Combination(n, i) * Math.Pow(t, i) * Math.Pow(1 - t, n - i);
px += b * x[i];
py += b * y[i];
}
}
static void Main()
{
double[] x = { 0.0, 2.0, 4.0, 6.0, 8.0 };
double[] y = { 1.0, 3.0, 4.0, 2.0, 0.0 };
double px, py;
CalculateBezierPoint(0.5, x, y, out px, out py);
Console.WriteLine("Bezier curve point at t=0.5: ({0}, {1})", px, py);
}
}
在上述代码中，Factorial函数用于计算阶乘，Combination函数用于计算组合数。CalculateBezierPoint函数用于计算贝塞尔曲线上的点，其中t为参数，x和y为控制点数组，px和py为输出的贝塞尔曲线上的点。Main函数用于测试，在控制点数组x和y中给定5个点，计算t=0.5时贝塞尔曲线上的点，并输出结果。
C#语言实现的白塞尔公式反算程序参考代码
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Drawing;
namespace BezierCurve
{
class Program
{
static void Main(string[] args)
{
List<Point> points = new List<Point>();
points.Add(new Point(20, 20));
points.Add(new Point(50, 100));
points.Add(new Point(120, 60));
points.Add(new Point(180, 180));
List<Point> controlPoints = CalculateControlPoints(points);
Console.WriteLine("控制点坐标：");
foreach (var point in controlPoints)
{
Console.WriteLine("({0}, {1})", point.X, point.Y);
}
Console.ReadLine();
}
static List<Point> CalculateControlPoints(List<Point> points)
{
int n = points.Count - 1;
List<Point> controlPoints = new List<Point>();
// 初始化m和A
int m = 2 * n;
double[,] A = new double[m + 1, m + 2];
for (int i = 0; i <= m; i++)
{
for (int j = 0; j <= n; j++)
{
if (i >= j && i < j + n + 1)
{
A[i, j] = Math.Pow(i - j, n);
}
}
}
// 初始化b
for (int i = 0; i <= n; i++)
{
A[i + n + 1, m + 1] = points[i].X;
A[i + n + 1, m + 2] = points[i].Y;
}
// 高斯-约旦消元求解
for (int i = 0; i <= m; i++)
{
double temp = A[i, i];
for (int j = i; j <= m + 2; j++)
{
A[i, j] /= temp;
}
for (int j = 0; j <= m; j++)
{
if (j != i)
{
temp = A[j, i];
for (int k = i; k <= m + 2; k++)
{
A[j, k] -= temp * A[i, k];
}
}
}
}
// 提取控制点
for (int i = 0; i <= n; i++)
{
controlPoints.Add(new Point(Convert.ToInt32(A[i + n + 1, m + 1]), Convert.ToInt32(A[i + n + 1, m + 2])));
}
return controlPoints;
}
}
}