8高斯投影正反算 C#代码

高斯投影正反算程序设计一．程序设计流程本程序(de)设计思路如下：（1）,程序采用VS08版本作为开发平台,并采用C语言作为开发语言,设计为WindowsForm窗体程序形式.（2）,程序主要(de)算法来自于教材.但是本程序为了更加实用,添加了更多(de)解算基准,包括：WGS-84,国际椭球1975,克氏椭球,和2000国家大地坐标系.（3）,程序为了更方便(de)读取数据和输出数据,故需要自己定义了固定(de)数据输入格式和数据输出格式或形式,请老师注意查看.二．代码using System;using Systusing SystemponentModel;using System.Data;using System.Drawing;using System.Text;namespace Gauss{public partial class Form1 : Form{//大地坐标//Geodetic Coordinatepublic struct CRDGEODETIC{public double dLongitude;public double dLatitude;public double dHeight;}//笛卡尔坐标//Cartesian Coordinatepublic struct CRDCARTESIAN{public double x;public double y;public double z;}public Form1(){InitializeComponent();}private void button1_Click(object sender, EventArgs e) {double ee = 0;double a = 0;string tt;try{}catch{MessageBox.Show("Gauss Inverse: Choose datum error");return;}if (ttpareTo("克氏椭球")==0){a = 6378245.00;}if (ttpareTo("WGS-84") == 0){a = 6378.00;}if (ttpareTo("1975国际椭球") == 0){a = 6378140.00;ee =}if (ttpareTo("2000国家大地坐标系") == 0){a = 6378.0;}const double pai = 3.1415926;double b = Math.Sqrt(a a (1 - ee ee));double c = a a / b;double epp = Math.Sqrt((a a - b b) / b / b);CRDGEODETIC pcrdGeo;CRDCARTESIAN pcrdCar;double midlong;//求纬度string[] temp;double[] tempradius = new double[3];for (int i = 0; i < 3; i++){tempradius[i] = Convert.ToDouble(temp[i]);}pcrdGeo.dLatitude = tempradius[0] / 180.0 pai + tempradius[1] / 180.0 / 60.0 pai + tempradius[2] / 180 / 60.0 / 60 pai;//求经度for (int i = 0; i < 3; i++){tempradius[i] = Convert.ToDouble(temp[i]);}pcrdGeo.dLongitude = tempradius[0] / 180.0 pai + tempradius[1] / 180.0 / 60.0 pai + tempradius[2] / 180 / 60.0 / 60 pai;int deglon = Convert.ToInt32(pcrdGeo.dLongitude180 / pai);//求中央经度int num; //带号midlong = 0; //默认值,需要制定分带try{}catch{MessageBox.Show("Choose 3/6 error");return;}if (ttpareTo("3度带") == 0){num = Convert.ToInt32(deglon / 6 + 1);midlong = (6 num - 3) / 180.0 pai;}if (ttpareTo("6度带") == 0){num = Convert.ToInt32((deglon + 1.5) / 3);midlong = num 3 pai / 180;}double lp=pcrdGeo.dLongitude - midlong;double N = c / Math.Sqrt(1 + epp eppMath.Cos(pcrdGeo.dLatitude) Math.Cos(pcrdGeo.dLatitude));double M = c / Math.Pow(1 + epp eppMath.Cos(pcrdGeo.dLatitude) Math.Cos(pcrdGeo.dLatitude), 1.5);double ita = epp Math.Cos(pcrdGeo.dLatitude);double t = Math.Tan(pcrdGeo.dLatitude);double Nscnb = N Math.Sin(pcrdGeo.dLatitude) Math.Cos(pcrdGeo.dLatitude);double Ncosb = N Math.Cos(pcrdGeo.dLatitude);double cosb = Math.Cos(pcrdGeo.dLatitude);double X;double m0, m2, m4, m6, m8;double a0, a2, a4, a6, a8;m0 = a (1 - ee ee);m2 = 3.0 / 2.0 m0 ee ee;m4 = 5.0 / 4.0 ee ee m2;m6 = 7.0 / 6.0 ee ee m4;m8 = 9.0 / 8.0 ee ee m6;a0 = m0 + m2 / 2.0 + 3.0 / 8.0 m4 + 5.0 / 16.0 m6 + 35.0 / 128.0 m8;a2 = m2 / 2 + m4 / 2 + 15.0 / 32.0 m6 + 7.0 / 16.0 m8;a4 = m4 / 8.0 + 3.0 / 16.0 m6 + 7.0 / 32.0 m8;a6 = m6 / 32.0 + m8 / 16.0;a8 = m8 / 128.0;double B = pcrdGeo.dLatitude;double sb = Math.Sin(B);double cb = Math.Cos(B);double s2b = sb cb 2;double s4b = s2b (1 - 2 sb sb) 2;double s6b = s2b Math.Sqrt(1 - s4b s4b) + s4b Math.Sqrt(1 - s2b s2b);X = a0 B - a2 / 2.0 s2b + a4 s4b / 4.0 - a6 / 6.0 s6b;pcrdCar.x = Nscnb lp lp / 2.0 + Nscnb cosb cosb Math.Pow(lp, 4) (5 - t t + 9 ita ita + 4 Math.Pow(ita, 4)) / 24.0 + Nscnb Math.Pow(cosb, 4) Math.Pow(lp, 6) (61 - 58 t t + Math.Pow(t, 4)) / 720.0 + X;pcrdCar.y = Ncosb Math.Pow(lp, 1) + Ncosb cosb cosb (1 - t t + ita ita) / 6.0 Math.Pow(lp, 3) + NcosbMath.Pow(lp, 5) Math.Pow(cosb, 4) (5 - 18 t t+ Math.Pow(t, 4) + 14 ita ita - 58 ita ita t t) / 120.0 ;if (pcrdCar.y < 0)pcrdCar.y += 500000;richTextBox1.Text = "Results:\nX:\t" +Convert.ToString(pcrdCar.x) +"\nY:\t"+ Convert.ToString(pcrdCar.y);}private void button2_Click(object sender, EventArgs e) {double ee = 0;double a = 0;string tt;int num; //带号string ytext; //利用y值求带号和中央经线try{}catch{MessageBox.Show("Gauss Inverse: Choose datumerror");return;}if (ttpareTo("克氏椭球") == 0){a = 6378245.00;}if (ttpareTo("WGS-84") == 0){a = 6378.00;}if (ttpareTo("1975国际椭球") == 0){a = 6378140.00;}if (ttpareTo("2000国家大地坐标系") == 0){a = 6378.0;}double b = Math.Sqrt(a a (1 - ee ee));double c = a a / b;double epp = Math.Sqrt((a a - b b) / b / b); CRDGEODETIC pcrdGeo;CRDCARTESIAN pcrdCar;double midlong = 0;//求X,Y和带号pcrdCar.x = Convert.ToDouble(textBox4.Text); ytext = textBox5.Text;string temp = ytext.Substring(0, 2);num = Convert.ToInt32(temp);ytext = ytext.Remove(0, 2);pcrdCar.y = Convert.ToDouble(ytext) - 500000; try{}catch{MessageBox.Show("Choose 3/6 error");return;}if (ttpareTo("3度带") == 0){midlong = num 3 pai / 180;}if (ttpareTo("6度带") == 0){midlong = (6 num - 3) pai / 180;}b = Math.Sqrt(a a (1 - ee ee));c = a a / b;epp = Math.Sqrt(a a - b b) / b;double m0, m2, m4, m6, m8;double a0, a2, a4, a6, a8;m0 = a (1 - ee ee);m2 = 3.0 / 2.0 m0 ee ee;m4 = 5.0 / 4.0 ee ee m2;m6 = 7.0 / 6.0 ee ee m4;m8 = 9.0 / 8.0 ee ee m6;a0 = m0 + m2 / 2.0 + 3.0 / 8.0 m4 + 5.0 / 16.0 m6 + 35.0 / 128.0 m8;a2 = m2 / 2 + m4 / 2 + 15.0 / 32.0 m6 + 7.0 / 16.0 m8;a4 = m4 / 8.0 + 3.0 / 16.0 m6 + 7.0 / 32.0 m8;a6 = m6 / 32.0 + m8 / 16.0;a8 = m8 / 128.0;double Bf, B;Bf = pcrdCar.x / a0;B = 0.0;while (Math.Abs(Bf - B) > 1E-10){B = Bf;double sb = Math.Sin(B);double cb = Math.Cos(B);double s2b = sb cb 2;double s4b = s2b (1 - 2 sb sb) 2;double s6b = s2b Math.Sqrt(1 - s4b s4b) + s4b Math.Sqrt(1 - s2b s2b);Bf = (pcrdCar.x - (-a2 / 2.0 s2b + a4 / 4.0 s4b - a6 / 6.0 s6b)) / a0;}double itaf, tf, Vf, Nf;itaf = epp Math.Cos(Bf);tf = Math.Tan(Bf);Vf = Math.Sqrt(1 + epp epp Math.Cos(Bf)Math.Cos(Bf));Nf = c / Vf;double ynf = pcrdCar.y / Nf;pcrdGeo.dLatitude = Bf - 1.0 / 2.0 Vf Vf tf (ynf ynf - 1.0 / 12.0 Math.Pow(ynf, 4) (5 + 3 tf tf + itaf itaf - 9 Math.Pow(itaf tf, 2)) +1.0 / 360.0 (61 + 90 tf tf + 45 Math.Pow(tf, 4)) Math.Pow(ynf, 6));pcrdGeo.dLongitude = (ynf / Math.Cos(Bf) - (1 + 2 tf tf + itaf itaf) Math.Pow(ynf, 3) / 6.0 / Math.Cos(Bf) +(5 + 28 tf tf + 24 Math.Pow(tf, 4) + 6 itaf itaf + 8 Math.Pow(itaf tf, 2)) Math.Pow(ynf, 5) / 120.0 /Math.Cos(Bf));pcrdGeo.dLongitude = pcrdGeo.dLongitude + midlong;//pcrdGeo.dLatitude = pcrdGeo.dLatitude;richTextBox2.Text = "Results:\nLatitude: " + Convert.ToString(pcrdGeo.dLatitude) + "\nLongtitude: " +Convert.ToString(pcrdGeo.dLongitude);}private void label13_Click(object sender, EventArgs e) {}}}三．程序运行结果分析通过选取书上(de)具体实例进行测试,本程序(de)精度大体满足要求,一般正算(de)精度在0.01米和0.001米之间,反算(de)精度在0.0001秒左右.以下是程序运行(de)截图.。

高斯投影正反算学院：资源与环境工程工程学院专业：测绘工程 学号：X51414012：超一、高斯投影概述想象有一个椭圆柱面横套在地球椭球体外面，并与某一条子午线相切，椭圆柱的中心轴通过椭球体的中心，然后用一定投影方法，将中央子午线两侧各一定经差围的地区投影到椭圆柱面上，再将此柱面展开即成为投影面。

高斯投影由于是正形投影，故保证了投影的角度不变性，图形的相似性以及在某点各方向上长度比的同一性。

由于采用了同样法则的分带投影，这即限制了长度变形，又保证了在不同投影带中采用相同的简便公式和数表进行变形引起的各项改正的计算，并且带与带间的互相换算也能用相同的公式和方法进行。

高斯投影的这些优点必将使它得到广泛的推广和具有国际意义。

二、高斯投影坐标正算公式1.高斯投影必须满足以下三个条件 1）中央子午线投影后为直线 2）中央子午线投影后长度不变 3）投影具有正形性质，即正形投影条件2.高斯正算公式推导1）由第一个条件可知，由于地球椭球体是一个旋转椭球体，所以高斯投影必然有这样一个性质，即中央子午线东西两侧的投影必然对称于中央子午线。

2）由于高斯投影是换带投影，在每带经差ｌ是不大的，ｌρ是一个微小量，所以可以将X=X （ｌ，q ），Y=Y （l ，q ）展开为经差为l 的幂级数，它可写成如下的形式X=m 0+m 2l 2+m 4l 4+…Y=m 1l+m 3l 2+m 5l 5+…式中m 0,m1,m2,…是待定系数，他们都是纬度B 的函数。

3）由第三个条件：∂y ∂l =∂x ∂q 和∂x ∂l =-∂y∂q ，将上式分别对l 和q 求偏导2340123423401234...........x m m l m l m l m l y n n l n l n l n l =+++++=+++++可得到下式0312123403121234111,,,, 234111,,,,234dm dm dm dm n n n n dq dq dq dq dn dn dn dn m m m m dq dq dq dq ⎧====⎪⎪⎨⎪=-=-=-=-⎪⎩经过计算可以得出232244524632235242225sin cos sin cos (594)224 sin cos (6158)720cos cos (1)6cos (5181458)120N N x X B B l B B t l NB B t t l Ny N B l B t l NB t t t l ηηηηη=+⋅+-+++-+=⋅+-++-++-三、高斯投影坐标反算公式推导1.思路：级数展开，应用高斯投影三个条件，待定系数法求解。

§8.3高斯投影坐标正反算公式任何一种投影①坐标对应关系是最主要的；②如果是正形投影，除了满足正形投影的条件外（C-R 偏微分方程），还有它本身的特殊条件。

8.3.1高斯投影坐标正算公式： B,l ⇒ x,y高斯投影必须满足以下三个条件：①中央子午线投影后为直线；②中央子午线投影后长度不变；③投影具有正形性质，即正形投影条件。

由第一条件知中央子午线东西两侧的投影必然对称于中央子午线，即(8-10)式中，x 为l 的偶函数，y 为l 的奇函数；0330'≤l ，即20/1/≈''''ρl ，如展开为l 的级数，收敛。

+++=++++=553316644220l m l m l m y l m l m l m m x （8-33）式中 ,,10m m 是待定系数，它们都是纬度B 的函数。

由第三个条件知：qy l x l y q x ∂∂-=∂∂∂∂=∂∂, (8-33)式分别对l 和q 求偏导数并代入上式----=++++++=+++5533156342442204523164253l dqdm l dq dm l dq dm l m l m l m l dqdm l dq dm dq dm l m l m m (8-34) 上两式两边相等，其必要充分条件是同次幂l 前的系数应相等，即dq dm m dqdm m dqdm m 2312013121⋅=⋅-==(8-35)(8-35)是一种递推公式，只要确定了0m 就可依次确定其余各系数。

由第二条件知:位于中央子午线上的点，投影后的纵坐标x 应等于投影前从赤道量至该点的子午线弧长X ，即(8-33)式第一式中，当0=l时有：0m X x == (8-36) 顾及(对于中央子午线)B V Mr M B N dq dB M dBdXcos cos 2==== 得：B V cB N r dq dB dB dX dq dX dq dm m cos cos 01===⋅===(8-37,38)B B Ndq dB dB dm dq dm m cos sin 22121112=⋅-=⋅-= (8-39)依次求得6543,,,m m m m 并代入(8-33)式，得到高斯投影正算公式6425644223422)5861(cos sin 720)495(cos 24cos sin 2lt t B B N lt B simB N l B B N X x ''+-''+''++-''+''⋅''+=ρηηρρ5222425532233)5814185(cos 120)1(cos 6cos l t t t B N lt B N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ (8-42) 8.3.2高斯投影坐标反算公式x,y⇒B,l投影方程：),(),(21y x l y x B ϕϕ== (8-43)满足以下三个条件：①x 坐标轴投影后为中央子午线是投影的对称轴；② x 坐标轴投影后长度不变；③投影具有正形性质，即正形投影条件。

高斯投影正反算原理高斯投影是一种常用于地图制图的投影方式，也被广泛应用于其他领域的空间数据处理。

高斯投影正反算是对于已知的地球坐标系上的位置（经纬度），通过计算得到该点的平面坐标（东、北坐标），或者对于已知的平面坐标（东、北坐标），通过计算得到该点的地球坐标系上的位置（经纬度）的过程。

本文将详细介绍高斯投影正反算的原理。

一、高斯投影简介高斯投影是一种圆锥投影，其投影面在地球表面的某个经线上，也就是说，投影面是以该经线为轴的圆锥面。

经过对圆锥体的调整后，使其切于地球椭球面，在该经线上进行投影，同时保持沿该经线方向的比例尺一致，从而达到地图上各点在包括该经线的垂直面上映射的目的。

这种投影方式在某一特定区域内得到高精度的结果，因此广泛应用于地图制图。

二、高斯投影数学模型对于高斯投影正反算，需要先建立高斯投影坐标系与地球坐标系的转换模型。

1.高斯投影坐标系的建立高斯投影坐标系的建立需要确定圆锥面的基本参数，首先需要确定其所处的中央子午线，再确定该子午线上的经度为零点，并利用该经线上某一点的经度和该点的高度来确定该点所在的圆锥体。

圆锥体的底面包括所有与地球椭球面相切的圆面，通过对这些圆面进行调整，使得圆锥体转动后能够在中央子午线上进行投影。

在此基础上，可建立高斯投影坐标系，其中投影面为圆锥面，且中央子午线与投影面的交点称为该投影坐标系的中心，投影面的上端点和下端点分别对应正北方向和正南方向。

2.地球坐标系的建立地球坐标系是以地球椭球体为基础建立的，其坐标系原点确定为地球椭球体上的一个特定点。

在已知该点经纬度和高度的前提下，可确定以该点为中心的地球椭球体，并可根据它与地球坐标系之间的转换关系得到平面坐标系。

3.高斯投影坐标系与地球坐标系之间的转换关系由于高斯投影坐标系与地球坐标系存在不同的坐标体系和基准面，因此需要通过数学关系式来建立它们之间的转换关系。

（1）高斯投影坐标系转地球坐标系：已知高斯投影坐标系中任意一点的东北坐标（N，E），以及所属的中央子午线经度λ0、椭球参数a和e，则可通过以下公式求出该点的地球坐标系经纬度（φ，λ）和高度H：A0为以地球椭球体中心为原点，高斯投影坐标系中心投影坐标为（0，0）的点到椭球面的距离。

高斯投影正反算LT2340123423401234...........q m m y m y m y m y l n n y n y n y n y '''''=+++++'''''=+++++3.引入高斯投影条件之一：正形条件4.由于可得到00n '=，带入上式可得到 2402435135...........q m m y m y l n y n y n y '''=+++'''=+++5.引入高斯投影条件之三：中央子午线投影后长度不变001122223322442255sec sec 122sec (12)6sec (56)24sec (5286)120f f f f f f f f f f f f f f f f f f f m q B dm n dx N t B dn m dx N B n t N t B m t N B n t N ηηη'=⎧⎪'⎪'==⎪⎪⎪''=-=-⎪⎪⎪⎨'=-++⎪⎪⎪'=++⎪⎪⎪⎪'=++⎪⎩四、高斯投影的特点1.当l 等于常数时，随着B 的增加x 的值增大，y 的值减小，无论B 值为正或为负，y 值不变。

这就是说，椭球面上除中央子午线外，其它子午线投影后，均向中央子午线弯曲，并向两极收敛，同时还对称于中央子午线和赤道。

2.当B 等于常数时，随着l 的增加，x 值和y 值都增大。

所以在椭球面上对称于赤道的纬圈，投影后仍成为对称的曲线，同时与子午线的投影曲线互相垂直凹向两极。

3.据中央子午线越远的子午线，投影后弯曲越厉害，长度变形越大。

五、MATLAB编程实现坐标正反算1.编写main函数function maindisp('欢迎使用高斯投影正反算及相邻带的坐标换算程序');disp('1：高斯正算 2：高斯反算 3：换带计算');K=0;while (K<1||K>3)K=input('请根据上列选择计算类型 K=');switch Kcase 1GSZS;case 2GSFS;case 3HDJS;otherwisedisp('K 值无效（1-3）');enddisp('程序作者：亚里士多墩');disp('指导老师：亚里士多德');end2.编写高斯正算GSZS函数function GSZS%GSZS 是将大地坐标换算为高斯坐标的子函数%此函数要调用 DHH 和 HHD 两个子函数%此函数包含子午线收敛角的计算disp('你选择的是高斯正算');B=input('输入大地坐标 B=');L=input('输入大地坐标 L=');L0=input('输入所用中央子午线 L0=');B=DHH(B);L=DHH(L);L0=DHH(L0);disp('1:克拉索夫斯基椭球 2:1975 年国际椭球 3:WGS-84 椭球');T=0;while (T<1||T>3)T=input('请根据上列选择椭球模型 T=');switch Tcase 1a=6378245.0000000000;b=6356863.0187730473;X=111134.861*(B*180/pi)-16036.480*sin(2*B)+16.828*sin(4*B)-0.022*sin (6*B);case 2a=6378140.0000000000;b=6356755.2881575287;X=111133.005*(B*180/pi)-16038.528*sin(2*B)+16.833*sin(4*B)-0.022*sin (6*B);case 3a=6378137.0000000000;f=1/298.257223563;b=a*(1-f);X=6367449.1458*B-32009.8185*cos(B)*sin(B)-133.9975*cos(B)*(sin(B))^3-0.6975*cos(B)*(sin(B))^5;otherwisedisp('T 值无效（1-3）');endende=(sqrt(a^2-b^2))/a;e1=(sqrt(a^2-b^2))/b;V=sqrt(1+(e1^2)*(cos(B))^2);c=(a^2)/b;M=c/(V^3);N=c/V;t=tan(B);n=sqrt((e1^2)*(cos(B))^2);l=L-L0;xp1=X;xp2=(N*sin(B)*cos(B)*l^2)/2;xp3=(N*sin(B)*((cos(B))^3)*(5-t^2+9*n^2+4*n^4)*l^4)/24; xp4=(N*sin(B)*((cos(B))^5)*(61-58*t^2+t^4)*l^6)/720;x=xp1+xp2+xp3+xp4;yp1=N*cos(B)*l;yp2=N*(cos(B))^3*(1-t^2+n^2)*l^3/6;yp3=N*(cos(B))^5*(5-18*t^2+t^4)*l^5/120;y=yp1+yp2+yp3;r1=l*sin(B);r2=(1/3)*sin(B)*(cos(B))^2*(l^3)*(1+3*n^2+2*n^4);r3=(1/15)*sin(B)*(cos(B))^2*(l^5)*(2-t^2);r=r1+r2+r3;format long gxyR=HHD(r)End3.编写高斯反算公式GSFS函数function GSZS%GSZS 是将大地坐标换算为高斯坐标的子函数%此函数要调用 DHH 和 HHD 两个子函数%此函数包含子午线收敛角的计算disp('你选择的是高斯正算');B=input('输入大地坐标 B=');L=input('输入大地坐标 L=');L0=input('输入所用中央子午线 L0=');B=DHH(B);L=DHH(L);L0=DHH(L0);disp('1:克拉索夫斯基椭球 2:1975 年国际椭球 3:WGS-84 椭球');T=0;while (T<1||T>3)T=input('请根据上列选择椭球模型 T=');switch Tcase 1a=6378245.0000000000;b=6356863.0187730473;X=111134.861*(B*180/pi)-16036.480*sin(2*B)+16.828*sin(4*B)-0.022*sin (6*B);case 2a=6378140.0000000000;b=6356755.2881575287;X=111133.005*(B*180/pi)-16038.528*sin(2*B)+16.833*sin(4*B)-0.022*sin (6*B);case 3a=6378137.0000000000;f=1/298.257223563;b=a*(1-f);X=6367449.1458*B-32009.8185*cos(B)*sin(B)-133.9975*cos(B)*(sin(B))^3 -0.6975*cos(B)*(sin(B))^5;otherwisedisp('T 值无效（1-3）');endende=(sqrt(a^2-b^2))/a;e1=(sqrt(a^2-b^2))/b;V=sqrt(1+(e1^2)*(cos(B))^2);c=(a^2)/b;M=c/(V^3);N=c/V;t=tan(B);n=sqrt((e1^2)*(cos(B))^2);l=L-L0;xp1=X;xp2=(N*sin(B)*cos(B)*l^2)/2;xp3=(N*sin(B)*((cos(B))^3)*(5-t^2+9*n^2+4*n^4)*l^4)/24;xp4=(N*sin(B)*((cos(B))^5)*(61-58*t^2+t^4)*l^6)/720;x=xp1+xp2+xp3+xp4;yp1=N*cos(B)*l;yp2=N*(cos(B))^3*(1-t^2+n^2)*l^3/6;yp3=N*(cos(B))^5*(5-18*t^2+t^4)*l^5/120;y=yp1+yp2+yp3;r1=l*sin(B);r2=(1/3)*sin(B)*(cos(B))^2*(l^3)*(1+3*n^2+2*n^4); r3=(1/15)*sin(B)*(cos(B))^2*(l^5)*(2-t^2);r=r1+r2+r3;format long gxyR=HHD(r)End六、代码测试1.L=111°47'24〞.8974，B=31°04'41〞.6832，L0=111°2.L=111.47532575，B=31.23484275，L0=1113.L=114.20，B=30.30，L0=1114.L=118.54152206，B=32.24576522，L0=117。

高斯投影坐标反算c语言代码#include<stdio.h>#include<math.h>#include<conio.h>main(){printf("#####################################################\n");printf("# 角度输入说明：如26°12′45.2″输入为26,12,45.2 #\n");printf("#####################################################\n");double x,y;int j,L0;printf("请输入高斯投影坐标（自然坐标），中间用逗号隔开：\n");scanf("%lf,%lf",&x,&y);//自然坐标输入printf("请输入中央子午线L0:\n");scanf("%d,%d,%lf",&L0); //中央子午线输入printf("请选择参考椭球：1.北京1954参考椭球。

\n 2.西安1980参考椭球。

\n");printf("选择的参考椭球为：");scanf("%d",&j);//选择椭球参数if(j==1){long double Bf0=0.157046064172*pow(10,-6)*x;long double Bf=Bf0+cos(Bf0)*(0.005051773759*sin(Bf0)-0.000029837302*pow(sin(Bf0),3)+0.00000023818 9*pow(sin(Bf0),5));long double t=tan(Bf);long double m=0.00673852541468*pow(cos(Bf),2);long double V=1+m;long double N=6378245.000/sqrt(1-0.00669342162297*pow(sin(Bf),2));long double B1=Bf-1.0/2*V*t*pow(y/N,2)+1.0/24*(5+3*pow(t,2)+m-9*m*pow(t,2))*V*t*pow(y/N,4)-1.0/72 0*(61+90*pow(t,2)+45*pow(t,4))*V*t*pow(y/N,6);long double l1=(1/cos(Bf))*(y/N)-1.0/6*(1+2*pow(t,2)+m)*(1/cos(Bf))*pow(y/N,3)+1.0/120*(5+28*pow(t,2)+24*pow(t,4)+6*m+8*m*pow(t,2))*(1/cos(Bf))*pow(y/N,5);long double B=B1*57.29577951;long double l=l1*57.29577951;long double L=L0+l;int d2=int(B);int f2=int((B-d2)*60);long double m2=((B-d2)*60-f2)*60;printf("B=%d %d %.12lf\n",d2,f2,m2);int d3=int(L);int f3=int((L-d3)*60);long double m3=((L-d3)*60-f3)*60;printf("L=%d %d %.12lf\n",d3,f3,m3); //北京1954参考椭球}if(j==2){long double Bf0=0.157048687473*pow(10,-6)*x;long double Bf=Bf0+cos(Bf0)*(0.005052505593*sin(Bf0)-0.000029847335*pow(sin(Bf0),3)+0.0000002416* pow(sin(Bf0),5)+0.0000000022*pow(sin(Bf0),7));long double t=tan(Bf);long double m=0.006739501819473*pow(cos(Bf),2);long double V=1+m;long double N=6378140.000/sqrt(1-0.006694384999588*pow(sin(Bf),2));long double B1=Bf-1.0/2*V*t*pow(y/N,2)+1.0/24*(5+3*pow(t,2)+m-9*m*pow(t,2))*V*t*pow(y/N,4)-1.0/72 0*(61+90*pow(t,2)+45*pow(t,4))*V*t*pow(y/N,6);long double l1=(1/cos(Bf))*(y/N)-1.0/6*(1+2*pow(t,2)+m)*(1/cos(Bf))*pow(y/N,3)+1.0/120*(5+28*pow(t,2)+24*pow(t,4)+6*m+8*m*pow(t,2))*(1/cos(Bf))*pow(y/N,5);long double B=B1*57.29577951;long double l=l1*57.29577951;long double L=L0+l;int d2=int(B);int f2=int((B-d2)*60);long double m2=((B-d2)*60-f2)*60;printf("B=%d %d %.12lf\n",d2,f2,m2);int d3=int(L);int f3=int((L-d2)*60);long double m3=((L-d2)*60-f2)*60;printf("L=%d %d %.12lf\n",d3,f3,m3); //西安1980参考椭球}getch();}。

《测绘程序设计（）》上机实验报告（Visual C++.Net）班级：学号：姓名：序号：二零一一年五月实验7 常用测量程序设计1.实验目的：1.1 巩固类的创建与使用;1.2 掌握数组参数的传递；1.3 掌握常用测绘程序设计的技巧。

1.42.实验内容：编写高斯投影正反算程序。

3.设计思路：这次的实验目的是实现高斯正反算。

需要考虑投影方式即分带的方式，又要考虑椭球参数的类型，所以我添加了两个函数来完成此功能。

分别是int SetProjectType(int m)和void SetParameter(int m,double &a,double &b)。

4.界面设计：界面设计很简单，具体见运行结果。

5.主要代码：文件名：GaussProjectDlg.cpp代码：const double PI=4*atan(1.0);//获得分带方式返回中央子午线经度int CGaussProjectDlg::SetProjectType(int m){UpdateData(TRUE);int n; //记录分带带号double L; //经度L=iDegreeL+iMinL/60+dSecondL/3600;if(m==1) //6度带{n=int(L/6)+1;L0=6*n-3;}else if(m==2) //3度带{n=int((L+1.5)/3);L0=3*n;}else if(m==3) //自主分带L0=L0;return L0;}//获取椭球参数void CGaussProjectDlg::SetParameter(int m,double &a,double &b) {if(m==1) //克拉索夫斯基椭球{a=6378245.0;b=6356863.0187730473;//e=sqrt(0.006693421622966);}else if(m==2) //1975国际协议椭球{a=6378140.0;b=6356755.2881575287;//e=sqrt(0.006694384999588);}else if(m==3) //WGS-84椭球{a=6378137.0;b=6356752.3142;//e=sqrt(0.0066943799013);}}void CGaussProjectDlg::OnBnClickedButtonpositivecal(){// TODO: 在此添加控件通知处理程序代码UpdateData(TRUE);double N;double t;double Eta;double X;double A0,A2,A4,A6,A8;double RadB;double Rou;Rou=180*3600/PI;double a,b,e1,e2; //椭球参数SetParameter(iParameterType,a,b);e1=sqrt(a*a-b*b)/a;e2=sqrt(a*a-b*b)/b;double l;L0=SetProjectType(iProjectType);double L;L=iDegreeL+double(iMinL)/60+dSecondL/3600;l=(L-L0)*3600;RadB=(iDegreeB+double(iMinB)/60+dSecondB/3600)*PI/180;N=a/sqrt(1-e1*e1*sin(RadB)*sin(RadB));t=tan(RadB);Eta=e2*cos(RadB);A0=1+3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8);A2=-1.0/2*(3.0/4*e1*e1+60.0/64*pow(e1,4)+525.0/512*pow(e1,6)+17640.0/16384*pow(e1,8));A4=1.0/4*(15.0/64*pow(e1,4)+210.0/512*pow(e1,6)+8820.0/16384*pow(e1,8));A6=-1.0/6*(35.0/512*pow(e1,6)+2520.0/16384*pow(e1,8));A8=1.0/8*(315.0/16384*pow(e1,8));X=a*(1-e1*e1)*(A0*RadB+A2*sin(2*RadB)+A4*sin(4*RadB)+A6*sin(6*RadB)+A8*sin(8*RadB));x=X+N/(2*Rou*Rou)*sin(RadB)*cos(RadB)*l*l+N/(24*pow(Rou,4))*sin(RadB)*pow(cos(RadB),3)*(5-t*t+9*Eta*Eta+4*pow(Eta,4))*pow(l,4)+ N/(720*pow(Rou,6))*sin(RadB)*pow(cos(RadB),5)*(61-58*t*t+pow(t,4))*pow(l,6);y=N/Rou*cos(RadB)*l+N/(6*pow(Rou,3))*pow(cos(RadB),3)*(1-t*t+Eta*Eta)*pow(l,3)+N/(120*pow(Rou,5))*pow(cos(RadB),5)*(5-18*t*t+pow(t,4)+14*Eta*Eta-58*Eta*Eta*t*t)*pow( l,5);UpdateData(FALSE);}void CGaussProjectDlg::OnBnClickedButtonantical(){// TODO: 在此添加控件通知处理程序代码UpdateData(TRUE);double t_f;double Eta_f;double B_f;double N_f;double M_f;double X=x;double B0;double K0,K2,K4,K6;double a,b,e1,e2; //椭球参数SetParameter(iParameterType,a,b);e1=sqrt(a*a-b*b)/a;e2=sqrt(a*a-b*b)/b;double A0;A0=1+3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8);B0=X/(a*(1-e1*e1)*A0);K0=1.0/2*(3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8));K2=-1.0/3*(63.0/64*pow(e1,4)+1108.0/512*pow(e1,6)+58239.0/16384*pow(e1,8));K4=1.0/3*(604.0/512*pow(e1,6)+68484.0/16384*pow(e1,8));K6=-1.0/3*(26328.0/16384*pow(e1,8));B_f=B0+sin(2*B0)*(K0+sin(B0)*sin(B0)*(K4+K6*sin(B0)*sin(B0)));t_f=tan(B_f);Eta_f=e2*cos(B_f);N_f=a/sqrt(1-e1*e1*sin(B_f)*sin(B_f));M_f=N_f/(1+e2*e2*cos(B_f)*cos(B_f));double B;B=B_f-t_f/(2*M_f*N_f)*y*y+t_f/(24*M_f*pow(N_f,3))*(5+3*t_f*t_f+Eta_f*Eta_f-9*Eta_f*Eta_f*t_f*t_f)*pow(y,4)- t_f/(720*M_f*pow(N_f,5))*(61+90*t_f*t_f+45*pow(t_f,4))*pow(y,6);double l;l=1.0/(N_f*cos(B_f))*y-1.0/(6*pow(N_f,3)*cos(B_f))*(1+2*t_f*t_f+Eta_f*Eta_f)*pow(y,3)+1.0/(120*pow(N_f,5)*cos(B_f))*(5+28*t_f*t_f+24*pow(t_f,4)+6*Eta_f*Eta_f+8*Eta_f*Eta_f* t_f*t_f)*pow(y,5);//将B转化为度分秒的形式double dDegB;dDegB=B*180/PI;iDegreeB=int(dDegB);iMinB=int((dDegB-iDegreeB)*60);dSecondB=((dDegB-iDegreeB)*60-iMinB)*60;double dDegL;dDegL=l*180/PI+L0;iDegreeL=int(dDegL);iMinL=int((dDegL-iDegreeL)*60);dSecondL=((dDegL-iDegreeL)*60-iMinL)*60;UpdateData(FALSE);}6.运行结果：实验的运行结果如下图所示：正算：反算：7.实验总结这次实验是实现高斯投影的正反算，方法很多，实现并不复杂，但是计算公式复杂，变量繁多，稍有不慎，就会造成计算错误。

高斯投影正反算程序设计一．程序设计流程本程序的设计思路如下：（1），程序采用VS08版本作为开发平台，并采用C#语言作为开发语言，设计为WindowsForm窗体程序形式。

（2），程序主要的算法来自于教材。

但是本程序为了更加实用，添加了更多的解算基准，包括：WGS-84，国际椭球（3{{//{p ublicdoubledLongitude;publicdoubledLatitude;publicdoubledHeight;}//笛卡尔坐标//CartesianCoordinatepublicstructCRDCARTESIAN{publicdoublex;publicdoubley;publicdoublez;}publicForm1(){InitializeComponent();}privatevoidbutton1_Click(objectsender,EventArgse) {try{tt=}{}{ee=}{ee=}if(pareTo("1975国际椭球")==0){a=6378140.00;ee=}if(pareTo("2000国家大地坐标系")==0){a=6378137.0;}constdoublepai=3.1415926;doubleb=Math.Sqrt(a*a*(1-ee*ee));doublec=a*a/b;doubleepp=Math.Sqrt((a*a-b*b)/b/b);CRDGEODETICpcrdGeo;CRDCARTESIANpcrdCar;doublemidlong;//{}ai;//{}pai;//intnum;//带号midlong=0;//默认值，需要制定分带try{tt=}catch{MessageBox.Show("Choose3/6error!");return;if(pareTo("3度带")==0){num=Convert.ToInt32(deglon/6+1);midlong=(6*num-3)/180.0*pai;}if(pareTo("6度带")==0){num=Convert.ToInt32((deglon+1.5)/3);}a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8;a6=m6/32.0+m8/16.0;a8=m8/128.0;doubleB=pcrdGeo.dLatitude;doublesb=Math.Sin(B);doublecb=Math.Cos(B);doubles2b=sb*cb*2;doubles4b=s2b*(1-2*sb*sb)*2;doubles6b=s2b*Math.Sqrt(1-s4b*s4b)+s4b*Math.Sqrt(1-s2b*s2b); X=a0*B-a2/2.0*s2b+a4*s4b/4.0-a6/6.0*s6b;pcrdCar.x=Nscnb*lp*lp/2.0+Nscnb*cosb*cosb*Math.Pow(lp,4)*(5-t*t+9*ita*ita+4*Math.Pow(ita,4))/24.0 +Nscnb*Math.Pow(cosb,4)*Math.Pow(lp,6)*(61-58*t*t+Math.Pow(t,4))/720.0+X;pcrdCar.y=Ncosb*Math.Pow(lp,1)+Ncosb*cosb*cosb*(1-t*t+ita*ita)/6.0*Math.Pow(lp,3)+Ncosb*Math.Pow(l p,5)*Math.Pow(cosb,4)*(5-18*t*t+Math.Pow(t,4)+14*ita*ita-58*ita*ita*t*t)/120.0;if(pcrdCar.y<0)pcrdCar.y+=500000;richTextBox1.Text="Results:\nX:\t"+Convert.ToString(pcrdCar.x)+"\nY:\t"+Convert.ToString(pcrdCar.y );}{try{tt=}{}{ee=}if(pareTo("WGS-84")==0){a=6378137.00;ee=}if(pareTo("1975国际椭球")==0){a=6378140.00;ee=}if(pareTo("2000国家大地坐标系")==0) {a=6378137.0;ee}constdoublepai=doubleb=Math.Sqrt(a*a*(1-ee*ee));//求Xtry{tt=}{}if(pareTo("3度带")==0){midlong=num*3*pai/180;}if(pareTo("6度带")==0){midlong=(6*num-3)*pai/180;}b=Math.Sqrt(a*a*(1-ee*ee));c=a*a/b;epp=Math.Sqrt(a*a-b*b)/b;doublem0,m2,m4,m6,m8;doublea0,a2,a4,a6,a8;m0=a*(1-ee*ee);m2=3.0/2.0*m0*ee*ee;m4=5.0/4.0*ee*ee*m2;m6=7.0/6.0*ee*ee*m4;m8=9.0/8.0*ee*ee*m6;{}tf=Math.Tan(Bf);Vf=Math.Sqrt(1+epp*epp*Math.Cos(Bf)*Math.Cos(Bf));Nf=c/Vf;doubleynf=pcrdCar.y/Nf;pcrdGeo.dLatitude=Bf-1.0/2.0*Vf*Vf*tf*(ynf*ynf-1.0/12.0*Math.Pow(ynf,4)*(5+3*tf*tf+itaf*itaf-9*Math.Pow(itaf*tf,2))+1.0/360.0*(61+90*tf*tf+45*Math.Pow(tf,4))*Math.Pow(ynf,6));pcrdGeo.dLongitude=(ynf/Math.Cos(Bf)-(1+2*tf*tf+itaf*itaf)*Math.Pow(ynf,3)/6.0/Math.Cos(Bf)+ (5+28*tf*tf+24*Math.Pow(tf,4)+6*itaf*itaf+8*Math.Pow(itaf*tf,2))*Math.Pow(ynf,5)/120.0/Math.Cos(Bf ));pcrdGeo.dLongitude=pcrdGeo.dLongitude+midlong;//pcrdGeo.dLatitude=pcrdGeo.dLatitude;richTextBox2.Text="Results:\nLatitude:"+Convert.ToString(pcrdGeo.dLatitude)+"\nLongtitude:"+Conver t.ToString(pcrdGeo.dLongitude);}privatevoidlabel13_Click(objectsender,EventArgse){}}}米之间，。

高斯投影正反算编程一．高斯投影正反算基本公式（1）高斯正算基本公式（2）高斯反算基本公式以上主要通过大地测量学基础课程得到.这不进行详细的推导.只是列出基本公式指导编程的进行。

二．编程的基本方法和流程图（1）编程的基本方法高斯投影正反算基本上运用了所有的编程基本语句.本文中是利用C++语言进行基本的设计。

高斯正算中对椭球参数和带宽的选择主要运用了选择语句。

而高斯反算中除了选择语句的应用.在利用迭代算法求底点纬度还应用了循环语句。

编程中还应特别注意相关的度分秒和弧度之间的相互转换.这是极其重要的。

（2）相关流程图1）正算2）反算三．编程的相关代码（1）正算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))int i;struct jin{double B;double L;double L0;};struct jin g[100];main(int argc, double *argv[]) {FILE *r=fopen("a.txt","r"); assert(r!=NULL);FILE *w=fopen("b.txt","w"); assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].B,&g[i].L,&g[i].L 0)!=EOF){double a,b;int zuobiao;printf("\n请输入坐标系：北京54=1.西安80=2.WGS84=3：");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=6378245;b=6356863.0187730473;}if(zuobiao==2){a=6378140;b=6356755.2881575287;}if(zuobiao==3){a=6378137;b=6356752.3142;} //选择坐标系//double f=(a-b)/a;double e,e2;e=sqrt(2*f-f*f);e2=sqrt((a/b)*(a/b)-1);//求椭球的第一.第二曲率//double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;m0=a*(1-e*e);m2=3*e*e*m0/2;m4=5*e*e*m2/4;m6=7*e*e*m4/6;m8=9*e*e*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;double Bmiao,Lmiao, L0miao;Bmiao=(int)(g[i].B)*3600.0+(int)((g[i].B-(int)(g[i].B) )*100.0)*60.0+(g[i].B*100-(int)(g[i].B*100))*100.0;Lmiao=(int)(g[i].L)*3600.0+(int)((g[i].L-(int)(g[i].L) )*100.0)*60.0+(g[i].L*100-(int)(g[i].L*100))*100.0;L0miao=(int)(g[i].L0)*3600.0+(int)((g[i].L0-(int)(g[i] .L0))*100.0)*60.0+(g[i].L0*100-(int)(g[i].L0*100))*100 .0;double db;db=pi/180.0/3600.0;double B1,L1,l;B1=Bmiao*db;L1= Lmiao*db;l=L1-L0miao*db;//角度转化为弧度//double T=tan(B1)*tan(B1);double n=e2*e2*cos(B1)*cos(B1);double A=l*cos(B1);double X,x,y;X=a0*(B1)-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6 +a8*sin(8*B1)/8;//求弧长//double N=a/sqrt(1-e*e*sin(B1)*sin(B1));int Zonewide;int Zonenumber;printf("\n请输入带宽：3度带或6度带Zonewide=");scanf("%d",&Zonewide);getchar();if(Zonewide==3){Zonenumber=(int)((g[i].L-Zonewide/2)/Zonewide+1);}else if(Zonewide==6){Zonenumber=(int)g[i].L/Zonewide+1;}else{printf("错误");exit(0);}//选择带宽//doubleFE=Zonenumber*1000000+500000;//改写为国家通用坐标//y=FE+N*A+A*A*A*N*(1-T*T+n*n)/6+A*A*A*A*A*N*(5-18*T*T+T *T*T*T+14*n*n-58*n*n*T*T)/120;x=X+tan(B1)*N*A*A/2+tan(B1)*N*A*A*A*A*(5-T*T+9*n*n+4*n *n*n*n)/24+tan(B1)*N*A*A*A*A*A*A*(61-58*T*T+T*T*T*T)/7 20;printf("\n所选坐标系的转换结果：x=%lf y=%lf\n",x,y);fprintf(w,"%lf %lf\n",x,y);//输出结果到文本文件//}fclose(r);fclose(w);system("pause");return 0;}（2）反算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))double X,Y,B1,B2,B3,F,t;double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8,a1,b1;double BB,LL,Bf;double e,e1;int d,m,s,i,zuobiao;double sort(double,double);struct jin{double x;double y;double L0;};struct jin g[100];//x,y,L0为输入量：x,y坐标和中央子午线经度//main(int argc, double *argv[]){FILE *r=fopen("c.txt","r");assert(r!=NULL);FILE *w=fopen("d.txt","w");assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].x,&g[i].y,&g[i].L 0)!=EOF)//文件为空.无法打开//{double a1=6378245.0000000000;//克拉索夫斯基椭球参数//double b1=6356863.0187730473;double a75=6378140.0000000000;//1975国际椭球参数//double b75=6356755.2881575287;double a84=6378137.0000000000;//WGS-84系椭球参数//double b84=6356752.3142000000;double M,N;//mouyou圈曲率半径.子午圈曲率半径// double t,n;double A,B,C;double BB,LL,Bf,LL0,BB0;double a,b;printf("\n选择参考椭球：1=克拉索夫斯基椭球.2=1975国际椭球.3=WGS-84系椭球:");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=a1;b=b1;}if(zuobiao==2){a=a75;b=b75;}if(zuobiao==3){a=a84;b=b84;}//选择参考椭球.求解第一偏心率e,第二偏心率e1// Bf=sort(a,b);//调用求解底点纬度的函数//double q=sqrt(1-e*e*sin(Bf)*sin(Bf));double G=cos(Bf);M=a*(1-e*e)/(q*q*q);N=a/q;double H,I;A=g[i].y/N;H=A*A*A;I=A*A*A*A*A;t=tan(Bf);n=e1*cos(Bf);B=t*t;C=n*n;BB0=Bf-g[i].y*t*A/(2*M)+g[i].y*t*H/(24*M)*(5+3*B+C-9*B *C)-g[i].y*t*I/(720*M)*(61+90*B+45*B*B);LL0=g[i].L0*pi/180.0+A/G-H/(6*G)*(1.0+2*B+C)+I/(120*G)*(5.0+28*B+24*B*B+6*C+8*B*C);//利用公式求解经纬度// int Bdu,Bfen,Ldu,Lfen;double Bmiao,Lmiao;Ldu=int(LL0/pi*180);Lfen=int((LL0/pi*180)*60-Ldu*60);Lmiao=LL0/pi*180*3600-Ldu*3600-Lfen*60;Bdu=int(BB0/pi*180);Bfen=int((BB0/pi*180)*60-Bdu*60);Bmiao=BB0/pi*180*3600-Bdu*3600-Bfen*60;//将弧度转化为角度//printf("\n所选坐标系的转换结果：%d度%d分%lf秒 %d 度%d分%lf秒 \n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);fprintf(w,"%d°%d’%lf”%d°%d’%lf”\n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);//将结果输出到文本文件//}fclose(r);fclose(w);system("pause");return 0;}double sort(double a,double b){double e,e1;e=sqrt(1-(b/a)*(b/a));e1=sqrt((a/b)*(a/b)-1);double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;m0=a*(1-e*e);m2=3*e*e*m0/2;m4=5*e*e*m2/4;m6=7*e*e*m4/6;m8=9*e*e*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;B1=g[i].x/a0;do{F=-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6+a8*s in(8*B1)/8;B2=(g[i].x-F)/a0;B3=B1;B1=B2;} while(fabs(B3-B2)>10e-10);//利用迭代算法求解底点纬度//return B2;}。