高斯投影坐标反算c语言代码

合集下载

高斯投影正反算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)截图.。

[转载]高斯正反算

[转载]⾼斯正反算⼤地坐标向笛卡尔坐标转换⾼斯正反算采⽤不同椭球实现⾼斯克⾥格投影，将经纬度坐标转换为⾼斯平⾯坐标：正算⾼斯平⾯坐标转换为不同椭球下的经纬度坐标：反算1void GaussProjectDirect(double a,double efang,double B,double L,double L0,double& x,double &y,double& R)//⾼斯投影正算克⽒2 {34double b=aefangtob(a,efang);5double e2=seconde(a,b);6double W=sqrt(1-efang*sin(B)*sin(B));printf("W=%f",W);7double N=a/W;printf("N=%f",N);8double M=a*(1-efang)/pow(W,3);printf("M=%f",M);9double t=tee(B);10double eitef=eitefang(a,b,B);11double l=L-L0;12//主曲率半径计算13double m0,m2,m4,m6,m8,n0,n2,n4,n6,n8;14 m0=a*(1-efang); n0=a;15 m2=3.0/2.0*efang*m0; n2=1.0/2.0*efang*n0;16 m4=5.0/4.0*efang*m2; n4=3.0/4.0*efang*n2;17 m6=7.0/6.0*efang*m4; n6=5.0/6.0*efang*n4;18 m8=9.0/8.0*efang*m6; n8=7.0/8.0*efang*n6;19//⼦午线曲率半径20double a0,a2,a4,a6,a8;21 a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;22 a2=m2/2+m4/2+15*m6/32+7.0/16.0*m8;23 a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8;24 a6=m6/32+m8/16;25 a8=m8/128;2627double X=a0*B-a2/2*sin(2*B)+a4/4*sin(4*B)-a6/6*sin(6*B)+a8/8*sin(8*B);28 x=X+N/2*t*cos(B)*cos(B)*l*l+N/24*t*(5-t*t+9*eitef+4*pow(eitef,2))*pow(cos(B),4)*pow(l,4)+N/720*t*(61-58*t*t+pow(t,4))*pow(cos(B),6)*pow(l,6);29 y=N*cos(B)*l+N/6*(1-t*t+eitef)*pow(cos(B),3)*pow(l,3)+N/120*(5-18*t*t+pow(t,4)+14*eitef-58*eitef*t*t)*pow(cos(B),5)*pow(l,5);30 R=sqrt(M*N);31 }323334//⾼斯投影反算353637void GaussProjectInvert(double a,double efang,double x,double y,double L0,double &B,double& L,double& R)38 {39double b=aefangtob(a,efang);404142double m0,m2,m4,m6,m8,n0,n2,n4,n6,n8;43 m0=a*(1-efang); n0=a;44 m2=3.0/2.0*efang*m0; n2=1.0/2.0*efang*n0;45 m4=5.0/4.0*efang*m2; n4=3.0/4.0*efang*n2;46 m6=7.0/6.0*efang*m4; n6=5.0/6.0*efang*n4;47 m8=9.0/8.0*efang*m6; n8=7.0/8.0*efang*n6;484950//⼦午线曲率半径51double a0,a2,a4,a6,a8;52 a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;53 a2=m2/2+m4/2+15*m6/32+7.0/16.0*m8;54 a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8;55 a6=m6/32+m8/16;56 a8=m8/128;575859double X=x;60double FBf=0;61double Bf0=X/a0,Bf1=0;62while((Bf0-Bf1)>=0.0001)63 { Bf1=Bf0;64 FBf=a0*Bf0-a2/2*sin(2*Bf0)+a4/4*sin(4*Bf0)-a6/6*sin(6*Bf0)+a8/8*sin(8*Bf0);65 Bf0=(X-FBf)/a0;66 }67double Bf=Bf0;68double Vf=bigv(a,b,Bf);69double tf=tee(Bf);70double Nf=bign(a,b,Bf);71double eiteffang=eitefang(a,b,Bf);72double Bdu=rad_deg(Bf)-1/2.0*Vf*Vf*tf*(pow((y/Nf),2)-1.0/12*(5+3*tf*tf+eiteffang-9*eiteffang*tf*tf)*pow((y/Nf),4)+1.0/360.0*(61+90*tf*tf+45*tf*tf)*pow((y/Nf),6))*180/PI; 73double ldu=1.0/cos(Bf)*(y/Nf+1.0/6.0*(1+2*tf*tf+eiteffang)*pow((y/Nf),3)+1.0/120.0*(5+28*tf*tf+24*tf*tf+6*eiteffang+8*eiteffang*tf*tf)*pow((y/Nf),5))*180.0/PI;747576 B=deg_int(Bdu);77 L=L0+deg_int(ldu);78double W=sqrt(1-efang*sin(B)*sin(B));printf("W=%f\n",W); 79double N=a/W;printf("N=%f\n",N);80double M=a*(1-efang)/pow(W,3);printf("M=%f\n",M);81 R=sqrt(M*N);828384 }。

(完整word版)高斯投影正反算代码

l=(1-(b3-b5*Z*Z)*Z*Z)*Z;
L=l+L01; ///
反算就出
L B
l=L-L02;
B=Bf-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2;
t=tan(B);
n2=e2*cos(B)*cos(B);
V=sqrt(1+n2);
cB2=pow(cos(B),2);
N=6399698.902-(21562.267-(108.973-0.612*cB2)*cB2)*cB2; //N=c/V;
l=L-111*3600/P // l=((m%6)*3600+n*60+h)/P;
t=tan(B);
n2=e2*cos(B)*cos(B);
V=sqrt(1+n2);
cB2=pow(cos(B),2);
N=6399698.902-(21562.267-(108.973-0.612*cB2)*cB2)*cB2; // N=c/V;
a0=32140.404-(135.3302-(0.7092-0.004*cB2)*cB2)*cB2;
a4=(0.25+0.00252*cB2)*cB2-0.04166;
a6=(0.166*cB2-0.084)*cB2;
a3=(0.3333333+0.001123*cB2)*cB2-0.1666667;
#include "stdafx.h"
#include "iostream.h"
#include "math.h"
#include "stdio.h"
#define P 206264.806247096355

C#高斯正算高斯反算高斯换带等

