一、FileData
/// <summary>
/// .txt椭球参数及经纬度
/// </summary>
class Dat_Geo
{
/// <summary>
/// 终点起点经纬度
/// </summary>
public static double B1;
public static double B2;
public static double L1;
public static double L2;
/// <summary>
/// 椭球参数
/// </summary>
public static double a;
public static double b;
public static double f;
public static double e2;
public static double e2_;
}
/// <summary>
/// 中间计算数据
/// </summary>
class Dat_Mid
{
public static double l;//初始经差
public static double a1;//规划经纬度
public static double a2;
public static double b1;
public static double b2;
public static double u1;
public static double u2;
public static double A1;//起始大地方位角
public static double sigema;
public static double sinA02;
public static double cosA02;
public static double S;//大地线长度
}
二、Cal_Aux
/// <summary>
/// 辅助计算
/// </summary>
class Cal_Aux
{
/// <summary>
/// 计算偏心率
/// </summary>
/// <param name="dat_Geo">.txt椭球数据类</param>
public static void Cal_Ellipsoid()
{
Dat_Geo.b = Dat_Geo.a - Dat_Geo.a * Dat_Geo.f;
Dat_Geo.e2 = (Math.Pow(Dat_Geo.a, 2) - Math.Pow(Dat_Geo.b, 2))/ Math.Pow(Dat_Geo.a, 2);
Dat_Geo.e2_ =Dat_Geo.e2/(1 - Dat_Geo.e2);
}
/// <summary>
/// dd.mmss2rad
/// </summary>
/// <param name="dms">dd.mmss</param>
/// <returns></returns>
public static double Dms2Rad(double dms)
{
double deg = (int)dms;
double min = (int)((dms - deg) * 100);
double sec = dms*10000-deg*10000-min*100;
double rad = (deg + min/60.0 +sec/3600.0)*Math.PI/180.0;
return rad;
}
/// <summary>
/// rad2dd.mmss
/// </summary>
/// <param name="rad"></param>
/// <returns></returns>
public static string Rad2dms(double rad)
{
double dd = rad * 180.0 / Math.PI;
double deg = (int)dd;
double min = (int)((dd - deg) * 60);
double sec = (dd - deg - min / 60.0) * 3600;
string dms = deg.ToString("f0") + '.' + min.ToString("f0") + Math.Floor(sec).ToString("f0") + ((sec-Math.Floor(sec))*10000).ToString("f0");
return dms;
}
/// <summary>
/// 判断大地方位角象限
/// </summary>
/// <param name="p"></param>
/// <param name="q"></param>
/// <param name="A">大地方位角</param>
/// <returns></returns>
public static double Cal_Quad_A(double p, double q, double A)
{
double A_;
//象限角判断
if (p > 0 && q > 0)
{
A_ = Math.Abs(A);
}
else if (p > 0 && q < 0)
{
A_ = Math.PI-Math.Abs(A);
}
else if (p < 0 && q < 0)
{
A_ = Math.PI + Math.Abs(A);
}
else
{
A_ = 2*Math.PI - Math.Abs(A);
}
//判断A_是否属于0-2Π
if (A_ < 0)
{
A_ += 2 * Math.PI;
}
else if (A_ > 2 * Math.PI)
{
A_ -= 2 * Math.PI;
}
return A_;
}
/// <summary>
/// sigema象限改正
/// </summary>
/// <param name="sigema"></param>
/// <returns></returns>
public static double Cal_Quad_sigema(double sigema)
{
double sig;
if (Math.Cos(sigema ) > 0)
{
sig = Math.Abs(sigema);
}
else
{
sig = Math.PI - Math.Abs(sigema);
}
return sig;
}
}
三、Geoline
/// <summary>
/// 大地线解算
/// </summary>
class GeoLine
{
/// <summary>
/// 辅助计算
/// </summary>
/// <param name="dat_Geo">椭球参数类</param>
/// <param name="dat_Mid">中间运算参数类</param>
public static void Geo_Aux()
{
Dat_Mid.u1 = Math.Atan(Math.Pow(1-Dat_Geo.e2,0.5)*Math.Tan(Dat_Geo.B1));
Dat_Mid.u2 = Math.Atan(Math.Pow(1 - Dat_Geo.e2,0.5) * Math.Tan(Dat_Geo.B2));
Dat_Mid.l = Dat_Geo.L2 - Dat_Geo.L1;
Dat_Mid.a1 = Math.Sin(Dat_Mid.u1)*Math.Sin(Dat_Mid.u2);
Dat_Mid.a2 = Math.Cos(Dat_Mid.u1) * Math.Cos(Dat_Mid.u2);
Dat_Mid.b1 = Math.Cos(Dat_Mid.u1) * Math.Sin(Dat_Mid.u2);
Dat_Mid.b2 = Math.Sin(Dat_Mid.u1) * Math.Cos(Dat_Mid.u2);
}
/// <summary>
/// 计算起点大地方位角
/// </summary>
/// <param name="dat_Geo"></param>
/// <param name="dat_Mid"></param>
public static void Cal_A1()
{
// 变量
double p;
double q;
double delta0 = 0;
double delta;
double lamda = Dat_Mid.l + delta0;
double sigema1;
double sinSig;
double cosSig;
double alpha;
double belta;
double galma;
/*
逐次逼近法计算起点大地方位角A1和经差lamda
*/
do
{
p = Math.Cos(Dat_Mid.u2 * Math.Sin(lamda));
q = Dat_Mid.b1 - Dat_Mid.b2 * Math.Cos(lamda);
Dat_Mid.A1 = Cal_Aux.Cal_Quad_A(p, q, Math.Atan(p / q)); //象限角度归化
sinSig = p * Math.Sin(Dat_Mid.A1) + q * Math.Cos(Dat_Mid.A1);
cosSig = Dat_Mid.a1 + Dat_Mid.a2 * Math.Cos(lamda);
Dat_Mid.sigema = Cal_Aux.Cal_Quad_sigema(Math.Atan(sinSig / cosSig)); //象限角度归化
Dat_Mid.sinA02 = Math.Pow(Math.Cos(Dat_Mid.u1) * Math.Sin(Dat_Mid.A1), 2);
Dat_Mid.cosA02 = 1 - Dat_Mid.sinA02;
sigema1 = Math.Atan(Math.Tan(Dat_Mid.u1) / Math.Cos(Dat_Mid.A1));
// alpha belta galma
alpha = (Dat_Geo.e2 / 2.0 + Math.Pow(Dat_Geo.e2, 2) / 8.0 + Math.Pow(Dat_Geo.e2, 3) / 16.0)
- (Math.Pow(Dat_Geo.e2, 2) / 16.0 + Math.Pow(Dat_Geo.e2, 3) / 16.0) * Math.Pow(Dat_Mid.cosA02, 2)
+ 3 * Math.Pow(Dat_Mid.cosA02, 4) * Math.Pow(Dat_Geo.e2, 3) / 128;
belta = (Math.Pow(Dat_Geo.e2, 2) / 16.0 + Math.Pow(Dat_Geo.e2, 3) / 16.0) * Dat_Mid.cosA02
- (Math.Pow(Dat_Geo.e2, 3) / 32.0) * Math.Pow(Dat_Mid.cosA02, 2);
galma = (Math.Pow(Dat_Geo.e2, 3) / 256.0) * Math.Pow(Dat_Mid.cosA02, 2);
// delta
delta = (alpha * Dat_Mid.sigema + belta * Math.Cos(2 * sigema1 + Dat_Mid.sigema) * Math.Sin(Dat_Mid.sigema)
+ galma * Math.Sin(2 * Dat_Mid.sigema) * Math.Cos(4 * sigema1 + 2 * Dat_Mid.sigema)) * (Math.Cos(Dat_Mid.u1) * Math.Sin(Dat_Mid.A1));
//lamda
lamda = Dat_Mid.l + delta;
if (Math.Abs(delta - delta0) < 1.0e-10)
break;
delta0 = delta;
}
while (true);
}
/// <summary>
/// 大地线解算
/// </summary>
/// <param name="dat_Geo"></param>
/// <param name="dat_Mid"></param>
public static void Cal_S_line()
{
double sigema1;
double x;
double k2;
double A, B, C;
sigema1 = Math.Atan(Math.Tan(Dat_Mid.u1) / Math.Cos(Dat_Mid.A1));
//计算系数ABC
k2 = Dat_Geo.e2_ * Dat_Mid.cosA02;
A = (1 - k2/4 + 7*k2*k2/64 - 15*k2*k2*k2/256) / Dat_Geo.b;
B = k2 / 4 - k2 * k2 / 8 + 37 * k2 * k2 * k2 / 512;
C = k2 * k2 / 128 - k2 * k2 * k2 / 128;
x = C * Math.Sin(2 * Dat_Mid.sigema) * Math.Cos(4*sigema1+2*Dat_Mid.sigema);
//大地线长度S
Dat_Mid.S = (Dat_Mid.sigema - B * Math.Sin(Dat_Mid.sigema) * Math.Cos(2 * sigema1 + Dat_Mid.sigema) - x) /A;
}
}
四、Form1
public partial class Form1 : Form
{
// 计算后才能进行报告操作
public static bool call =false;
public static bool call_1 = false;
// 储存报告
public static string report;
public Form1()
{
InitializeComponent();
}
// 初始化窗口显示
private void Form1_Load(object sender, EventArgs e)
{
//定义datagridview初始状态
dataGridView_out.RowCount = 1;
dataGridView_input.RowCount = 2;
dataGridView_input[0, 0].Value = "P1";
dataGridView_input[0, 1].Value = "P2";
}
// 文件-打开 读取数据
private void 打开ToolStripMenuItem_Click(object sender, EventArgs e)
{
// 临时储存读取文件数据
string data_all_lines = null;
// 打开文件选择器
OpenFileDialog op = new OpenFileDialog();
op.Filter = "txt|*.txt";
// 选择文件
if(op.ShowDialog() == DialogResult.OK)
{
try
{
using (StreamReader reader = new StreamReader(op.FileName))
{
data_all_lines = reader.ReadToEnd();
}
}
catch
{
MessageBox.Show("文件选取错误，请重新选取！");
状态toolStripStatusLabel2.Text = "文件选取错误，请重新选取！";
}
}
else
{
return;
}
// 处理读取数据
try
{
string[] tem = data_all_lines.Split(new char[] { ' ', '\r','\n'},StringSplitOptions.RemoveEmptyEntries);
foreach (string item in tem)
{
if (string.IsNullOrEmpty(item))
throw new Exception();
}
// 接收数据元素
// 椭球参数 a,f
Dat_Geo.a = double.Parse(tem[0].Trim());
Dat_Geo.f = 1.0/double.Parse(tem[1].Trim());
// 经纬度BL读取
Dat_Geo.B1 = double.Parse(tem[2].Trim());dataGridView_input[1,0].Value = double.Parse(tem[2].Trim());
Dat_Geo.L1 = double.Parse(tem[3].Trim()); dataGridView_input[2, 0].Value = double.Parse(tem[3].Trim());
Dat_Geo.B2 = double.Parse(tem[4].Trim()); dataGridView_input[1, 1].Value = double.Parse(tem[4].Trim());
Dat_Geo.L2 = double.Parse(tem[5].Trim()); dataGridView_input[2, 1].Value = double.Parse(tem[5].Trim());
call_1 = true;
状态toolStripStatusLabel2.Text = "打开文件成功！";
}
catch(Exception)
{
MessageBox.Show("文件格式错误！");
状态toolStripStatusLabel2.Text = "文件格式错误！";
}
}
// 文件-保存
// 制作报告string
public string reporter(bool call)
{
string report = null;
report = "-------------------------计算报告-------------------------\r\n";
report += "P1点坐标（B,L,dd.mmss）:"+dataGridView_input[1,0].Value.ToString()+"\t "+ dataGridView_input[2,0].Value.ToString()+"\t\r\n";
report += "P2点坐标（B,L,dd.mmss）:"+ dataGridView_input[1, 1].Value.ToString() + "\t " + dataGridView_input[2, 1].Value.ToString() + "\t\r\n";
report += "-------------------------计算结果-------------------------\r\n";
report += "大地线S长度(m)；"+dataGridView_out[0,0].Value.ToString()+" \r\n";
return report;
}
private void 保存报告ToolStripMenuItem_Click(object sender, EventArgs e)
{
if (call_1 == true)
{
SaveFileDialog sa = new SaveFileDialog();
sa.Filter = "txt|*.txt";
if (sa.ShowDialog() == DialogResult.OK)
{
try
{
using (StreamWriter sw = new StreamWriter(sa.FileName))
{
sw.WriteLine(reporter(call));
MessageBox.Show("保存成功！");
状态toolStripStatusLabel2.Text = "保存成功！";
}
}
catch
{
MessageBox.Show("保存失败！");
状态toolStripStatusLabel2.Text = "保存失败！";
}
}
else
{
return;
}
}
else
{
MessageBox.Show("请计算数据！");
}
}
// 清空数据
private void 清除数据ToolStripMenuItem_Click(object sender, EventArgs e)
{
dataGridView_input.Rows.Clear();
dataGridView_out.Rows.Clear();
call = false;
report = null;
dataGridView_out.RowCount = 1;
dataGridView_input.RowCount = 2;
dataGridView_input[0, 0].Value = "P1";
dataGridView_input[0, 1].Value = "P2";
状态toolStripStatusLabel2.Text = "清除成功！";
}
private void 计算ToolStripMenuItem_Click(object sender, EventArgs e)
{
if (call_1 == true)
{
// 计算偏心率
Cal_Aux.Cal_Ellipsoid();
// 辅助计算
GeoLine.Geo_Aux();
// 计算大地方位角
GeoLine.Cal_A1();
// 计算大地线长度
GeoLine.Cal_S_line();
// 输出结果
dataGridView_out[0, 0].Value = Dat_Mid.S;
// 保存报告
call = true;
report = reporter(call);
// 显示结果
MessageBox.Show("计算成功！");
状态toolStripStatusLabel2.Text = "计算成功！";
// 初始化
call_1 = false;
}
else
{
MessageBox.Show("请打开文件！");
}
}
private void 退出ToolStripMenuItem_Click(object sender, EventArgs e)
{
System.Environment.Exit(0);
}
private void 显示报告ToolStripMenuItem_Click(object sender, EventArgs e)
{
if(call == true)
{
Report rp = new Report();
rp.ShowDialog();
}
else
{
MessageBox.Show("请先计算数据！");
}
}
}
该代码实现了一个计算两点间大地线长度的程序，涉及到椭球参数的计算、经纬度转换以及大地方位角的确定。程序首先读取.txt文件中的椭球数据和经纬度，然后进行偏心率计算、辅助计算和大地线解算，最终得出大地线长度并生成计算报告。
扫一扫
843
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
