大地主题解算深度干货超精
白赛尔大地主题解算
dS
M
cos A c
dL
dS
sin A N cos B
V c
sec B sin B
dA
dS
tan B sin N
AV c
tan B sin
A
二阶导数：
d2B dS 2
B
( dB dS
)
dB dS
A
( dB dS
)
dA dS
V4 c2
t (3 2
10
Fundation of Geodesy
4.7.3 高斯平均引数正算公式
高斯平均引数正算公式推导 的基本思想：
首先把勒让德级数在 P１点展 开改在大地线长度中点M展开，以 使级数公式项数减少，收敛快， 精度高；其次，考虑到求定中点 M 的复杂性，将 M 点用大地线两 端点平均纬度及平均方位角相对 应的 m 点来代替，并借助迭代计 算便可顺利地实现大地主题正解。
f (BM , AM )
F (Bm
BM
Bm
,
Am
AM
Am )
+
dB ( dS )M
f f (Bm , Am ) ( Bm )(BM
f Bm ) ( Am )( AM
Am )
+
dB
dB
dB ( dS )M
( )
f (Bm , Am ) (
dS Bm
S sin Am 24 N m2
[S 2tm2
sin 2
Am
S 2 cos2 Am (1m2 9m2tm2 )]
用C语言实现大地主题解算
４
６
４
（ 争 寺＋
‘ ） ｃ ｏ ５ ４ Ａ ’
裴连磊 Ｐ Ｅ Ｉ Ｌ ｉ ａ ｎ — ｌ ｅ ｉ
（ 新疆 地 矿 局 测 绘 大 队 ， 乌鲁木齐 ８ ３ ０ ０ １ ７ ） （ Ｘ ｉ ｎ ｊ ｉ ａ ｎ ｇ Ｇ ｅ ｏ ｌ ｏ ｇ ｙ ａ ｎ ｄ Ｍｉ ｎ ｅ ｒ ａ ｌ Ｂ ｕ ｒ ｅ ａ ｕ ， Ｓ ｕ ｒ ｖ ｅ ｙ ｉ ｎ ｇ ａ ｎ ｄ Ｍａ ｐ ｐ ｉ ｎ ｇ Ｂ ｒ ｉ ｇ ａ ｄ ｅ ， Ｕ ｒ ｕ ｍｑ ｉ ８ ３ ０ ０ １ ７ ， Ｃ ｈ ｉ ｎ ａ ）
１ － ５计 算经 度 差 改 正 数
本 文的程序 实现 是采用 的 白塞 尔大地 主题 解 算的方 法， 根 据 白塞 尔大地主题解 算的 方法 ， 分析得 出适 合计 算 机编程 的大地主题正反 算的具体 实现步骤。 １ 白塞尔法大地主题正算思想 已知量 ： 大 地线 起点 的纬 度 Ｂ ， 经度 Ｌ ， 大地 方位 角 Ａ １ 及大地线长度 Ｓ 。 求解量 ： 大地线终点 的纬度 Ｂ ： ， 经度 Ｌ ２
ｃ 番 一 番 ＋ ． ．
ｃ ＿ ｂ （ 薷 ４ 一 ６ 一）
２
＝
ｃ ０ ｓ ｕ ２ ＝
Ｌ ａ １ ＝ ｓ ｉ ｎ ｕ ｌ ｓ ｉ ｎ ｕ ２ Ｉ ａ ２ ＿ ｃ ｏ ｓ ｕ ｌ ｃ ０ ｓ ｕ ２ ｂ＝
Ｃ Ｏ Ｓ Ｕｌ ｓ ｉ ｎ ｕ ２ ， ｂ ２ ＝ ｓ ｉ ｎ ｕｌ Ｃ Ｏ Ｓ Ｕ ２
之 值 。 即下 式 ：
求解 量 ： 大地 线长度 Ｓ及起、 终点 处的大地 方位角 Ａ
及 Ａ＇ 。
Ａ ： ｂ （ １ 一 鲁十 击ｋ ６ ＿ －
２ ． １辅 助计算 ｗ ：
大地主题解算深度干货超精
大地主题解算深度干货超精在我们生活的世界上，大地是我们所处的基本环境。
它是被压迫、被掠夺的，也是我们努力保护和改善的对象。
深度解算是一种精确的测量和计算方法，可以帮助我们更好地了解和利用大地资源。
本文将介绍大地主题解算和其带来的深度干货。
一、什么是大地主题解算？大地主题解算是一种基于卫星观测数据和大地测量学原理的技术，用于精确计算大地物体的三维位置和形态变化。
通过对不同时刻的卫星图像进行分析和计算，可以得到地面物体的精确坐标、高程等信息。
大地主题解算的核心思想是通过计算大地物体的三维位置和形态变化，来揭示地球表面的变化过程。
这种技术可以用于监测地表的沉降、地壳的变形、建筑物的结构和变化等。
它在测绘、地质灾害监测、城市规划等领域具有广泛的应用前景。
二、大地主题解算的应用领域1. 地质灾害监测地质灾害是世界各地普遍存在的问题，如地震、山体滑坡、地裂缝等。
通过大地主题解算技术，可以实时监测和预测地质灾害的发生和演变趋势。
这对于提前采取相应的防灾措施、保护人民的生命财产安全具有重要意义。
2. 城市规划随着城市化进程的不断加快，城市规划也变得愈发重要。
大地主题解算可以为城市规划提供精确的地面信息，包括土地利用、道路交通、建筑物布局等。
这有助于提高城市规划的科学性和有效性，减少资源浪费和环境污染。
3. 地表沉降监测地表沉降是地下水开采、矿井开采等人类活动造成的一个重要问题。
通过大地主题解算技术，可以实时监测地表的沉降情况，并对地下水开采和矿井开采等活动进行合理调整和管理。
这有助于减少地表沉降对城市建设和生态环境的不利影响。
4. 环境保护大地主题解算可以为环境保护提供准确的数据支持。
例如，通过对森林面积、湿地面积等进行监测和计算，可以及时发现和预防环境变化和破坏。
这对于保护生态环境和维护生态平衡具有重要意义。
三、大地主题解算的优势和挑战大地主题解算作为一种先进的测量和计算技术，具有很多优势。
首先，它可以实时获取和分析大地物体的位置和形态变化，具有高精度和高时效性。
分段累加法大地主题解算及高斯投影
大地测量编程实习报告班号：XXXX 学号：XXXXXXXXXX 姓名：XXX大地主题解算（分段累加法正算）结果截图：主要代码：privatevoid button1_Click(object sender, EventArgs e){double bx, by, bz, B1, lx, ly, lz, L1, ax, ay, az, A1, S;double dB, dL;double e2 = 0.006693421622966, a = 6378245.0000000000;double M, N, W, C;double B2 = 0, L2 = 0, A2 = 0;int Bx, By, Lx, Ly, Ax, Ay;double Bz, Lz, Az;bx = Convert.ToDouble(du1.Text); by = Convert.ToDouble(fen1.Text); bz = Convert.ToDouble(miao1.Text);lx = Convert.ToDouble(du2.Text); ly = Convert.ToDouble(fen2.Text); lz = Convert.ToDouble(miao2.Text); ax = Convert.ToDouble(du3.Text); ay = Convert.ToDouble(fen3.Text); az = Convert.ToDouble(miao3.Text); S = Convert.ToDouble(m.Text);double PI = Math.PI;B1 = (bx + by / 60 + bz / 3600) * PI / 180;L1 = (lx + ly / 60 + lz / 3600) * PI / 180;A1 = (ax + ay / 60 + az / 3600) * PI / 180;W = Math.Sqrt(1 - e2 * (Math.Sin(B1) * Math.Sin(B1)));M = a * (1 - e2) / (W * W * W);N = a / W;C = N * Math.Cos(B1) * Math.Sin(A1);int n ,i; double s;n = (int)(S / 4000 + 1);s = S / n;for (i = 1; i <= n; i++){dB = Math.Cos(A1) * s / M;dL = Math.Sin(A1) * s / N / Math.Cos(B1);B2 = B1 + dB;L2 = L1 + dL;W = Math.Sqrt(1 - e2 * (Math.Sin(B2) * Math.Sin(B2)));M = a * (1 - e2) / (W * W * W);N = a / W;A2 = Math.Asin(C / N / Math.Cos(B2));B1 = B2; L1 = L2; A1 = A2;}Bx =(int)( B2 / PI * 180);By = (int)((B2 / PI * 180 - Bx) * 60);Bz = Math.Abs(((B2 / PI * 180 - Bx) * 60 - By) * 60);By = Math.Abs(By);Lx = (int)(L2 / PI * 180);Ly = (int)((L2 / PI * 180 - Lx) * 60);Lz = Math.Abs(((L2 / PI * 180 - Lx) * 60 - Ly) * 60);Ly = Math.Abs(Ly);Ax = (int)(A2 / PI * 180);Ay = (int)((A2 / PI * 180 - Ax) * 60);Az = Math.Abs(((A2 / PI * 180 - Ax) * 60 - Ay) * 60);Ay = Math.Abs(Ay);du4.Text = Bx.ToString(); fen4.Text = By.ToString(); miao4.Text = Bz.ToString("0.0000"); du5.Text = Lx.ToString(); fen5.Text = Ly.ToString(); miao5.Text = Lz.ToString("0.0000"); du6.Text = Ax.ToString(); fen6.Text = Ay.ToString(); miao6.Text = Az.ToString("0.0000"); }高斯投影正、反算结果截图：正算：反算：主要代码：正算：privatevoid button1_Click(object sender, EventArgs e){double a = 0 , e12 = 0, e2 = 0;double PI = Math.PI;string tt;try{tt = boBox1.Items[boBox1.SelectedIndex].ToString();}catch{MessageBox.Show("请选择坐标系！"); return;}if (pareTo("1954年北京坐标系") == 0){a = 6378245;e2 = 0.006693421622966;e12 = 0.006738525414683;}if (pareTo("1980年国家大地坐标系") == 0){a = 6378140;e2 = 0.006694384999588;e12 = 0.006739501819473;}if (pareTo("WGS-84") == 0){a = 6378137;e2 = 0.00669437999013;e12 = 0.00673949674227;}if (pareTo("CGCS2000") == 0){a = 6378137;e2 = 0.00669438002290;e12 = 0.00673949677548;}double m0, m2, m4, m6, m8, a0, a2, a4, a6, a8, N, B, L, X, Bb, cosB ,sinB , ρ, η2, t, l, d, x, y;B = Convert.ToDouble(du1.Text) + (Convert.ToDouble(fen1.Text) / 60) + (Convert.ToDouble(miao1.Text) / 3600); L = Convert.ToDouble(du2.Text) + (Convert.ToDouble(fen2.Text) / 60) +(Convert.ToDouble(miao2.Text) / 3600);Bb = B * PI / 180;m0 = a * (1 - e2);m2 = 3 * e2 * m0 / 2;m4 = 5 * e2 * m2 / 4;m6 = 7 * e2 * m4 / 6;m8 = 9 * e2 * m6 / 8;a0 = m0 + (m2 / 2) + ((3 * m4) / 8) + ((5 * m6) / 16) + ((35 * m8) / 128);a2 = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16;a4 = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32;a6 = m6 / 32 + m8 / 16;a8 = m8 / 128;ρ = 180 * 3600 / PI;cosB = Math.Cos(Bb);sinB = Math.Sin(Bb);η2 = e12 * cosB * cosB;t = Math.Tan(Bb);d = Math.Floor(L / 6) + 1;l = Math.Abs(L - (6 * d - 3)) * 3600;N = a / Math.Sqrt(1 - e2 * sinB * sinB);X = a0 * Bb - sinB * cosB * ((a2 - a4 + a6) + (2 * a4 - 16 * a6 / 3) * sinB * sinB + 16 * a6 * sinB * sinB * sinB * sinB / 3);x = X + N * sinB * cosB * l * l / (2 * ρ * ρ) + N * sinB * cosB * cosB * cosB * (5 - t * t + 9 * η2) * l * l * l * l / (24 * ρ * ρ * ρ * ρ) + N * sinB * cosB * cosB * cosB * cosB * cosB * (61 - 58 * t * t + t * t * t * t) * l * l * l * l * l * l / (720 * ρ * ρ * ρ * ρ * ρ * ρ);y = N * cosB * l / ρ + N * cosB * cosB * cosB * (1 - t * t + η2) * l * l * l / (6 * ρ * ρ* ρ)+ N * cosB * cosB * cosB * cosB * cosB * (5 - 18 * t * t + t * t * t * t) * l * l * l * l * l/ (120 * ρ * ρ * ρ * ρ * ρ);x1.Text = x.ToString("0.0000")+"m";y1.Text = y.ToString("0.0000")+"m";}反算：privatevoid button2_Click(object sender, EventArgs e){double a = 0, e12 = 0, e2 = 0;double PI = Math.PI;string tt;try{tt = boBox2.Items[boBox2.SelectedIndex].ToString();}catch{MessageBox.Show("请选择坐标系！"); return;}if (pareTo("1954年北京坐标系") == 0){a = 6378245;e2 = 0.006693421622966;e12 = 0.006738525414683;}if (pareTo("1980年国家大地坐标系") == 0){a = 6378140;e2 = 0.006694384999588;e12 = 0.006739501819473;}if (pareTo("WGS-84") == 0){a = 6378137;e2 = 0.00669437999013;e12 = 0.00673949674227;}if (pareTo("CGCS2000") == 0){a = 6378137;e2 = 0.00669438002290;e12 = 0.00673949677548;}double x, y, m0, m2, m4, m6, m8, a0, a2, a4, a6, a8, Bf, Bf1, B , Mf , W, Nf, tf, ηf2, n2, n4, n6, n8 , l;x = Convert.ToDouble(x2.Text);y = Convert.ToDouble(y2.Text);m0 = a * (1 - e2);m2 = 3 * e2 * m0 / 2;m4 = 5 * e2 * m2 / 4;m6 = 7 * e2 * m4 / 6;m8 = 9 * e2 * m6 / 8;a0 = m0 + (m2 / 2) + ((3 * m4) / 8) + ((5 * m6) / 16) + ((35 * m8) / 128);a2 = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16;a4 = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32;a6 = m6 / 32 + m8 / 16;a8 = m8 / 128;Bf1 = 0;Bf = x / a0;while (Math.Abs(Bf - Bf1) > 1E-10){Bf1 = Bf;double sinB = Math.Sin(Bf1);double cosB = Math.Cos(Bf1);double sin2B = sinB * cosB * 2;double sin4B = sin2B * (1 - 2 * sinB * sinB) * 2;double sin6B = sin2B * Math.Sqrt(1 - sin4B * sin4B) + sin4B * Math.Sqrt(1 - sin2B * sin2B);Bf = (x - (-a2 / 2.0 * sin2B + a4 / 4.0 * sin4B - a6 / 6.0 * sin6B)) / a0;}double sinBf = Math.Sin(Bf);double cosBf = Math.Cos(Bf);ηf2 = e12 * cosBf * cosBf;tf = Math.Tan(Bf);W = Math.Sqrt(1 - e2 * sinBf * sinBf);Mf = a * (1 - e2) / (W * W * W);Nf = a / W;B = Bf - tf * y * y / (2 * Mf * Nf) + tf * (5 + 3 * tf * tf + ηf2 - 9 * ηf2 * tf * tf) * y * y * y * y / (24 * Mf * Nf * Nf * Nf) - tf * (61 + 90 * tf * tf + 45 * tf * tf * tf * tf) * y * y * y * y * y * y / (720 * Mf * Nf * Nf * Nf * Nf * Nf);l = y / (Nf * cosBf) - (1 + 2 * tf * tf + ηf2) * y * y * y / (6 * Nf * Nf * Nf * cosBf) + (5 + 28 * tf * tf + 24 * tf * tf * tf * tf + 6 * ηf2 + 8 * ηf2 * tf * tf) * y * y * y * y * y / (120 * Nf * Nf * Nf * Nf * Nf * cosBf);int Bx, By, lx, ly;double Bz, lz;Bx = (int)(B / PI * 180);By = (int)((B / PI * 180 - Bx) * 60);Bz = Math.Abs(((B / PI * 180 - Bx) * 60 - By) * 60);By = Math.Abs(By);lx = (int)(l / PI * 180);ly = (int)((l / PI * 180 - lx) * 60);lz = Math.Abs(((l / PI * 180 - lx) * 60 - ly) * 60);ly = Math.Abs(ly);du3.Text = Bx.ToString();fen3.Text = By.ToString();miao3.Text = Bz.ToString("0.0000");du4.Text = lx.ToString(); fen4.Text = ly.ToString();miao4.Text = lz.ToString("0.0000");}。
大地主题解算
大地主题解算一、实验目的：1.提高运用计算机语言编程开发的能力；2.加深对大地主题解算计算公式及辅助参数的理解并掌握计算步骤；3.通过编程语言实现大地主题解算。
二、工具：Windows XP Mode 环境下的Microsoft Visual C++ 6.0三、注意事项：1.计算所需变量多，容易混淆；2.正反算函数的编写；3.函数调用；4.弧度与角度之间的转化。
四、实验要求：1.提交报告，实验总结，编写代码；2.独立编程，调试运行；3.上交成果：编写思想，编写过程，问题分析，源代码，计算结果；五、编程过程实现：1.对白塞尔法大地主题解算有一定的了解，并参考教材P148-P150；2.由于参数较多，而在C语言环境下很多符号无法定义，需要符合要求的定义符号替代书本上那些无法直接在C语言环境下定义的符号来达到实现实验的目的；3.程序中采用弧度与度分秒之间转换的函数定义与调用，减轻一定的实验麻烦；4.在C语言环境下，数学函数fabs代替abs起绝对值作用，atan代替arctan起反函数作用；5.程序中尤其注意弧度与角度之间转换，在C语言环境下电脑默认为弧度。
六、源程序代码：#include<stdio.h>#include<math.h>double hudu(double,double,double); /*度分秒转换为弧度*/double du(double); /*弧度转换为度*/double fen(double); /*弧度转换为分*/double miao(double); /*弧度转换为秒*/#define PI 3.1415926void main (void){int k;printf("请选择执行正算或者反算，若执行正算，请输入1；若执行反算，请输入2。
\n");scanf("%d",&k);/*正算*/if(k==1){double bz,lz,az,S,bz2,lz2,az2,B1,L1,A1,B2,L2,A2,bx,by,lx,ly,ax,ay;int bx2,by2,lx2,ly2,ax2,ay2;doublee2,W1,sinu1,cosu1,sinA0,coto1,sin2o1,cos2o1,sin2o,cos2o,A,B,C,r,t,o0,o,g,si nu2,q;/*以度分秒顺序输入数据*/printf("请输入大地线起点纬度度分秒\n");scanf("%lf%lf%lf",&bx,&by,&bz);printf("请输入大地线起点经度度分秒\n");scanf("%lf%lf%lf",&lx,&ly,&lz);printf("请输入大地方位角度分秒\n");scanf("%lf%lf%lf",&ax,&ay,&az);printf("请输入大地线长度\n");scanf("%lf",&S);/*调用函数*/B1=hudu(bx,by,bz);L1=hudu(lx,ly,lz);A1=hudu(ax,ay,az);/*白塞尔大地主题解算*/e2=0.006693421622966;W1=sqrt(1-e2*sin(B1)*sin(B1));sinu1=sin(B1)*(sqrt(1-e2))/W1;cosu1=cos(B1)/W1;sinA0=cosu1*sin(A1);coto1=cosu1*cos(A1)/sinu1;sin2o1=2*coto1/(coto1*coto1+1);cos2o1=(coto1*coto1-1)/(coto1*coto1+1);A=6356863.020+(10718.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);B=(5354.469-8.978*(1-sinA0*sinA0))*(1-sinA0*sinA0);C=(2.238*(1-sinA0*sinA0))*(1-sinA0*sinA0)+0.006;r=691.46768-(0.58143-0.00144*(1-sinA0*sinA0))*(1-sinA0*sinA0);t=(0.2907-0.0010*(1-sinA0*sinA0))*(1-sinA0*sinA0);o0=(S-(B+C*cos2o1)*sin2o1)/A;sin2o=sin2o1*cos(2*o0)+cos2o1*sin(2*o0);cos2o=cos2o1*cos(2*o0)-sin2o1*sin(2*o0);o=o0+(B+5*C*cos2o)*sin2o/A;g=(r*o+t*(sin2o-sin2o1))*sinA0;/*求B2*/sinu2=sinu1*cos(o)+cosu1*cos(A1)*sin(o);B2=atan(sinu2/(sqrt(1-e2)*sqrt(1-sinu2*sinu2)));/*求L2*/q=atan(sin(A1)*sin(o)/(cosu1*cos(o)-sinu1*sin(o)*cos(A1)));/*判断q*/if(sin(A1)>0 && tan(q)>0)q=fabs(q);else if(sin(A1)>0 && tan(q)<0)q=PI-fabs(q);else if(sin(A1)<0 && tan(q)<0)q=-fabs(q);elseq=fabs(q)-PI;L2=L1+q-g/3600/180*PI;/*求A2*/A2=atan(cosu1*sin(A1)/(cosu1*cos(o)*cos(A1)-sinu1*sin(o))); /*判断A2*/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);elseA2=2*PI-fabs(A2);/*调用函数*/bx2=(int)(du(B2));by2=(int)(fen(B2));bz2=miao(B2);lx2=(int)(du(L2));ly2=(int)(fen(L2));lz2=miao(L2);ax2=(int)(du(A2));ay2=(int)(fen(A2));az2=miao(A2);printf("大地线终点纬度度分秒分别为：\n%d\n%d\n%lf\n",bx2,by2,bz2);printf("大地线终点经度度分秒分别为：\n%d\n%d\n%lf\n",lx2,ly2,lz2);printf("终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",ax2,ay2,az2);}/*反算*/else if(k==2){doublebz,lz,bz2,lz2,az,az2,B1,L1,B2,L2,S,A1,A2,bx,by,lx,ly,bx2,by2,lx2,ly2;int ax,ay,ax2,ay2;doublee2,W1,W2,sinu1,sinu2,cosu1,cosu2,L,a1,a2,b1,b2,g,g2,g0,r,p,q,sino,coso,o,si nA0,x,t1,t2,A,B,C,y;/*以度分秒顺序输入数据*/printf("请输入大地线起点纬度度分秒\n");scanf("%lf%lf%lf",&bx,&by,&bz);printf("请输入大地线起点经度度分秒\n"); scanf("%lf%lf%lf",&lx,&ly,&lz);printf("请输入大地线终点纬度度分秒\n"); scanf("%lf%lf%lf",&bx2,&by2,&bz2); printf("请输入大地线终点经度度分秒\n"); scanf("%lf%lf%lf",&lx2,&ly2,&lz2);/*调用函数*/B1=hudu(bx,by,bz);L1=hudu(lx,ly,lz);B2=hudu(bx2,by2,bz2);L2=hudu(lx2,ly2,lz2);/*白塞尔大地主题解算*/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;sinu2=sin(B2)*sqrt(1-e2)/W2;cosu1=cos(B1)/W1;cosu2=cos(B2)/W2;L=L2-L1;a1=sinu1*sinu2;a2=cosu1*cosu2;b1=cosu1*sinu2;b2=sinu1*cosu2;/*逐次趋近法求解A1*/g0=-662.904266/3600*PI/180; g=0;r=L;while(1){p=cosu2*sin(r);q=b1-b2*cos(r);A1=atan(p/q);/*判断A1*/if(p>0 && q>0)A1=fabs(A1);else if(p>0 && q<0)A1=PI-fabs(A1);else if(p<0 && q<0)A1=PI+fabs(A1);elseA1=2*PI-fabs(A1);sino=p*sin(A1)+q*cos(A1);coso=a1+a2*cos(r);o=atan(sino/coso);/*判断o*/if(coso>0)o=fabs(o);elseo=PI-fabs(o);sinA0=cosu1*sin(A1);x=2*a1-(1-sinA0*sinA0)*coso;t1=(33523299-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*0.0000000001;t2=(28189-94*(1-sinA0*sinA0))*0.0000000001;g2=(t1*o-t2*x*sino)*sinA0;/*检验循环次数*/printf("\ng2=%lf\ng0=%lf\n",g2,g0);if(g2<=g0)break;elser=L+g2;}/*求解S*/A=6356863.020+(10708.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);B=10708.938-17.956*(1-sinA0*sinA0);C=4.487;y=((1-sinA0*sinA0)*(1-sinA0*sinA0)-2*x*x)*coso; S=A*o+(B*x+C*y)*sino;/*求解A2*/A2=atan(cosu1*sin(r)/(b1*cos(r)-b2));/*判断A2*/if(p<0 && q<0)A2=fabs(A2);else if(p<0 && q>0)A2=PI-fabs(A2);else if(p>0 && q>0)A2=PI+fabs(A2);elseA2=2*PI-fabs(A2);/*调用函数*/ax=(int)(du(A1));ay=(int)(fen(A1));az=miao(A1);ax2=(int)(du(A2));ay2=(int)(fen(A2));az2=miao(A2);printf("起点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",ax,ay,az);printf("终点大地方位角度分秒分别为：\n%d\n%d\n%lf\n",ax2,ay2,az2);printf("大地线长度为：%lf\n",S);}/*数据错误*/elseprintf("数据错误，请重新输入\n");}/*度分秒转换为弧度*/double hudu(double a0,double b0,double c0){double A0;A0=(a0+b0/60+c0/3600)*PI/180;return A0;}/*弧度转换为度*/double du(double B0){double x0;x0=(int)(B0*180/PI);return x0;}/*弧度转换为分*/double fen(double C0){double _y,y0;_y=(int)(C0*180/PI);y0=(fabs)((int)((C0*180/PI-_y)*60));return y0;}/*弧度转换为秒*/double miao(double D0){double _z1,_z2,z0;_z1=(int)(D0*180/PI);_z2=(int)((D0*180/PI-_z1)*60);z0=(fabs)((double)(((D0*180/PI-_z1)*60-_z2)*60));return z0;}大地主题解算正算：大地主题解算反算：以上检验数据来自书本P151，P152。
大地主题解算方法综述
测绘科学 Sc ience o f Survey ing and M app ing
V o l132 N o14 Ju l1
大地主题解算方法综述
周振宇, 郭广礼, 贾新果
(中国矿业大学 环测学院, 江苏徐 州 221008)
【摘 要】大地主题解算是大地测量中的 重要问 题, 由于 椭球的 复杂性, 随 之产生 的解算 方法也 是多种 多样。本
和方位角差。
像这类方法 还有 很多, 就 不再一 一赘 述, 其 基本 原理
大致相同。
21 2 以白塞尔 大地投影为基础
可以近似认 为地球 椭球 的形状 与圆 球区 别不大。 在椭
球面上解算大地主 题问题 时可 借助于 球面 三角 学公式 简单
而严密的进 行。因此, 如 将椭 球面上 的大 地线长 度投 影到
在于: 解算精 度与距 离有 关, 距离 越长, 收敛 越慢, 因此
只适用于较短的距离。当公式 为四次项时 , 在中纬度 地区,
公式可用于 200km 以下的 大地主 题解算, 经纬 度计算 精度
可达 01 0001", 方位 角计 算精度 达 01001", 计算 边长 可精
确至厘 米。是 短 距离 解 算 的理 想 公式。 当距 离 小于 70km
椭球的过渡。
白塞尔 ( Besse l) 首先提出并解决了投影条件, 使这一解
法得以实现。
这类公式的特点是计算 公式展 开 e2 或 的幂级 数, 解算
精度与距离 长短无 关。因此 它既 适用于 短距 离解 算, 也适 用于长距离 解算 [ 5] 。其主 要缺 点在 于: 由 S 求 R、由 L 求 K, 或相反的运算, 需要进行迭代。同时还要预先 计算辅助
大地主题解算-C#
大地主题解算-正算-C#大地主题解算-正算-程序using System;using System.Collections.Generic;using ponentModel;using System.Data;using System.Drawing;using System.Linq;using System.Text;using System.Windows.Forms;namespace WindowsFormsApplication2{public partial class Form1 : Form{public Form1(){InitializeComponent();}private void button1_Click(object sender, EventArgs e) {int n,du,fen,miao;double B1, L1, A1, S, dS, dB1, dL1, M1, N1, W1;double E = 0.006694384999588, a = 6378140, C;B1 =( Convert.ToDouble(textBox1.Text) + Convert.ToDouble(textBox9.Text)/60 + Convert.ToDouble(textBox10.Text)/3600 )* Math.PI / 180;L1 = (Convert.ToDouble(textBox2.Text) + Convert.ToDouble(textBox11.Text)/60 + Convert.ToDouble(textBox12.Text)/3600) * Math.PI / 180; ;A1 = (Convert.ToDouble(textBox3.Text) + Convert.ToDouble(textBox13.Text)/60 + Convert.ToDouble(textBox14.Text)/3600) * Math.PI / 180; ;S = Convert.ToDouble(textBox5.Text);n = Convert.ToInt16(textBox4.Text);double[] L = new double[100000];double[] B = new double[100000];double[] dL = new double[100000];double[] dB = new double[100000];double[] A = new double[100000];double[] W = new double[100000];double[] N = new double[100000];double[] M = new double[100000];W1 = Math.Pow(1 - E * Math.Pow(Math.Sin(B1), 2), 0.5);N1 = a / W1;M1 = a * (1 - E) / (W1 * W1 * W1);C = N1 * Math.Sin(A1) * Math.Cos(B1);dS = S / n;dB1 = dS * Math.Cos(A1) / M1;dL1 = dS * Math.Sin(A1) / N1 / Math.Cos(B1);B[0] = B1;L[0] = L1;dB[0] = dB1;dL[0] = dL1;A[0] = A1;M[0] = M1;N[0] = N1;W[0] = W1;int i;for (i = 1; i < n; i++){dB[i - 1] = dS * Math.Cos(A[i - 1]) / M[i - 1];B[i] = B[i - 1] + dB[i - 1];W[i] = Math.Pow(1 - E * Math.Pow(Math.Sin(B[i]), 2), 0.5);N[i] = a / W[i];M[i] = a * (1 - E) / (W[i] * W[i] * W[i]);dL[i - 1] = dS * Math.Sin(A[i - 1]) / (N[i - 1] * Math.Cos(B[i - 1])); L[i] = L[i - 1] + dL[i - 1];A[i] = Math.Asin(C / (Math.Cos(B[i]) * N[i]));}du = Convert.ToInt16(Math.Floor(B[i - 1] * 180 / Math.PI));//dufen = Convert.ToInt16(Math.Floor((B[i - 1] * 180 / Math.PI - Math.Floor(B[i - 1] * 180 / Math.PI)) * 60)); //fenmiao = Convert.ToInt16(Convert.ToInt16(((B[i - 1] * 180 / Math.PI - Math.Floor(B[i - 1] * 180 / Math.PI)) * 60 - Math.Floor((B[i - 1] * 180 / Math.PI - Math.Floor(B[i - 1] * 180 / Math.PI)) * 60)) * 60));//秒?textBox6.Text = Convert.ToString(du) + "度è" + Convert.ToString(fen) + "分?" + Convert.ToString(miao) + "秒?";du = Convert.ToInt16(Math.Floor(L[i - 1] * 180 / Math.PI));//dufen = Convert.ToInt16(Math.Floor((L[i - 1] * 180 / Math.PI - Math.Floor(L[i - 1] * 180 / Math.PI)) * 60)); //fenmiao = Convert.ToInt16(Convert.ToInt16(((L[i - 1] * 180 / Math.PI - Math.Floor(L[i - 1] * 180 / Math.PI)) * 60 - Math.Floor((L[i - 1] * 180 / Math.PI - Math.Floor(L[i - 1] * 180 / Math.PI)) * 60)) * 60));//秒?textBox7.Text = Convert.ToString(du) + "度è" + Convert.ToString(fen) + "分?" + Convert.ToString(miao) + "秒?";du = Convert.ToInt16(Math.Floor((180 + A[i - 1] * 180 / Math.PI)));//dufen = Convert.ToInt16(Math.Floor(((180 + A[i - 1] * 180 / Math.PI) - Math.Floor((180 + A[i - 1] * 180 / Math.PI))) * 60)); //fenmiao = Convert.ToInt16(Convert.ToInt16((((180 + A[i - 1] * 180 / Math.PI) - Math.Floor((180 + A[i - 1] * 180 / Math.PI))) * 60 - Math.Floor(((180 + A[i - 1] * 180 / Math.PI) - Math.Floor((180 + A[i - 1] * 180 / Math.PI))) * 60)) * 60));//秒?textBox8.Text = Convert.ToString(du) + "度è" + Convert.ToString(fen) + "分?" + Convert.ToString(miao) + "秒?";}private void Form1_Load(object sender, EventArgs e){} }}。
白塞尔大地主题解算(正算和反算)
大地测量实验报告实验名称：白塞尔大地主题解算（正算和反算）实验目的：1.通过编写白塞尔大地主题电算程序进一步掌握白塞尔法解算大地主题的基本思想。
2.熟练掌握将椭球面上的大地元素按照白塞尔投影条件投影到辅助球面上，继而在球面上进行大地主题解算，最后再将球面上的计算结果换算到椭球面上的基本方法和步骤。
3.学会掌握计算机编程的基本能力。
实验环境：Microsoft Visual C++注意事项：1.在编写程序的过程当中要注意代码的前后统一和重复。
2.注意数值类型的转换和度分秒的换算。
实验步骤：正算：1.计算起点的规划纬度2.计算辅助函数值3.计算系数A,B,C及d,e.4计算球面长度5.计算经度差改正数6.计算终点大地坐标及大地方位角。
反算:1.进行计算辅助函数值2.用逐次趋近法同时计算起点大地方位角、球面长度及经差。
3.计算系数A,B,C及大地线长度S.4.计算反方位角及确定符号。
程序源代码：正算：#include<stdio.h>#include<math.h>#define ee 0.006694384999588#define I 3.141592653double F(double,double,double);void main(void){double A1,B1,L1,S,A2,B2,L2; double x1,x2,x3,y1,y2,y3,z1,z2,z3; double W1,sinu1,sinu2,cosu1,sinA0;doublecota1,cos2a1,sin2a1,cosA0A0;double A,B,C,d,e,a0,a1,m;double n,a,Q,R;printf("请输入数据B1= "); scanf("%lf %lf %lf",&x1,&x2,&x3);B1=F(x1,x2,x3);printf("请输入数据L1= "); scanf("%lf %lf %lf",&y1,&y2,&y3);L1=F(y1,y2,y3);printf("请输入A1= ");scanf("%lf %lf %lf",&z1,&z2,&z3);A1=F(z1,z2,z3);printf("请输入S= ");scanf("%lf",&S);printf("B1=%f\n",B1);printf("L1=%f\n",L1);printf("A1=%f\n",A1);printf("S=%f\n",S);/*计算起点的规划纬度*/W1=sqrt(1-ee*sin(B1)*sin(B1));sinu1=sin(B1)*sqrt(1-ee)/W1;cosu1=cos(B1)/W1;printf("W1=%f\n",W1);printf("sinu1=%f\n",sinu1);printf("cosu1=%f\n",cosu1);/*计算辅助函数值*/sinA0=cosu1*sin(A1);cota1=cosu1*cos(A1)/sinu1;sin2a1=2*cota1/(cota1*cota1+1);cos2a1=(cota1*cota1-1)/(cota1*cota1+1);printf("sinA0=%f\n",sinA0);printf("cota1=%f\n",cota1);printf("sin2a1=%f\n",sin2a1);printf("cos2a1=%f\n",cos2a1);/*计算系数ABC及de*/cosA0A0=1-sinA0*sinA0;A=6356755.288+(10710.341-(13.534*cosA0A0))*cosA0A0;B=(5355.171-9.023*cosA0A0)*cosA0A0;C=(2.256*(cosA0A0))*cosA0A0+0.006;d=691.46768-(0.58143-0.00144*cosA0A0)*cosA0A0;e=(0.2907-cosA0A0*0.0010)*cosA0A0;printf("cosA0A0=%f\n",cosA0A0);printf("A=%f\n",A);printf("B=%f\n",B);printf("C=%f\n",C);printf("d=%f\n",d);printf("e=%f\n",e);/*计算球面长度*/a0=(S-(B+C*cos2a1)*sin2a1)/A;m=sin2a1*cos(2*a0)+cos2a1*sin(2*a0);n=(cos2a1)*(cos(2*a0))-(sin2a1)*(sin(2*a0));a=a0+((B+5*C*n))*m/A;printf("a0=%f\n",a0);printf("m=%f\n",m);printf("n=%f\n",n);printf("a=%f\n",a);/*计算经度差改正数*/Q=(d*a+(e*(m-sin2a1))/3600/180*I)*sinA0;printf("Q=%f\n",Q);/*计算终点大地坐标及大地方位角*/sinu2=sinu1*cos(a)+cosu1*cos(A1)*sin(a);B2=180*atan(sinu2/((sqrt(1-ee))*(sq rt(1-sinu2*sinu2))))/I;R=180*atan(sin(A1)*sin(a)/(cosu1*co s(a)-sinu1*sin(a)*cos(A1)))/I;printf("sinu2=%f\n",sinu2);printf("B2=%f\n",B2);printf("R=%f\n",R*180/I);/*确定R的值*/if(sin(A1)>0 && tan(R)>0)R=abs(R);else if(sin(A1)>0 && tan(R)<0)R=I-abs(R);else if(sin(A1)<0 && tan(R)<0)R=-abs(R);elseR=abs(R)-I;/*确定L2A2的值*/L2=(L1*180/I+R-(Q/206265*180/I));A2=atan(cosu1*sin(A1)/(cosu1*cos(a) *cos(A1)-sinu1*sin(a)));if(sin(A1)<0&&tan(A2)>0)A2=(fabs(A2))*180/I;else if(sin(A1)<0&&tan(A2)<0)A2=(I-fabs(A2))*180/I;else if(sin(A1)>0&&tan(A2)>0)A2=(I+fabs(A2))*180/I;elseA2=(2*I-fabs(A2))*180/I;printf("A2=%3f\n B2=%3f\nL2=%3f\n",A2,B2,L2);}double F(double a2,double b2,doublec2){double d2;d2=(double)(a2+1.0*b2/60+1.0*c2/3600);d2=(d2/180)*I;return (d2);}注：A1，B1，L1，S分别为大地线起点的大地方位角、纬度、经度、大地线长;B2,L2,A2为大地线终点纬度、经度及方位角。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
大地主题解算(正算)代码与白塞尔大地主题解算大地主题解算(正算)代码：根据经纬度和方向角以及距离计算另外一点坐标新建模块->拷贝下面的大地主题(正算)代码，调用方法示例：起点经度：116.235(度)终点纬度：37.435(度)方向角：50（度）长度：500（米）终点经纬度（"经度,纬度"）=Computation(37.435,116.235,50,500)Const Pi = 3.14Private a, b, c, alpha, e, e2, W, V As Double'a 长轴半径'b 短轴'c 极曲率半径'alpha 扁率'e 第一偏心率'e2 第二偏心率'W 第一基本纬度函数'V 第二基本纬度函数Private B1, L1, B2, L2 As Double'B1 点1的纬度'L1 点1的经度'B2 点1的纬度'L2 点2的经度Private S As Double '''''大地线长度Private A1, A2 As Double'A1 点1到点2的方位角'A2 点2到点1的方位角Function Computation(STARTLAT, STARTLONG, ANGLE1, DISTANCE As Double) As StringB1 = STARTLATL1 = STARTLONGA1 = ANGLE1S = DISTANCEa = 6378245b = 6356752.3142c = a ^ 2 / balpha = (a - b) / ae = Sqr(a ^ 2 - b ^ 2) / ae2 = Sqr(a ^ 2 - b ^ 2) / bB1 = rad(B1)L1 = rad(L1)A1 = rad(A1)W = Sqr(1 - e ^ 2 * (Sin(B1) ^ 2))V = W * (a / b)Dim W1 As DoubleE1 = e ''''第一偏心率'// 计算起点的归化纬度W1 = W ''Sqr(1 - e1 * e1 * Sin(B1 ) * Sin(B1 )) sinu1 = Sin(B1) * Sqr(1 - E1 * E1) / W1cosu1 = Cos(B1) / W1'// 计算辅助函数值sinA0 = cosu1 * Sin(A1)cotq1 = cosu1 * Cos(A1)sin2q1 = 2 * cotq1 / (cotq1 ^ 2 + 1)cos2q1 = (cotq1 ^ 2 - 1) / (cotq1 ^ 2 + 1)'// 计算系数AA,BB,CC及AAlpha, BBeta的值。
cos2A0 = 1 - sinA0 ^ 2e2 = Sqr(a ^ 2 - b ^ 2) / bk2 = e2 * e2 * cos2A0Dim aa, BB, CC, EE22, AAlpha, BBeta As Doubleaa = b * (1 + k2 / 4 - 3 * k2 * k2 / 64 + 5 * k2 * k2 * k2 / 256)BB = b * (k2 / 8 - k2 * k2 / 32 + 15 * k2 * k2 * k2 / 1024)CC = b * (k2 * k2 / 128 - 3 * k2 * k2 * k2 / 512)e2 = E1 * E1AAlpha = (e2 / 2 + e2 * e2 / 8 + e2 * e2 * e2 / 16) - (e2 * e2 / 16 + e2 * e2 * e2 / 16) * cos2A0 + (3 * e2 * e2 * e2 / 128) * cos2A0 * cos2A0BBeta = (e2 * e2 / 32 + e2 * e2 * e2 / 32) * cos2A0 - (e2 * e2 * e2 / 64) * cos2A0 * cos2A0'// 计算球面长度q0 = (S - (BB + CC * cos2q1) * sin2q1) / aasin2q1q0 = sin2q1 * Cos(2 * q0) + cos2q1 * Sin(2 * q0)cos2q1q0 = cos2q1 * Cos(2 * q0) - sin2q1 * Sin(2 * q0)q = q0 + (BB + 5 * CC * cos2q1q0) * sin2q1q0 / aa'// 计算经度差改正数theta = (AAlpha * q + beta * (sin2q1q0 - sin2q1)) * sinA0'// 计算终点大地坐标及大地方位角sinu2 = sinu1 * Cos(q) + cosu1 * Cos(A1) * Sin(q)B2 = Atn(sinu2 / (Sqr(1 - E1 * E1) * Sqr(1 - sinu2 * sinu2))) * 180 / Pilamuda = Atn(Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1))) * 180 / PiIf (Sin(A1) > 0) ThenIf (Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1)) > 0) Thenlamuda = Abs(lamuda)Elselamuda = 180 -Abs(lamuda)End IfElseIf (Sin(A1) * Sin(q) / (cosu1 * Cos(q) - sinu1 * Sin(q) * Cos(A1)) > 0) Thenlamuda = Abs(lamuda) - 180Elselamuda = -Abs(lamuda)End IfEnd IfL2 = L1 * 180 / Pi + lamuda - theta * 180 / PiA2 = Atn(cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q))) * 180 / PiIf (Sin(A1) > 0) ThenIf (cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q)) > 0) ThenA2 = 180 + Abs(A2)ElseA2 = 360 - Abs(A2)End IfElseIf (cosu1 * Sin(A1) / (cosu1 * Cos(q) * Cos(A1) - sinu1 * Sin(q)) > 0) ThenA2 = Abs(A2)ElseA2 = 180 - Abs(A2)End IfEnd IfComputation = L2 & "," & B2End FunctionPrivate Function rad(ByVal angle_d As Double) As Doublerad = angle_d * Pi / 180End Function白塞尔大地主题解算一、正算代码：#include<stdio.h>#include<math.h>#define ee 0.2966#define I 3.141592653double F(double,double,double);void main(void){double A1,B1,L1,S,A2,B2,L2;double x1,x2,x3,y1,y2,y3,z1,z2,z3;double W1,sinu1,sinu2,cosu1,sinA0;double cota1,cos2a1,sin2a1,cosA0A0;double A,B,C,d,e,a0,a1,m;double n,a,Q,R;printf("请输入数据B1= ");scanf("%lf %lf %lf",&x1,&x2,&x3);B1=F(x1,x2,x3);printf("请输入数据L1= ");scanf("%lf %lf %lf",&y1,&y2,&y3);L1=F(y1,y2,y3);printf("请输入A1= ");scanf("%lf %lf %lf",&z1,&z2,&z3);A1=F(z1,z2,z3);printf("请输入S= ");scanf("%lf",&S);printf("B1=%.9f\n",B1);printf("L1=%.9f\n",L1);printf("A1=%.9f\n",A1);printf("S=%f\n",S);/*计算起点的规划纬度*/W1=sqrt(1-ee*sin(B1)*sin(B1));sinu1=sin(B1)*sqrt(1-ee)/W1;cosu1=cos(B1)/W1;printf("W1=%.9f\n",W1);printf("sinu1=%.9f\n",sinu1);printf("cosu1=%.9f\n",cosu1);/*计算辅助函数值*/sinA0=cosu1*sin(A1);cota1=cosu1*cos(A1)/sinu1;sin2a1=2*cota1/(cota1*cota1+1);cos2a1=(cota1*cota1-1)/(cota1*cota1+1);printf("sinA0=%.9f\n",sinA0);printf("cota1=%.9f\n",cota1);printf("sin2a1=%.9f\n",sin2a1);printf("cos2a1=%.9f\n",cos2a1);/*计算系数ABC及de*/cosA0A0=1-sinA0*sinA0;A=6356755.288+(10710.341-(13.534*cosA0A0))*cosA0A0;B=(5355.171-9.023*cosA0A0)*cosA0A0;C=(2.256*(cosA0A0))*cosA0A0+0.006;d=691.46768-(0.58143-0.00144*cosA0A0)*cosA0A0;e=(0.2907-cosA0A0*0.0010)*cosA0A0;printf("cosA0A0=%.9f\n",cosA0A0);printf("A=%.3f\n",A);printf("B=%.6f\n",B);printf("C=%.9f\n",C);printf("d=%.7f\n",d);printf("e=%.9f\n",e);/*计算球面长度*/a0=(S-(B+C*cos2a1)*sin2a1)/A;m=sin2a1*cos(2*a0)+cos2a1*sin(2*a0);n=(cos2a1)*(cos(2*a0))-(sin2a1)*(sin(2*a0));a=a0+((B+5*C*n))*m/A;printf("a0=%.9f\n",a0);printf("m=%.9f\n",m);printf("n=%.9f\n",n);printf("a=%.9f\n",a);/*计算经度差改正数*/Q=(d*a+(e*(m-sin2a1))/3600/180*I)*sinA0;printf("Q=%.7f\n",Q);/*计算终点大地坐标及大地方位角*/sinu2=sinu1*cos(a)+cosu1*cos(A1)*sin(a);B2=180*atan(sinu2/((sqrt(1-ee))*(sqrt(1-sinu2*sinu2))))/I;R=180*atan(sin(A1)*sin(a)/(cosu1*cos(a)-sinu1*sin(a)*cos(A1)))/I;printf("sinu2=%.9f\n",sinu2);printf("B2=%f\n",B2);printf("R=%f\n",R*180/I);/*确定R的值*/if(sin(A1)>0 && tan(R)>0)R=abs(R);else if(sin(A1)>0 && tan(R)<0)R=I-abs(R);else if(sin(A1)<0 && tan(R)<0)R=-abs(R);elseR=abs(R)-I;/*确定L2A2的值*/L2=(L1*180/I+R-(Q/206265*180/I));A2=atan(cosu1*sin(A1)/(cosu1*cos(a)*cos(A1)-sinu1*sin(a)));if(sin(A1)<0&&tan(A2)>0)A2=(fabs(A2))*180/I;else if(sin(A1)<0&&tan(A2)<0)A2=(I-fabs(A2))*180/I;else if(sin(A1)>0&&tan(A2)>0)A2=(I+fabs(A2))*180/I;elseA2=(2*I-fabs(A2))*180/I;printf("A2=%3f\n B2=%3f\n L2=%3f\n",A2,B2,L2);}double F(double a2,double b2,double c2){double d2;d2=(double)(a2+1.0*b2/60+1.0*c2/3600);d2=(d2/180)*I;return (d2);}运行结果：二、反算代码：#include<stdio.h>#include<math.h>#define ee 0.2966#define I 3.141592653double F(double,double,double);void main(void){double B1,B2,L1,L2,A1,A2,S,Y;double W1,W2,L,Q,R,A,B,C;double x,y,z,p,q;double x1,x2,x3,y1,y2,y3,z1,z2,z3,w1,w2,w3;double a1,a2,b1,b2,m,n;double sinp,cosp,sinu1,sinu2,cosu1,cosu2,sinA0,cosA0;Q=0;q=0;printf("请输入起始数据B1= ");scanf("%lf %lf %lf",&x1,&x2,&x3);B1=F(x1,x2,x3);printf("请输入起始数据L1= ");scanf("%lf %lf %lf",&y1,&y2,&y3);L1=F(y1,y2,y3);printf("请输入起始数B2= ");scanf("%lf %lf %lf",&z1,&z2,&z3);B2=F(z1,z2,z3);printf("请输入起始数据L2= ");scanf("%lf %lf %lf",&w1,&w2,&w3);L2=F(w1,w2,w3);printf("B1=%f\n",B1);printf("L1=%.9f\n",L1);printf("B2=%.9f\n",B2);printf("L2=%.9f\n",L2);/*辅助计算*/W1=sqrt(1-ee*sin(B1)*sin(B1));W2=sqrt(1-ee*sin(B2)*sin(B2));sinu1=sin(B1)*sqrt(1-ee)/W1;sinu2=sin(B2)*(sqrt(1-ee))/W2;cosu1=cos(B1)/W1;cosu2=cos(B2)/W2;L=L2-L1;a1=sinu1*sinu2;a2=cosu1*cosu2;b1=cosu1*sinu2;b2=sinu1*cosu2;printf("W1=%.9f\n",W1);printf("W2=%.9f\n",W2);printf("sinu1=%.9f\n",sinu1);printf("sinu2=%.9f\n",sinu2);printf("cosu1=%.9f\n",cosu2);printf("L=%.9f\n",L);printf("a1=%.9f\n",a1);printf("a2=%.9f\n",a2);printf("b1=%.9f\n",b1);printf("b2=%.9f\n",b2);/*用逐次趋近法同时计算起点大地方位角、球面长度及经差R*/ do{R=L+Q;x=cosu2*sin(R);y=b1-b2*cos(R);A1=atan(x/y);if(x>0&&y>0)A1=fabs(A1);else if(x>0&&y<0)A1=I-fabs(A1);else if(x<0&&y<0)A1=I+fabs(A1);elseA1=2*I-fabs(A1);sinp=x*sin(A1)+y*cos(A1);cosp=a1+a2*cos(R);p=atan(sinp/cosp);if(cosp>0)p=fabs(p);sinA0=cosu1*sin(A1);cosA0=sqrt(1-sinA0*sinA0);p=I-fabs(p);z=2*a1-cosA0*cosA0*cos(p);m=(33523299-(28189-70*cosA0*cosA0))*(1e-10);n=(28189-94*cosA0*cosA0)*(1e-10);Q=q;q=(m*p-m*z*sin(p))*sinA0;}while(fabs(206265*(q-Q))>(1e-4));printf("Q=%f\n",Q);/*计算系数ABC及大地线长度*/A=6356863.020+(10708.949-13.474*cosA0*cosA0)*cosA0*cosA0;B=10708.938-17.956*cosA0*cosA0;C=4.487;Y=(cosA0*cosA0*cosA0*cosA0-2*z*z)*cos(p);S=A*p+(B*z+C*Y)*sinp;/*计算反方位角*/A2=atan((cosu1*sin(R))/(b1*cos(R)-b2));A2=A2*180/I;A1=A1*180/I;if(A1>0)A2=abs(A2);elseA2=-abs(A2);printf("A1=%3f\n A2=%3f\n S=%3f\n",A1,A2,S);}double F(double a,double b,double c){double d;if(a<0)d=(double)(a-1.0*b/60-1.0*c/3600);elsed=(double)(a+1.0*b/60+1.0*c/3600);d=(d/180)*I;return d; }运行结果：。