C#⾼斯正算⾼斯反算⾼斯换带等⾸先，你要确定椭球参数：C#代码1. a = 6378140; //西安80椭球 IGA752. e2 = 0.006694384999588;3. m0 = a * (1 - e2);4. m2 = 3.0 / 2 * e2 * m0;5. m4 = 5.0 / 4 * e2 * m2;6. m6 =7.0 / 6 * e2 * m4;7. m8 = 9.0 / 8 * e2 * m6;8. a0 = m0 + m2 / 2 + (3.0 / 8.0) * m4 + (5.0 / 16.0) * m6 + (35.0 / 128.0) * m8;9. a2 = m2 / 2 + m4 / 2 + 15.0 / 32 * m6 + 7.0 / 16 * m8;10. a4 = m4 / 8 + 3.0 / 16 * m6 + 7.0 / 32 * m8;11. a6 = m6 / 32 + m8 / 16;12. a8 = m8 / 128;13. xx = 0;14. yy = 0;15. _x = 0;16. _y = 0;17. BB = 0;18. LL = 0;下⾯才开始正题：⾼斯正算：把经纬度坐标转换为平⾯坐标C#代码1. void GaussPositive(double B, double L, double L0)2. {3. double X, t, N, h2, l, m, Bmiao, Lmiao;4. int Bdu, Bfen, Ldu, Lfen;5. Bdu = (int)B;6. Bfen = (int)(B * 100) % 100;7. Bmiao = (B - Bdu - Bfen * 0.01) * 10000.0;8. B = Bdu * PI / 180.0 + (Bfen / 60.0) * PI / 180.0 + Bmiao / 3600.0 * PI / 180.0;9. Ldu = (int)L;10. Lfen = (int)(L * 100) % 100;11. Lmiao = (L - Ldu - Lfen * 0.01) * 10000.0;12. L = Ldu * PI / 180.0 + (Lfen / 60.0) * PI / 180 + Lmiao / 3600.0 * PI / 180.0;13. l = L - L0 * PI / 180;14. X = a0 * B - Math.Sin(B) * Math.Cos(B) * ((a2 - a4 + a6) + (2 * a4 - 16.0 / 3.0 * a6) * Math.Sin(B) * Math.Sin(B) + 16.0 / 3.0 * a6 * Math.Pow(Math.Sin(B), 4)) + a8 / 8.0 * Math.Sin(8 * B);15. t = Math.Tan(B);16. h2 = e2 / (1 - e2) * Math.Cos(B) * Math.Cos(B);17. N = a / Math.Sqrt(1 - e2 * Math.Sin(B) * Math.Sin(B));18. m = Math.Cos(B) * l;19. xx = X + N * t * ((0.5 + (1.0 / 24.0 * (5 - t * t + 9 * h2 + 4 * h2 * h2) + 1.0 / 720.0 * (61 - 58 * t * t + Math.Pow(t, 4)) * m * m) * m * m) * m * m);20. yy = N * ((1 + (1.0 / 6.0 * (1 - t * t + h2) + 1.0 / 120.0 * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * h2 - 58 * h2 * t * t) * m * m) * m * m) * m);21. yy = yy + 500000;22. }⾼斯反算：把平⾯坐标转换为经纬度坐标：C#代码1. void GaussNegative(double x, double y, double L0)2.3. double Bf, Vf, l, tf, hf2, Nf, Bmiao, Lmiao;4. int Bdu, Bfen, Ldu, Lfen;4. int Bdu, Bfen, Ldu, Lfen;5. y = y - 500000;6. Bf = hcfansuan(x);7. Vf = Math.Sqrt(1 + e2 / (1 - e2) * Math.Cos(Bf) * Math.Cos(Bf));8. tf = Math.Tan(Bf);9. hf2 = e2 / (1 - e2) * Math.Cos(Bf) * Math.Cos(Bf);10. Nf = a / Math.Sqrt(1 - e2 * Math.Sin(Bf) * Math.Sin(Bf));11. BB = (Bf - 0.5 * Vf * Vf * tf * (Math.Pow(y / Nf, 2) - 1.0 / 12 * (5 + 3 * tf * tf + hf2 - 9 * hf2 * tf * tf) * Math.Pow(y / Nf, 4) + 1.0 / 360 * (61 + 90 * tf * tf + 45 * tf * tf) * Math.Pow(y / Nf, 6))) * 180 / PI;12. Bdu = (int)BB;13. Bfen = (int)((BB - Bdu) * 60);14. Bmiao = ((BB - Bdu) * 60 - Bfen) * 60;15. BB = Bdu + 0.01 * Bfen + 0.0001 * Bmiao;16. l = 1.0 / Math.Cos(Bf) * (y / Nf - 1.0 / 6.0 * (1 + 2 * tf * tf + hf2) * Math.Pow(y / Nf, 3) + 1.0 / 120.0 * (5 + 28 * tf * tf + 24 * Math.Pow(tf, 4) + 6 * hf2 + 8 * hf2 * tf * tf) * Math.Pow(y / Nf, 5)) * 180.0 / PI;17. LL = L0 + l;18. Ldu = (int)LL;19. Lfen = (int)((LL - Ldu) * 60);20. Lmiao = ((LL - Ldu) * 60 - Lfen) * 60;21. LL = Ldu + 0.01 * Lfen + 0.0001 * Lmiao;⾥⾯涉及到的弧长反算：C#代码1. double hcfansuan(double pX)2. {3. double Bf0 = pX / a0;4. double Bf1, Bf2;5. Bf1 = Bf0;6. Bf2 = (pX - F(Bf1)) / a0;7. while ((Bf2 - Bf1) > 1.0E-11)8. {9. Bf1 = Bf2;10. Bf2 = (pX - F(Bf1)) / a0;11. }12. return Bf1;13. }⾼斯换带就⽐较简单了：C#代码1. void GaussZone(double x, double y, double L0, double L0new)2. {3. GaussNegative(x, y, L0);4. GaussPositive(BB, LL, L0new);5. }。

高斯转换大地坐标C代码

C代码：/* 参数说明所有投影转换公式都基于椭球体椭球体长半轴北京54 a = 6378245西安80 a = 6378140WGS 84 a = 6378137椭球体短半轴北京54 b = 6356863.0188西安80 b = 6356755.2882WGS 84 b = 6356752.3142扁率 f = (a - b) / a第一偏心率 e = sqrt(1 - (b / a) * (b / a))第二偏心率 e' = sqrt((a / b) * (a / b) - 1)卯酉圈曲率半径N = (a * a / b) / sqrt(1 + e' * e' * cosβ * cosβ子午圈曲率半径 R = a * (1 - e * e) / sqrt((1 - e * e * sinβ * sinβ) * (1 - e * e * sinβ * sinβ) * (1 - e * e * sinβ * sinβ))B -- 纬度，L -- 经度，单位弧度(RAD)N X -- 纵直角坐标, E Y -- 横直角坐标，单位米(M)//*/// 高斯投影转换void GaussCalc(double L, double B, double& X, double& Y){// 椭球体短半轴//double a = 6378245; // 北京54//double b = 6356863.0188; // 北京54double a = 6378140; // 西安80double b = 6356755.2882; // 西安80//double a = 6378137; // WGS 84//double b = 6356752.3142; // WGS 84double f = (a - b) / a; // 扁率double e = 1 - (b / a) * (b / a); // 第一偏心率平方double e2 = (a / b) * (a / b) - 1; // 第二偏心率平方double dblPI = 3.1415926535898 / 180.0; // 角度弧度转换参数// 求经度所在带号int ZoneWide = 3; // 带宽3度带或6度带int ZoneNumber; // 带号double L0; // 中央经度if (3 == ZoneWide){// 3度带ZoneNumber = int(L - ZoneWide / 2) / ZoneWide + 1;}else// if (6 == ZoneWide){// 6度带ZoneNumber = int(L) / ZoneWide + 1;}// 中央经度L0 = (ZoneNumber - 1) * ZoneWide + ZoneWide / 2;L0 = 120;// 角度转换成弧度double L1 = L * dblPI;double B1 = B * dblPI;L0 = L0 * dblPI;// 子午圈曲率半径double R = a * (1 - e) / sqrt((1 - e * sin(B1) * sin(B1)) * (1 - e * sin(B1) * sin(B1)) * (1 - e * sin(B1) * sin(B1)));double T = tan(B1) * tan(B1);double C = e2 * cos(B1) * cos(B1);double A = (L1 - L0) * cos(B1);// 子午线弧长double M = a * ((1 - e / 4 - 3 * e * e / 64 - 5 * e * e * e / 256) * B1- (3 * e / 8 + 3 * e * e / 32 + 45 * e * e * e / 1024) * sin(2 * B1)+ (15 * e * e / 256 + 45 * e * e * e / 1024) * sin(4 * B1)- (35 * e * e * e / 3072) * sin(6 * B1));// 卯酉圈曲率半径double N = a / sqrt(1 - e * sin(B1) * sin(B1));// 东纬度偏移double FE = ZoneNumber * 1000000L + 500000L;// 高斯- 克吕格投影比例因子double k0 = 1;// 根据投影公式计算点在投影坐标系中的位置X = FE + k0 * N * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 14 * C - 58 * T * C) * A * A * A * A * A / 120);Y = k0 * (M + N * tan(B1) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24)+ (61 - 58 * T + T * T + 270 * C - 330 * T * C) * A * A * A * A * A * A / 720);}#define RAD 57.295779 //rad = degree * 3.1415926 / 180 ;=>rad = degree / RADvoid Gauss2( double L , double B , double& x , double& y ){// 椭球体//double a = 6378245; // 北京54//double b = 6356863.0188; // 北京54double a = 6378140; // 西安80double b = 6356755.2882; // 西安80//double a = 6378137; // WGS 84//double b = 6356752.3142; // WGS 84double f = (a - b) / a; // 扁率double e = 1 - (b / a) * (b / a); // 第一偏心率平方double e2 = (a / b) * (a / b) - 1; // 第二偏心率平方// 求经度所在带号int ZoneWide = 3; // 带宽3度带或6度带int ZoneNumber; // 带号double L0; // 中央经度if (3 == ZoneWide){// 3度带ZoneNumber = int(L - ZoneWide / 2) / ZoneWide + 1;L0 = ZoneNumber * ZoneWide ;}else// if (6 == ZoneWide){// 6度带ZoneNumber = int(L) / ZoneWide + 1;L0 = ZoneNumber*ZoneWide - 3;}//L、B转化成度L = int( L ) + (int( L*100 ) - int( L )*100)/double(60) + (L*10000 - int(L*100)*100)/3600;B = int( B ) + (int( B*100 ) - int( B )*100)/double(60) + (B*10000 - int(B*100)*100)/3600;double A = (L - L0)/RAD;double TB = tan( B/RAD );double TB2 = TB*TB;double CB = cos( B/RAD );double C = e2*CB*CB;double C2 = 1 + C;double N = (a*a/b) / sqrt( C2 );//卯酉曲率半径double P = A*A*CB*CB;double P2 = sin( B/RAD );double Q = P2*P2;double R = 32005.78006+Q*(133.92133+Q*0.7031);x = 6367558.49686*B/RAD - P2*CB*R + ( ( ( (TB2-58)*TB2 + 61)*P/30+(4*N+5)*C2-TB2)*P/12+1)*N*CB*P/2;y = 500000+( ( ( (TB2-18)*TB2-(58*TB2-14)*N+5)*P/20+C2-TB2)*P/6+1)*N*(A*CB);//为保证为正数，向西便偏移500公里，即＋500,000；}。

高斯投影坐标正反算编程报告

高斯投影坐标正反算编程报告1.编程思想高斯投影坐标正反算的编程需要涉及大量公式。

为了使程序更清晰，各模块的数据重用性更强，本文采用结构化编程思想。

程序由四大块组成。

大地测量家庭作业。

Cpp文件用于存储main（）函数，是整个程序的入口。

尝试通过结构化编程简化main（）函数。

myfunction.h和myfunction.cpp用于存放计算过程中进行角度弧度换算时所要用到的一些自定的转换函数。

正蒜。

H和zhengsuan CPP用于存储zhengsuan类。

在正算类中，声明了用于高斯投影坐标正演计算的所有变量。

成员变量的初始化和正向计算在类的构造函数中执行。

通过get函数获得相应的阳性结果。

fansuan.h和fansuan.cpp用于存放fansuan类，类似于zhengsuan类，fansuan类中声明了高斯投影坐标反算所要用到的所有变量，在类的构造函数中进行成员变量的初始化及反算计算。

通过get函数获得相应的反算结果。

2.计算模型高斯投影正算公式十、十、nn232244？？？？sinbcosb？Lsimbcosb（5？t？9？4？）l2？？？224??? 4n5246？？sinbcosb（61？58t？t）l720？？？6y？？nn3223？？cosb？Lcosb（1？t？？）L6.3n5242225？？cosb（5？18吨？14？58吨）l120？？？五tf2mfnf5f高斯投影反算公式B男朋友？？tfy？2tf24mfn3f？5.3t2f224？？2f？9? ftfy？720mfn46y61？90t2？45泰夫？？yy32l？？1.2t2f？？f3nfcosbf6nfcosbf？？？y524222?5?28t?24t?6??8?fffftf5120nfcosbf?3.程序框图一开始输入b，l求定带号n，中央纬度l0，纬度差l正算按照实用公式计算x,y换算为国家统一坐标x，y输出x,y输入国家统一坐标x,y由y取定带号n，并换算出x，y求出中央经线l0反算按照实用公式计算b，ll=l0+l求出大地经度l输出b，l结束4.计算结果25.附录：程序代码/////主函数入口大地测量家庭作业。

高斯正反算

#include <stdio.h>#include <math.h>int main (void){ int o;printf ("\n============================欢迎使用高斯投影计算程序============================\n");printf(" 说明\n");printf("==此程序可实现克拉索夫斯基椭球和1975国际椭球上进行高斯三度带和六度带正算和反算==\n");printf ("\n");printf ("=========================请选择正反算。

1为正算，2为反算。

=======================\n");scanf ("%d",&o);if (o==1){int L0,i,e,a,b,c,f,g,L1,u,w;float d,h;double pi = 3.141592653;double P = 206264.806247;double x,y,l,N,B,a0,a4,a6,a3,a5,B2,C,L;printf ("=======================请输入分带方法。

1为3度带，2为6度带=======================\n");scanf ("%d",&u);printf ("=====================请输入经纬度。

例：1，00，00；31，00，00====================\n");scanf ("%d,%d,%d;%d,%d,%d",&a,&b,&c,&e,&f,&g);if (u==1){L1 = (e/3)*3;};if (u==2){L1 = (((e+f/60+g/3600)/6)+1)*6-3;};d = a+(b/60)+((c/60)/60);B = (d * pi)/180;C = cos(B)*cos(B);L = e*3600+f*60+g;B2 = a*3600+b*60+c;L0 = L1*3600;l = (L-L0)/P;printf ("================请选择椭球。

高斯投影正算

高斯投影正、反算代码//高斯投影正、反算//////6度带宽 54年北京坐标系//高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres)void GaussProjCal(double longitude, double latitude, double *X, double *Y){int ProjNo=0; int ZoneWide; ////带宽double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;double a,f, e2,ee, NN, T,C,A, M, iPI;iPI = 0.0174532925199433; ////3.1415926535898/180.0;ZoneWide = 6; ////6度带宽a=6378245.0; f=1.0/298.3; //54年北京坐标系参数////a=6378140.0; f=1/298.257; //80年西安坐标系参数ProjNo = (int)(longitude / ZoneWide) ;longitude0 = ProjNo * ZoneWide + ZoneWide / 2;longitude0 = longitude0 * iPI ;latitude0=0;longitude1 = longitude * iPI ; //经度转换为弧度latitude1 = latitude * iPI ; //纬度转换为弧度e2=2*f-f*f;ee=e2*(1.0-e2);NN=a/sqrt(1.0-e2*sin(latitude1)*sin(latitude1));T=tan(latitude1)*tan(latitude1);C=ee*cos(latitude1)*cos(latitude1);A=(longitude1-longitude0)*cos(latitude1);M=a*((1-e2/4-3*e2*e2/64-5*e2*e2*e2/256)*latitude1-(3*e2/8+3*e2*e2 /32+45*e2*e2*e2/1024)*sin(2*latitude1)+(15*e2*e2/256+45*e2*e2*e2/1024)*sin(4*latitude1)-(35*e2*e2*e2/30 72)*sin(6*latitude1));xval = NN*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/120);yval = M+NN*tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24+(61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720);X0 = 1000000L*(ProjNo+1)+500000L;Y0 = 0;xval = xval+X0; yval = yval+Y0;*X = xval;*Y = yval;}//高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)void GaussProjInvCal(double X, double Y, double *longitude, double *latitude){int ProjNo; int ZoneWide; ////带宽double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;iPI = 0.0174532925199433; ////3.1415926535898/180.0;a = 6378245.0; f = 1.0/298.3; //54年北京坐标系参数////a=6378140.0; f=1/298.257; //80年西安坐标系参数ZoneWide = 6; ////6度带宽ProjNo = (int)(X/1000000L) ; //查找带号longitude0 = (ProjNo-1) * ZoneWide + ZoneWide / 2;longitude0 = longitude0 * iPI ; //中央经线X0 = ProjNo*1000000L+500000L;Y0 = 0;xval = X-X0; yval = Y-Y0; //带内大地坐标e2 = 2*f-f*f;e1 = (1.0-sqrt(1-e2))/(1.0+sqrt(1-e2));ee = e2/(1-e2);M = yval;u = M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256));fai = u+(3*e1/2-27*e1*e1*e1/32)*sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*si n(4*u)+(151*e1*e1*e1/96)*sin(6*u)+(1097*e1*e1*e1*e1/512)*sin(8*u);C = ee*cos(fai)*cos(fai);T = tan(fai)*tan(fai);NN = a/sqrt(1.0-e2*sin(fai)*sin(fai));R = a*(1-e2)/sqrt((1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai))*(1-e2 *sin(fai)*sin(fai)));D = xval/NN;//计算经度(Longitude) 纬度(Latitude)longitude1 = longitude0+(D-(1+2*T+C)*D*D*D/6+(5-2*C+28*T-3*C*C+8*ee+24*T*T)*D*D*D*D*D/120)/cos(fai);latitude1 = fai -(NN*tan(fai)/R)*(D*D/2-(5+3*T+10*C-4*C*C-9*ee)*D*D*D*D/24+(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720);//转换为度 DD*longitude = longitude1 / iPI;*latitude = latitude1 / iPI;}。

大地高斯正反算C

大地高斯正反算-C#代码如下using System;using System.Collections.Generic;using SystemponentModel;using System.Data;using System.Drawing;using System.Linq;using System.Text;using System.Windows.Forms;namespace WindowsFormsApplication3public partial class Form1 : Formdouble a, E,E1,bbb;public Form1()InitializeComponent();private void menuStrip1_ItemClicked(object sender, ToolStripItemClickedEventArgs e)private void toolStripTextBox1_Click(object sender, EventArgs e)private void toolStripMenuItem1_Click(object sender, EventArgs e)private void button1_Click(object sender,EventArgs e)if (a != 6378245.00 && a != 6378137 && a != 6378137 && a!=6378140)textBox3.Text = "请?选?择?坐?标括?系μ";textBox4.Text = "请?选?择?坐?标括?系μ";elseint du, fen, miao;double B, L, X;double N, x, y, t, n, m, m0, m2, m4, m6, m8,l;B = Convert.ToDouble(textBox1.Text) * Math.PI / 180 + Convert.ToDouble(textBox2.Text)/60 * Math.PI / 180 + Convert.ToDouble(textBox9.Text)/3600 * Math.PI / 180;L = Convert.ToDouble(textBox10.Text) * Math.PI / 180 + Convert.ToDouble(textBox11.Text)/60 * Math.PI / 180 + Convert.ToDouble(textBox12.Text)/3600 * Math.PI / 180; ;l = L - Convert.ToInt16((L-3.0) / 6) * 6-3.0 ;l = l * Math.PI / 180;N = a / Math.Pow(1 - E * Math.Sin(B) * Math.Sin(B), 0.5);m0 = a * (1 - E);m2 = 1.5 * E * m0;m4 = 5 / 4 * E * m2;m6 = 7 / 6 * E * m4;m8 = 9 / 8 * E * m6;X = (m0 + 0.5 * m2 + 3 / 8 * m4 + 5 / 16 * m6 + 35 / 128 * m8) * B - (0.5 * m2 + 0.5 * m4 + 15 / 32 * m6 + 7 / 16 * m8) / 2 * Math.Sin(2 * B) + (m4 / 8 + 3 / 16 * m6 + 7 / 32 * m8) / 4 * Math.Sin(4 * B) - (m6 / 32 + m8 / 16) / 6 * Math.Sin(6 * B) + m8 / 128 / 8 * Math.Sin(8 * B);t = Math.Tan(B);n = Math.Cos(B) * Math.Cos(B) * E1;m = Math.Cos(B) * l;x = X + N / 2 * t * Math.Cos(B) * Math.Cos(B) * l * l + N / 24 * t * (5 - t * t + 9 * n + 4 * n * n) * Math.Pow(Math.Cos(B) * l, 4) + N / 720 * t * (61 - 58 * t * t + t * t * t * t) * Math.Pow(Math.Cos(B) * l, 6);y = N * ((1 + (1 / 6 * (1 - t * t + n) + 1 / 120 * (5 - 18 * t * t + t * t * t * t + 14 * n - 58 * n * t * t) * m * m)* m * m) * m);textBox3.Text = Convert.ToString(x);textBox4.Text = Convert.ToString(y);private void button2_Click(object sender, EventArgs e)int du, fen, miao;if (a != 6378245.00 && a != 6378137 && a != 6378137 && a != 6378140)textBox7.Text = "请?选?择?坐?标括?系μ";textBox8.Text = "请?选?择?坐?标括?系μ";elsedouble Bf, bb, Vf, tf, xf, yf, nf, Nf, B0, l0;xf = Convert.ToDouble(textBox5.Text);yf = Convert.ToDouble(textBox6.Text);bb = xf / bbb;Bf = bb + (50221746 + (293622 + (2350 + 22 * Math.Cos(bb) * Math.Cos(bb)) * Math.Cos(bb) * Math.Cos(bb)) * Math.Cos(bb) * Math.Cos(bb)) * Math.Sin(bb) * Math.Cos(bb) * Math.Pow(10, -10);tf = Math.Tan(Bf);nf = Math.Cos(Bf) * Math.Cos(Bf) * E1;Vf = Math.Pow(1 + E1, 0.5);Nf = a / Math.Pow(1 - E * Math.Sin(Bf) * Math.Sin(Bf), 0.5);B0 = Bf * 180 / Math.PI - 0.5 * Vf * Vf * tf * (yf * yf / Nf / Nf - 1 / 12.00 * (5 + 3 * tf * tf + nf - 9 * nf * tf * tf) * Math.Pow(yf / Nf, 4) + 1 / 360.00 * (61 + 90 * tf * tf + 45 * tf * tf) * Math.Pow(yf / Nf, 6)) * 180 / Math.PI;l0 = 1 / Math.Cos(Bf) * (yf / Nf - 1 / 6 * (1 + 2 * tf * tf + nf) * Math.Pow(yf / Nf, 3) + 1 / 120 * (5 + 28 * tf * tf + 24 * tf * tf + 6 * nf + 8 * nf * tf * tf) * Math.Pow(yf / Nf, 5)) * 180 / Math.PI;du = Convert.ToInt16( Math.Floor(B0));//dufen = Convert.ToInt16( Math.Floor((B0 - Math.Floor(B0)) * 60)); //fenmiao = Convert.ToInt16( Convert.ToInt16( ( (B0 - Math.Floor(B0)) * 60 - Math.Floor((B0 - Math.Floor(B0))* 60) ) * 60));//秒?textBox7.Text = Convert.ToString(du) + "度è"+ Convert.ToString(fen)+"分?"+Convert.ToString(miao)+"秒?";du =Convert.ToInt16(Math.Floor(l0));//dufen = Convert.ToInt16(Math.Floor((l0 - Math.Floor(l0)) * 60)); //fenmiao = Convert.ToInt16(Convert.ToInt16(((l0 - Math.Floor(l0)) * 60 - Math.Floor((l0 - Math.Floor(l0)) * 60)) * 60));//秒?textBox8.Text = Convert.ToString(du) + "度è"+ Convert.ToString(fen) + "分?"+ Convert.ToString(miao) + "秒?";private void toolStripMenuItem2_Click(object sender, EventArgs e)a = 6378245;E = 0.966;E1 = 0.4683;bbb = 6367588.4969;private void toolStripMenuItem3_Click(object sender, EventArgs e)a = 6378137;E = 0.90;E1 = 0.548;bbb = 6367452.133;private void toolStripMenuItem4_Click(objectsender, EventArgs e)a = 6378137;E = 0.13;E1 = 0.227;bbb = 6367452.133;//待鋣定¨private void toolStripMenuItem5_Click(object sender, EventArgs e)a = 6378140;E = 0.9588;E1 = 0.468;//待鋣定¨private void Form1_Load(object sender, EventArgs e)。

C++高斯正反算代码

CCoordinate CEarthPlaneDialog::EarthToPlane(const CEarth & earth){double B = earth.m_dB.m_realRadians;double L = earth.m_dL.m_realRadians;L = SetRad_PIPI(L);double L0 = SetL0(L);double X, N, t, t2, m, m2, ng2;double sinB, cosB;X = A1 * B * 180.0 / PI + A2 * sin(2 * B) + A3 * sin(4 * B) + A4 * sin(6 * B);sinB = sin(B);cosB = cos(B);t = tan(B);t2 = t * t;N = a / sqrt(1 - e2 * sinB * sinB);m = cosB * (L - L0);m2 = m * m;ng2 = cosB * cosB * e2 / (1 - e2);double x = X + N * t * ((0.5 + ((5 - t2 + 9 * ng2 + 4 * ng2 * ng2) / 24.0 + (61 - 58 * t2 + t2 * t2) * m2 / 720.0) * m2) * m2);double y = N * m * ( 1 + m2 * ( (1 - t2 + ng2) / 6.0 + m2 * ( 5 - 18 * t2 + t2 * t2 + 14 * ng2 - 58 * ng2 * t2 ) / 120.0));y += 500000;CCoordinate coor(x,y);return coor;}//将大地坐标系转换为平面直角坐标系CEarth CEarthPlaneDialog::PlaneToEarth(const CCoordinate & coor, UINT Sign){double x = coor.m_dX;double y = coor.m_dY;double L0 = (6.0 * Sign - 3.0) / 180 * PI;double sinB, cosB, t, t2, N ,ng2, V, yN;double preB0, B0;double eta;y -= 500000;B0 = x / A1;do{preB0 = B0;B0 = B0 * PI / 180.0;B0 = (x - (A2 * sin(2 * B0) + A3 * sin(4 * B0) + A4 * sin(6 * B0))) / A1;eta = fabs(B0 - preB0);}while(eta > 0.0000000000000001);B0 = B0 * PI / 180.0;sinB = sin(B0);cosB = cos(B0);t = tan(B0);t2 = t * t;N = a / sqrt(1 - e2 * sinB * sinB);ng2 = cosB * cosB * e2 / (1 - e2);V = sqrt(1 + ng2);yN = y / N;double B = B0 - (yN * yN - (5 + 3 * t2 + ng2 - 9 * ng2 * t2) * yN * yN * yN * yN / 12.0 + (61 + 90 * t2 + 45 * t2 * t2) * yN * yN * yN * yN * yN * yN / 360.0) * V * V * t / 2;double L = L0 + (yN - (1 + 2 * t2 + ng2) * yN * yN * yN / 6.0 + (5 + 28 * t2 + 24 * t2 * t2 + 6 * ng2 + 8 * ng2 * t2) * yN * yN * yN * yN * yN / 120.0) / cosB;CEarth earth(B,L,0,'r');return earth;}//将平面直角坐标系转换为大地坐标系void CEarthPlaneDialog::OnRADIOPlane1954(){// TODO: Add your control notification handler code herea = 6378245;f = 298.3;e2 = 1 - ((f - 1) / f) * ((f - 1) / f);e12 = (f / (f - 1)) * (f / (f - 1)) - 1;A1 = 111134.8611;A2 = -16036.4803;A3 = 16.8281;A4 = -0.0220;m_iReferenceFrame = 0;}void CEarthPlaneDialog::OnRADIOPlane1980(){// TODO: Add your control notification handler code herea = 6378140;f = 298.257;e2 = 1 - ((f - 1) / f) * ((f - 1) / f);e12 = (f / (f - 1)) * (f / (f - 1)) - 1;A1 = 111133.0047;A2 = -16038.5282;A3 = 16.8326;A4 = -0.0220;m_iReferenceFrame = 1;}void CEarthPlaneDialog::OnRADIOPlaneWGS84(){// TODO: Add your control notification handler code herea = 6378137;f = 298.257223563;e2 = 1 - ((f - 1) / f) * ((f - 1) / f);e12 = (f / (f - 1)) * (f / (f - 1)) - 1;A1 = 111133.0047;A2 = -16038.5282;A3 = 16.8326;A4 = -0.0220;m_iReferenceFrame = 2;}double CEarthPlaneDialog::SetL0(double L){L = L / PI * 180;if( L > 0 ){for(m_iCincture = 1 ; !(( L >= (6 * m_iCincture - 6 )) & ( L < ( 6 * m_iCincture ))); m_iCincture++){}L0 = 6 * m_iCincture - 3;}else if( L < 0 ){L = -L;for(m_iCincture = 1 ; !(( L > (6 * m_iCincture - 6 )) & ( L <= ( 6 * m_iCincture ))); m_iCincture++){}m_iCincture = 61 - m_iCincture;L0 = m_iCincture * 6 - 363;}else{L0 = 3;m_iCincture = 1;}double temp = ((double)L0) / 180 * PI;return temp;}double CEarthPlaneDialog::SetRad_PIPI(double Rad) {if(Rad > 0){double i = floor(Rad / PI / 2);Rad -= 2 * PI * i;}if(Rad < 0){Rad = -Rad;double i = floor(Rad / PI / 2);Rad -= 2 * PI * i;Rad = -Rad + PI * 2;}if(Rad > PI)Rad -= 2 * PI;return Rad;}。

高斯投影正反算程序

{
printf("高斯投影反算\nunder CGCS2000,六度带\n测试数据：x=2562083.2708,y=19512837.2851\n");
double a=6378137;
double f=1/298.257222101;
double e2=-f*f+2*f;
double x=2562083.2708;
double y=19512837.2851;
int Y=int(y);
int count=0;
while (Y>=100)
{
Y=Y/10;
count++;
}
y=y-Y*pow(10.0,count)-500000;
{
int model;
printf("高斯投影正算请输入1，反算请输入2\n");
scanf("%d",&model);
if(model==1)
{
printf("高斯投影正算\nunder WGS-84,六度带\n测试数据：B=22.155898294，L=111.285215387\n");
double W=sqrt(1-e2*sin(B)*sin(B));
double N=a/W;
double ep2=e2/(1-e2);
double ita2=ep2*cos(B)*cos(B);
double t=tan(B);
double temp1=5-t*t+9*ita2+4*ita2*ita2;
fBf=-a2*sin(2*Bf)/2+a4*sin(4*Bf)/4-a6*sin(6*Bf)/6+a8*sin(8*Bf)/8;

高斯投影正反算——包括3度和6度带的选择

// guass coordinateDlg.cpp : implementation file//#include "stdafx.h"#include "guass coordinate.h"#include "guass coordinateDlg.h"#include "math.h"#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif///////////////////////////////////////////////////////////////////////////// // CAboutDlg dialog used for App Aboutclass CAboutDlg : public CDialog{public:CAboutDlg();// Dialog Data//{{AFX_DATA(CAboutDlg)enum { IDD = IDD_ABOUTBOX };//}}AFX_DATA// ClassWizard generated virtual function overrides//{{AFX_VIRTUAL(CAboutDlg)protected:virtual void DoDataExchange(CDataExchange* pDX); // DDX/DDV support //}}AFX_VIRTUAL// Implementationprotected://{{AFX_MSG(CAboutDlg)//}}AFX_MSGDECLARE_MESSAGE_MAP()};CAboutDlg::CAboutDlg() : CDialog(CAboutDlg::IDD){//{{AFX_DATA_INIT(CAboutDlg)//}}AFX_DATA_INIT}void CAboutDlg::DoDataExchange(CDataExchange* pDX){CDialog::DoDataExchange(pDX);//{{AFX_DATA_MAP(CAboutDlg)//}}AFX_DATA_MAP}BEGIN_MESSAGE_MAP(CAboutDlg, CDialog)//{{AFX_MSG_MAP(CAboutDlg)// No message handlers//}}AFX_MSG_MAPEND_MESSAGE_MAP()///////////////////////////////////////////////////////////////////////////// // CGuasscoordinateDlg dialogCGuasscoordinateDlg::CGuasscoordinateDlg(CWnd* pParent /*=NULL*/) : CDialog(CGuasscoordinateDlg::IDD, pParent){//{{AFX_DATA_INIT(CGuasscoordinateDlg)m_bdu = 0;m_bfen = 0;m_bmiao = 0.0;m_x = 0.0;m_y = 0.0;m_ldu = 0;m_lfen = 0;m_lmiao = 0.0;m_ZoneWidth = -1;//}}AFX_DATA_INIT// Note that LoadIcon does not require a subsequent DestroyIcon in Win32 m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);}void CGuasscoordinateDlg::DoDataExchange(CDataExchange* pDX){CDialog::DoDataExchange(pDX);//{{AFX_DATA_MAP(CGuasscoordinateDlg)DDX_Text(pDX, IDC_B_DU, m_bdu);DDV_MinMaxInt(pDX, m_bdu, 0, 90);DDX_Text(pDX, IDC_B_FEN, m_bfen);DDV_MinMaxInt(pDX, m_bfen, 0, 60);DDX_Text(pDX, IDC_B_MIAO, m_bmiao);DDV_MinMaxDouble(pDX, m_bmiao, 0., 60.);DDX_Text(pDX, IDC_x, m_x);DDX_Text(pDX, IDC_y, m_y);DDX_Text(pDX, IDC_L_DU, m_ldu);DDV_MinMaxInt(pDX, m_ldu, 0, 180);DDX_Text(pDX, IDC_L_FEN, m_lfen);DDV_MinMaxInt(pDX, m_lfen, 0, 60);DDX_Text(pDX, IDC_L_MIAO, m_lmiao);DDV_MinMaxDouble(pDX, m_lmiao, 0., 60.);DDX_Radio(pDX, IDC_RADIO1, m_ZoneWidth);//}}AFX_DATA_MAP}BEGIN_MESSAGE_MAP(CGuasscoordinateDlg, CDialog)//{{AFX_MSG_MAP(CGuasscoordinateDlg)ON_WM_SYSCOMMAND()ON_WM_PAINT()ON_WM_QUERYDRAGICON()ON_COMMAND(ID_ABOUT, OnAbout)ON_COMMAND(ID_ZHENG, OnZheng)ON_COMMAND(ID_FAN, OnFan)ON_BN_CLICKED(ID_CAL, OnCal)ON_BN_CLICKED(IDC_RADIO1, OnRadio1)ON_BN_CLICKED(IDC_RADIO2, OnRadio2)//}}AFX_MSG_MAPEND_MESSAGE_MAP()///////////////////////////////////////////////////////////////////////////// // CGuasscoordinateDlg message handlersBOOL CGuasscoordinateDlg::OnInitDialog(){CDialog::OnInitDialog();// Add "About..." menu item to system menu.// IDM_ABOUTBOX must be in the system command range.ASSERT((IDM_ABOUTBOX & 0xFFF0) == IDM_ABOUTBOX);ASSERT(IDM_ABOUTBOX < 0xF000);CMenu* pSysMenu = GetSystemMenu(FALSE);if (pSysMenu != NULL){CString strAboutMenu;strAboutMenu.LoadString(IDS_ABOUTBOX);if (!strAboutMenu.IsEmpty()){pSysMenu->AppendMenu(MF_SEPARATOR);pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);}}// Set the icon for this dialog. The framework does this automatically // when the application's main window is not a dialogSetIcon(m_hIcon, TRUE); // Set big iconSetIcon(m_hIcon, FALSE); // Set small icon// TODO: Add extra initialization hereCMenu menu;menu.LoadMenu(IDR_MENU); //IDR_MENU在资源中创建的菜单SetMenu(&menu);return TRUE; // return TRUE unless you set the focus to a control}void CGuasscoordinateDlg::OnSysCommand(UINT nID, LPARAM lParam){if ((nID & 0xFFF0) == IDM_ABOUTBOX){CAboutDlg dlgAbout;dlgAbout.DoModal();}else{CDialog::OnSysCommand(nID, lParam);}}// If you add a minimize button to your dialog, you will need the code below // to draw the icon. For MFC applications using the document/view model, // this is automatically done for you by the framework.void CGuasscoordinateDlg::OnPaint(){if (IsIconic()){CPaintDC dc(this); // device context for paintingSendMessage(WM_ICONERASEBKGND, (WPARAM) dc.GetSafeHdc(), 0);// Center icon in client rectangleint cxIcon = GetSystemMetrics(SM_CXICON);int cyIcon = GetSystemMetrics(SM_CYICON);CRect rect;GetClientRect(&rect);int x = (rect.Width() - cxIcon + 1) / 2;int y = (rect.Height() - cyIcon + 1) / 2;// Draw the icondc.DrawIcon(x, y, m_hIcon);}else{CDialog::OnPaint();}}// The system calls this to obtain the cursor to display while the user drags // the minimized window.HCURSOR CGuasscoordinateDlg::OnQueryDragIcon(){return (HCURSOR) m_hIcon;}void CGuasscoordinateDlg::OnAbout(){// TODO: Add your command handler code hereCAboutDlg Dlg;Dlg.DoModal();}void CGuasscoordinateDlg::OnZheng(){// TODO: Add your command handler code hereGetDlgItem(IDC_B_DU)->EnableWindow(TRUE);GetDlgItem(IDC_B_FEN)->EnableWindow(TRUE);GetDlgItem(IDC_B_MIAO)->EnableWindow(TRUE);GetDlgItem(IDC_L_DU)->EnableWindow(TRUE);GetDlgItem(IDC_L_FEN)->EnableWindow(TRUE);GetDlgItem(IDC_L_MIAO)->EnableWindow(TRUE);GetDlgItem(IDC_x)->EnableWindow(FALSE);GetDlgItem(IDC_y)->EnableWindow(FALSE);SetDlgItemText(IDC_STATIC_INPUT,"请输入B,L坐标值：");m_x=m_y=0;//m_bdu=m_bfen=m_ldu=m_lfen=0;// m_bmiao=m_lmiao=0;UpdateData(FALSE);}void CGuasscoordinateDlg::OnFan(){// TODO: Add your command handler code hereGetDlgItem(IDC_x)->EnableWindow(TRUE);GetDlgItem(IDC_y)->EnableWindow(TRUE);GetDlgItem(IDC_B_DU)->EnableWindow(FALSE);GetDlgItem(IDC_B_FEN)->EnableWindow(FALSE);GetDlgItem(IDC_B_MIAO)->EnableWindow(FALSE);GetDlgItem(IDC_L_DU)->EnableWindow(FALSE);GetDlgItem(IDC_L_FEN)->EnableWindow(FALSE);GetDlgItem(IDC_L_MIAO)->EnableWindow(FALSE);SetDlgItemText(IDC_STATIC_INPUT,"请输入x,y坐标值：");m_bdu=m_bfen=m_ldu=m_lfen=0;m_bmiao=m_lmiao=0;// m_x=m_y=0;UpdateData(FALSE);}void CGuasscoordinateDlg::OnCal(){// TODO: Add your control notification handler code hereUpdateData(TRUE);double B,L,l,N,a0,a4,a6,a3,a5;int n,Lo;double b,Bf,Nf,b2,b3,b4,b5,Z;if(!GetDlgItem(IDC_x)->IsWindowEnabled()&&(!GetDlgItem(IDC_B_DU)->IsWindow Enabled())){MessageBox("请选择计算方式！","提示",MB_OK|MB_ICONWARNING);return;}if (m_ZoneWidth==-1){MessageBox("请选择投影分带方式！","提示",MB_OK|MB_ICONWARNING);return;}if(GetDlgItem(IDC_x)->IsWindowEnabled()){b=m_x/6367558.4969;n=int(m_y/1000000.0);Bf=b+(50221746+(293622+(2350+22*cos(b)*cos(b))*cos(b)*cos(b))*cos(b)*cos(b))*p ow(10,-10)*sin(b)*cos(b);Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf)) *cos(Bf)*cos(Bf);Z=(m_y-500000-n*1000000.0)/(Nf*cos(Bf));b2=(0.5+0.00336975*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b3=0.333333-(0.1666667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b4=0.25+(0.161612+0.005617*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b5=0.2-(0.16667-0.00878*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);B=(Bf-(1-(b4-0.147*Z*Z)*Z*Z)*Z*Z*b2)*p;l=(1-(b3-b5*Z*Z)*Z*Z)*Z*p;if (m_ZoneWidth==0)Lo=(n*6-3)*60*60;else if(m_ZoneWidth=1)Lo=n*3*60*60;L=Lo+l;m_bdu=int(B/60/60);m_bfen=int(B/60-m_bdu*60);m_bmiao=B-m_bdu*60*60-m_bfen*60;m_ldu=int(L/60/60);m_lfen=int(L/60-m_ldu*60);m_lmiao=L-m_ldu*60*60-m_lfen*60;}if(GetDlgItem(IDC_B_DU)->IsWindowEnabled()){B=(m_bdu*60*60+m_bfen*60+m_bmiao)/p;L=m_ldu*60*60+m_lfen*60+m_lmiao;if (m_ZoneWidth==0){n=int((m_ldu+m_lfen/60.0+m_lmiao/60.0/60.0)/6.0)+1;Lo=(6*n-3)*60*60;}else if (m_ZoneWidth==1){n=int((m_ldu+m_lfen/60.0+m_lmiao/60.0/60.0)/3.0);Lo=3*n*60*60;}l=(L-Lo)/p;N=6399698.902-(21562.267-(108.973-0.612*cos(B)*cos(B))*cos(B)*cos(B))*cos( B)*cos(B);a0=32140.404-(135.3302-(0.7092-0.0040*cos(B)*cos(B))*cos(B)*cos(B))*cos(B) *cos(B);a4=(0.25+0.00252*cos(B)*cos(B))*cos(B)*cos(B)-0.04166;a6=(0.166*cos(B)*cos(B)-0.084)*cos(B)*cos(B);a3=(0.3333333+0.001123*cos(B)*cos(B))*cos(B)*cos(B)-0.1666667;a5=0.0083-(0.1667-(0.1968+0.0040*cos(B)*cos(B))*cos(B)*cos(B))*cos(B)*cos( B);m_x=6367558.4969*B-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B)*cos(B);m_y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B)+500000.0+n*1000000.0;}UpdateData(FALSE);}void CGuasscoordinateDlg::OnRadio1(){// TODO: Add your control notification handler code herem_ZoneWidth=0;}void CGuasscoordinateDlg::OnRadio2(){// TODO: Add your control notification handler code herem_ZoneWidth=1;}。

一个老师给的高斯投影正、反算c++源码

一个老师给的高斯投影正、反算c++源码//高斯投影正、反算//////6度带宽 54年北京坐标系//高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres)void GaussProjCal(double longitude, double latitude, double *X, double *Y){int ProjNo=0; int ZoneWide; ////带宽double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval; double a,f, e2,ee, NN, T,C,A, M, iPI;iPI = 0.0174532925199433; ////3.1415926535898/180.0;ZoneWide = 6; ////6度带宽a=6378245.0; f=1.0/298.3; //54年北京坐标系参数////a=6378140.0; f=1/298.257; //80年西安坐标系参数ProjNo = (int)(longitude / ZoneWide) ;longitude0 = ProjNo * ZoneWide + ZoneWide / 2;longitude0 = longitude0 * iPI ;latitude0=0;longitude1 = longitude * iPI ; //经度转换为弧度latitude1 = latitude * iPI ; //纬度转换为弧度e2=2*f-f*f;ee=e2*(1.0-e2);NN=a/sqrt(1.0-e2*sin(latitude1)*sin(latitude1));T=tan(latitude1)*tan(latitude1);C=ee*cos(latitude1)*cos(latitude1);A=(longitude1-longitude0)*cos(latitude1);M=a*((1-e2/4-3*e2*e2/64-5*e2*e2*e2/256)*latitude1-(3*e2/8+3*e2*e2/32 +45*e2*e2*e2/1024)*sin(2*latitude1) +(15*e2*e2/256+45*e2*e2*e2/1024)*sin(4*latitude1)-(35*e2*e2*e2/3072) *sin(6*latitude1));xval = NN*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/120); yval = M+NN*tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24+(61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720);X0 = 1000000L*(ProjNo+1)+500000L;Y0 = 0;xval = xval+X0; yval = yval+Y0;*X = xval;*Y = yval;}//高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)void GaussProjInvCal(double X, double Y, double *longitude, double *l atitude){int ProjNo; int ZoneWide; ////带宽double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval; double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;iPI = 0.0174532925199433; ////3.1415926535898/180.0;a = 6378245.0; f = 1.0/298.3; //54年北京坐标系参数////a=6378140.0; f=1/298.257; //80年西安坐标系参数ZoneWide = 6; ////6度带宽ProjNo = (int)(X/1000000L) ; //查找带号longitude0 = (ProjNo-1) * ZoneWide + ZoneWide / 2;longitude0 = longitude0 * iPI ; //中央经线X0 = ProjNo*1000000L+500000L;Y0 = 0;xval = X-X0; yval = Y-Y0; //带内大地坐标e2 = 2*f-f*f;e1 = (1.0-sqrt(1-e2))/(1.0+sqrt(1-e2));ee = e2/(1-e2);M = yval;u = M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256));fai = u+(3*e1/2-27*e1*e1*e1/32)*sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1 /32)*sin(4*u)+(151*e1*e1*e1/96)*sin(6*u)+(1097*e1*e1*e1*e1/512)*sin(8 *u);C = ee*cos(fai)*cos(fai);T = tan(fai)*tan(fai);NN = a/sqrt(1.0-e2*sin(fai)*sin(fai));R = a*(1-e2)/sqrt((1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai))* (1-e2*sin(fai)*sin(fai)));D = xval/NN;//计算经度(Longitude) 纬度(Latitude)longitude1 = longitude0+(D-(1+2*T+C)*D*D*D/6+(5-2*C+28*T-3*C*C+8*ee+ 24*T*T)*D*D*D*D*D/120)/cos(fai);latitude1 = fai -(NN*tan(fai)/R)*(D*D/2-(5+3*T+10*C-4*C*C-9*ee)*D*D* D*D/24+(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720);//转换为度 DD*longitude = longitude1 / iPI;*latitude = latitude1 / iPI;}如果有需要程序的，可以直接跟我联系，呵呵附：高斯正反算参数pi=0.0174532925 ※※0.0174532925199433 //π长半轴a=6378245.0; 扁率f=1.0/298.3; //54年北京坐标系参数长半轴a=6378140.0; 扁率f=1/298.257; //80年西安坐标系参数长半轴a=6378137m；扁率f=1:298.257223563。

高斯投影的正算反算C++源代码-课程设计

高斯投影的正算反算C++源代码-课程设计高斯投影的正算反算C++源代码一、设计目的：加深理解高斯—克吕格投影的实质，掌握高斯—克吕格投影直角坐标的运用，理解通用坐标和自然坐标值的关系。

二、设计内容：1、编程实现由经纬度到直角坐标的转换；2、正确计算点所在的分度带；3、准确计算出点的通用坐标值。

三、方法与步骤：1、理解公式中每个字母的含义；2、编程。

3、编程的程序如下：// 高斯投影正反算公式.cpp : Defines the entry point for the console application.//#include "iostream.h"#include "math.h"int main(int argc, char* argv[]){int n,a,d,c;int L[4],B[4],C[3],D[3],L0;double l,N,a0,a4,a6,a3,a5,x,y,BB,P,b,Z,Bf,b2,b3,b4,b5,Nf,LL,xx,yy;cout<<"1----正算(L,B-->x,y) 2----反算(x,y-->L,B)"<<endl;cout<<"请选择正算（1）或反算（2）:"<<endl;cin>>n;/////////////////正算if(n==1){cout<<"请按顺序输入L,B:"<<endl;cout<<"输入L(度,分,秒):"<<endl;cin>>L[1]>>L[2]>>L[3];cout<<"请输入B(度,分,秒):"<<endl;cin>>B[1]>>B[2]>>B[3];BB=3.141592654*(B[1]+B[2]/60.0+B[3]/3600.0)/180;P=3.141592654*(B[1]*3600+B[2]*60+B[3])/(180*3600);////////////////计算点所在投影带的中央子午线L[0]=6*(int((L[1]+L[2]/60.0+L[3]/3600.0)/6)+1)-3;cout<<"点所在投影带的中央子午线为:"<<endl;cout<<L[0]<<"度"<<endl;l=3.141592654*(L[1]+L[2]/60.0+L[3]/3600-L[0])/180;N=6399698.902-(21562.267-(108.973-0.612*cos(BB)*cos(BB))*cos(BB)*co s(BB))*cos(BB)*cos(BB);a0=32140.404-(135.3302-(0.7092-0.0040*cos(BB)*cos(BB))*cos(BB)*cos( BB))*cos(BB)*cos(BB);a4=(0.25+0.00252*cos(BB)*cos(BB))*cos(BB)*cos(BB)-0.04166;a6=(0.166*cos(BB)*cos(BB)-0.084)*cos(BB)*cos(BB);a3=(0.3333333+0.001123*cos(BB)*cos(BB))*cos(BB)*cos(BB)-0.1666667; a5=0.0083-(0.1667-(0.1968+0.0040*cos(BB)*cos(BB))*cos(BB)*cos(BB))* cos(BB)*cos(BB);///////////代号计算a=L[0]/6+1;d=L[0]-1.5;c=d/3+1;///////////////计算坐标x=6367558.4969*P-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(BB)*cos(BB); y=(1+(a3+a5*l*l)*l*l)*l*N*cos(BB);/////////////格式控制cout.setf(ios::fixed);cout.setf(ios::showpoint);cout.precision(3);cout<<"正算结果为:"<<endl;cout<<"平面坐标:"<<endl;cout<<"x="<<x<<endl;cout<<"y="<<y<<endl;cout<<"分度带号为："<<endl;cout<<"六度带号为："<<a<<endl;cout<<"三度带号为："<<c<<endl;cout<<"通用坐标为（米）："<<endl;cout<<"x="<<x<<endl;cout<<"y="<<y+500000<<endl; //通用坐标计算}else if(n==2) /////////////////////反算{cout<<"请依次输入坐标x,y(m):"<<endl;cin>>xx>>yy;cout<<"请输入点所在投影带的中央子午线(度):"<<endl;cin>>L0;b=xx/6367558.4969;Bf=b+(50221746+(293622+(2350+22*cos(b)*cos(b))*cos(b)*cos(b))*cos(b)*cos(b))*(1.0E-10)*sin(b)*cos(b);Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*cos(Bf)*c os(Bf))*cos(Bf)*cos(Bf);Z=yy/(Nf*cos(Bf));b2=(0.5+0.00336975*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b3=0.333333-(0.1666667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf); b4=0.25+(0.161612+0.005617*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b5=0.2-(0.1667-0.00878*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);BB=Bf*3600*180/3.141592654-(1-(b4-0.147*Z*Z)*Z*Z)*Z*Z*b2*180*3600/3.141592654;LLC[1]=int((BB-C[0]*3600)/60);C[2]=int(BB-C[0]*3600-C[1]*60);cout<<"B="<<C[0]<<"度"<<C[1]<<"分"<<C[2]<<"秒"<<endl;D[0]=int(LL/3600);D[1]=int((LL-D[0]*3600)/60);D[2]=int(LL-D[0]*3600-D[1]*60);cout<<"L="<<D[0]<<"度"<<D[1]<<"分"<<D[2]<<"秒"<<endl;}return 0;}四、成果：1、程序运行结果如下图：1260。

基于Visual C++ 的高斯投影坐标反算

文章编号院1674-9146渊圆园15冤08原园059原园2
术 Applied Technology 应用技
基于 Visual C++的高斯投影坐标反算
许大亮
渊安徽理工大学测绘学院袁安徽淮南 232001冤
摘要院文中利用 Visual C+ + 语言编写高斯投影坐标反算的程序,采用的数据是西安 80 椭球参数和高斯坐标曰首先根据子午线弧长公式迭代计算求出垂足纬度袁其次算出大地纬度和经差袁最后可以得到大地坐标遥程序计算出的坐标精度为 0.000 1义遥结果表明院用转换程序解算出的坐标存在一定的误差袁但是可以满足低精度的工程建设遥关键词院Visual C++语言曰坐标反算曰高斯坐标曰垂足纬度中图分类号院TP311 文献标志码院A DOI院10.3969/j.issn.1674-9146.2015.08.059
根据式渊3冤可写出迭代式
Bi+1 = X/c 茁0 -
蓸蔀 3
5
7
茁2cosBi + 茁4cos Bi + 茁6cos Bi + 茁8cos Bi
sinBi 茁0
.
渊4冤
1.3 Mf袁 Nf袁 tf袁浊f等关于Bf 的公式
地球椭球的长半径为 a袁短半径为 b遥
第一偏心率
e
=
姨a2
-
2
b
a.
换算后坐标的精度遥在垂足纬度的计算中袁常采用
迭代计算的方法遥笔者参考了一些文献尧著作袁利
用 Visual C++语言设计了一个程序袁实现了高斯投
影坐标转换成西安 80 大地坐标袁计算结果精确可

高斯投影正反算编程一.高斯投影正反算基本公式

高斯投影正反算编程一．高斯投影正反算基本公式（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].L0)!=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)/720;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].L0)!=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*sin(8*B1 )/8;B2=(g[i].x-F)/a0;B3=B1;B1=B2;} while(fabs(B3-B2)>10e-10);//利用迭代算法求解底点纬度//return B2; }。

大地测量_高斯正反算程序实现_C#

界面如下：输入数据计算：C#代码如下：using System;using System.Collections.Generic; using ponentModel;using System.Data;using System.Drawing;using System.Text;using System.Windows.Forms;namespace大地{public partial class Form1 : Form {double B, L;double x, y;double X, Y;double L0, l1;double p=206264.80625;public Form1(){InitializeComponent();}private void Form1_Load(object sender, EventArgs e){//自动化控件显示初始值radioButton2.Checked=true;radioButton3.Checked = true;textBox1.Focus();textBox1.Text = textBox2.Text = textBox3.Text = textBox4.Text = textBox5.Text = textBox6.Text = "0";this.richTextBox2.Text = "说明：输入的坐标需为按6°带投影且采用克氏椭球参数所得的国家统一坐标";}private void button1_Click(object sender, EventArgs e){//获取输入数据double bb1=Convert.ToDouble(this.textBox1.Text);double bb2=Convert.ToDouble(this.textBox2.Text);double bb3=Convert.ToDouble(this.textBox3.Text);double ll1=Convert.ToDouble(this.textBox4.Text);double ll2=Convert.ToDouble(this.textBox5.Text);double ll3=Convert.ToDouble(this.textBox6.Text);//检查输入格式的正确性if (bb1 >= 0 && bb1 <90 && bb2 >= 0 && bb2 < 60 && bb3 >= 0 && bb3 < 60){B = bb1 * 3600 + bb2 * 60 + bb3;}else{MessageBox.Show("纬度输入格式不正确！", "警告");return;}if (ll1 >= 0 && ll1<360 && ll2 >= 0 && ll2 < 60 && ll3 >= 0 && ll3 < 60)L = ll1 * 3600 + ll2 * 60 + ll3;}else{MessageBox.Show("经度输入格式不正确！", "警告");return;}double b = B/p;//获取用户选项值int tyd=radioButton1.Checked?3:6;int tq=radioButton3.Checked?1:2;//按带求带号以及中央子午线经度if (tyd == 6){N =(int)( L / (6*3600.0)) + 1;L0 = (6 * N - 3)*3600.0;l1 = (L-L0)/p;}if(tyd==3){N =(int)(L / (3*3600.0)+0.5);L0 = 3 * N * 3600.0;l1 =( L-L0)/p;}//采用克氏椭球参数的计算公式if (tq == 1){double c = Math.Pow( Math.Cos(b),2);double c1= Math.Sin(b) * Math.Cos(b);double c2=Math.Cos(b);double l2 = Math.Pow(l1,2);double n = 6399698.902 - (21562.267 - (108.973 - 0.612 * c) * c) * c;double a0 = 32140.404 - (135.3302 - (0.7092 - 0.0040 * c) * c) * c;double a4 = (0.25 + 0.00252 * c) * c - 0.04166;double a6 = (0.166 * c - 0.084) * c;double a3 = (0.3333333 + 0.001123 * c) * c - 0.1666667;double a5 = 0.0083 - (0.1667 - (0.1968 + 0.0040 * c) * c) * c; x = 6367558.4969 * b - (a0 - (0.5 + (a4 + a6 * l2) * l2) * l2 * n)*c1;y = (1 + (a3 + a5 * l2) * l2) * l1 * n * c2;X = x;double y1 = y + 500000.0;for (int i = 1; y1 / i > 1; i = i * 10){Y = N * i*10 + y1;}string tuoqiu = "采用克氏椭球参数,";this.richTextBox1.Text = "按经差" + tyd + "°进行投影分带，" + tuoqiu + "其计算结果为：" + "\r\n" + "x=" + x + "\r\n" + "y=" + y + "\r\n" +"国家统一坐标为：" + "\r\n" + "X=" + X + "\r\n" + "Y=" + Y; }//采用1975国际椭球参数的计算公式if (tq == 2){double c = Math.Cos(b) * Math.Cos(b);double c1 = Math.Sin(b) * Math.Cos(b);double c2 = Math.Cos(b);double l2 = Math.Pow(l1, 2);double n = 6399596.652 - (21565.045 - (108.996 - 0.603 * c) * c) * c;double a0 = 32144.5189 - (135.3646 - (0.7034 - 0.0041 * c) * c) * c;double a4 = (0.25 + 0.00253 * c) * c - 0.04167;double a6 = (0.167 * c - 0.083) * c;double a3 = (0.3333333 + 0.001123 * c) * c - 0.1666667;double a5 = 0.00878 - (0.1702 - 0.20382 * c) * c;x = 6367452.1328 * b-(a0 - (0.5 + (a4 + a6 * l2) * l2) * l2 * n) * c1;y = (1 + (a3 + a5 * l2) * l2) * l1 * n * c2;X = x;double y1 = y + 500000.0;for (int i = 1; y1 / i > 1; i = i * 10)//Y坐标加代号{Y = N * i * 10 + y1;}string tuoqiu = "采用1975国际椭球参数,";this.richTextBox1.Text = "按经差" + tyd + "°进行投影分带，" + tuoqiu + "其计算结果为：" + "\r\n" + "x=" + x + "\r\n" + "y=" + y+"\r\n"+"国家统一坐标为：" + "\r\n" + "X=" + X + "\r\n"+"Y="+Y;}}private void button2_Click(object sender, EventArgs e){X = Convert.ToDouble(this.textBox7.Text);//获取输入数据Y = Convert.ToDouble(this.textBox8.Text);for (int i = 1; Y/i >=10; i = i * 10)//对Y坐标处理并求出中央子午线经度{y = Y - (int)(Y / i) * i-500000;L0 = 6 * (int)(Y / i) - 3;}//按6°带克氏椭球反算double bt = x / 6367558.4969*p;double BT = x / 6367558.4969;double c3=Math.Cos(BT)*Math.Cos(BT);double c4=Math.Sin(BT)*Math.Cos(BT);doubleBf=(bt+(50221746+(293622+(2350+22*c3)*c3)*c3)*c4*Math.Pow(10,-10)*p)/p;double c5=Math.Pow(Math.Cos(Bf),2);double c6=Math.Sin(Bf)*Math.Cos(Bf);double Nf=6399698.902-(21562.267-(108.973-0.612*c5)*c5)*c5;double Z=y/(Nf*Math.Cos(Bf));double b2 = (0.5 + 0.003369 * c5) * c6;double b3 = 0.333333 - (0.166667 - 0.001123 * c5) * c5;double b4 = 0.25 + (0.16161 + 0.00562 * c5) * c5;double b5=0.2-(0.1667-0.0088*c5)*c5;double z2=Math.Pow(Z,2);B = (Bf*p - (1 - (b4 - 0.12 *z2) * z2) * z2 * b2 * p)/3600.0;L = L0+((1 - (b3 - b5 * z2) * z2) * Z * p)/3600.0;this.richTextBox2.Text = "按6°带克氏椭球反算后，结果为："+"\r\n"+"B=" + B + "\r\n" + "L=" + L;} }}。