高斯投影正反算及换带计算VB程序设计

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

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

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

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

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

（2）相关流程图1）正算选择带宽 3/6 度带计算带号输入大地坐标 B ，L和经差 L06 度带3 度带选择椭球参数计算带号计算弧长计算平面坐标 x,y打印 x,y开始计算平面坐标 x,y计算弧长打印 x,y2）反算开始输入自然值坐标x,y和经差L0选择椭球参数利用迭代算法求解底点纬度利用公式计算B和L打印B 和L三．编程的相关代码（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;}。

大地测量高斯投影正反算程序代码课程：大地测量学基础姓名：林江伟学号：2008301610045班级： 0804 界面如下：输入数据计算：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;int N;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);double Bf=(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;}}。

高斯投影正、反算及换带程序执行条件※数组投影选择T、换算点个数“Z=0 F≠0”、=0正算0、≠0反算※坐标系选择“54 ≠54”、=54换算为1954年北京坐标系输入54、≠54换算为1988年西安坐标系M、中央子午线经度（°′″）输入※大地坐标I、序列号B、L：大地纬度和经度（地理坐标）（°′″）※高斯平面坐标轴子午线I、序列号X、Y：高斯平面坐标（m） Z、轴子午线（°）输出※大地坐标子午收敛角N、序列号B、L：大地纬度和经度（地理坐标）（°′″） R、子午收敛角（°′″）※高斯平面坐标子午收敛角N、序列号X、Y：高斯平面坐标（m） R、子午收敛角（°′″）注：1、程序执行前必须进行数组定位。

如：Defm 10 T×2=5×2=102、Y坐标值要去掉带号及避免出现负值的500公里；4、本程序运算时，各已知数据、观测变量不会随之变化，可非常方便地进行各数据的核对;5、本程序在进行换带计算时采用的是间接换带计算法。

Prog GSXYDefm 10:TA“Z=0 F≠0”G“54 ≠54”Z:Fixm:I=0:「b」0:I=I+1◢J=2I-1:M=Z[J:L=Z[J+1:A=0=>Prog“3”:B=M:M=L+Z:Prog“3”:L=M:{BL}:M=B:Prog“2”: B=M:M=L:Prog“2”:L=M-Z:≠>X=M:Y=L:{XY}:B=X:L=Y⊿Z[J]=B:Z[J+1]=L:I<T=>Goto 0⊿G=54=>C=6399698.90178271:E=.006738525414684:≠>C=6399596.65198801:E=.006 739501819473⊿I=0:「b」0:I“N”=I+1◢J=2I-1:B=Z[J:L=Z[J+1:A≠0=>X=B:Y=L:Goto 2⊿S=sin B:G=54=>F=111134.8611B-(32 005.7799S+133.9238S∧3+.6973S∧5+.0039S∧7)cos B:≠>F=111133.0047B-(32009.857 S+133.9602S∧3+.6976S∧5+.0039S∧7)cos B⊿U=√Ecos B:V=√(1+U2:N=C÷V:W=tan B: M=cos B(Lπ÷180:X=F+NW(.5M2+1┛24(5-W2+9U2+4U∧4)M∧4+1┛720(61-58W2+W∧4)M∧6◢Y=N(M+1┛6(1-W 2+U 2)M ∧3+1┛120(5-18W 2+W ∧4+14U 2-58U 2W 2)M ∧5◢M=W ┛π(180M+60(1+3U 2+2U ∧4)M ∧3+12(2-W 2)M ∧5:Goto 3:「b 」2:W=E ﹣6X-3:G=54=>F=27.11115372595+9.024********W-.00579740442W 2-4.3532572E ﹣4W ∧3+4.857285E ﹣5W ∧4+2.15727E ﹣6W ∧5-1.9399E ﹣7W ∧6:≠>F=27.11162289465+9.024********W-.00579850656W2-4.3540029E ﹣4W ∧3+4.858357E ﹣5W ∧4+2.15769E ﹣6W ∧5-1.9404E ﹣7W ∧6⊿U=√Ecos F:V=√(1+U 2:Q=YV ÷C:W=tan F:M=F-(1+U 2)W ┛π(90Q 2-7.5(5+3W 2+U 2-9U 2W 2)Q ∧4+.25(61+90W 2+45W ∧4)Q ∧6:Prog “3”:B=M ◢M=Z+1┛(πcos F)(180Q-30(1+2W 2+U 2)Q ∧3+1.5(5+28W 2+24W ∧4)Q ∧5:Prog “3”:L=M ◢M=W ┛π(180Q-60(1+W 2-U 2)Q ∧3+12(2+5W 2+3W ∧4)Q ∧5:「b 」3:Prog “3”:R=M ◢ I<T=>Goto 1⊿“END ”概要说明：我国的经度范围西边自73°起，东边至135°，可分成6°带共11带或3°共22带。

#include "stdafx.h"#include "iostream.h"#include "math.h"#include "stdio.h"#define P 206264.806247096355#define PI 3.141592653589793void GaosZ_fun(){printf("高斯投影的正算\n");double l,L,B,n2,x,y,N,t,V,c,e2;double i,j,k,n,h,a0,a4,a6,a3,a5,cB2;int m;e2=0.006738525414683;c=6399698.901782711;B=17.33557339*3600/P;L=119.15521159*3600/P;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;a5=0.0083-(0.1667-(0.1968+0.0040*cB2)*cB2)*cB2;// x=X+N*sin(B)*cos(B)*l*l/2+N*sin(B)*pow(cos(B),3)*(5-t*t+9*n2+4*n2*n2)*pow(l, 4)/24+N*sin(B)*pow(cos(B),5)*(61-58*t*t+pow(t,4))*pow(l,6)/720;// y=N*cos(B)*l+N*pow(cos(B),3)*(1-t*t+n2)*pow(l,5)/6+N*pow(cos(B),5)*(5-18*t*t +pow(t,4)+14*n2-58*n2*t*t)*pow(l,5)/120;x=6367558.4969*B-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B)*cos(B);y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B);printf("x=%f\ny=%f\n",x,y);}void GaosF_fun(){printf("高斯投影的反算\n");double B,Bf,Nf,b,b2,b3,b4,b5,Z,x,y,L0,l;// printf("输入x y 的值\n");// scanf("x=%f",&x);// scanf("y=%f",&y);x=3380330.875;y=320089.976;L0=111;b=x/6367558.4969;Bf=b+(50221746+ ( 293622+ (2350+22*cos(b)*cos(b))*cos(b)*cos(b) )*cos(b)*cos(b)) *sin(b)*cos(b)*1e-10;Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf))*cos( Bf)*cos(Bf);Z=y/Nf/cos(Bf);b2=(0.5+0.003369*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b3=0.333333-(0.166667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b5=0.2-(0.167-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);B=Bf-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2;l=(1-(b3-b5*Z*Z)*Z*Z)*Z;printf("B=%f\nL=%f\n",B*180/PI,(L0+l*180/PI));}void GaosLT_fun(){double B,Bf,Nf,b,b2,b3,b4,b5,Z,x,y,L01,L02,n2,l,V,L,N,c,t,e2;double a0,a4,a6,a3,a5,cB2;e2=0.006738525414683;c=6399698.901782711;// printf("输入x y 的值\n");// cout<<"x=";// cin>>x;// cout<<"y=";// cin>>y;// scanf("x=%f",&x); ///////////// ?????// scanf("y=%f",&y);// x=3380330.875;// y=320089.976;x=1944359.607;y=240455.4563;L01=117*3600/P;L02=120*3600/P;b=x/6367558.4969;Bf=b+(50221746+ ( 293622+ (2350+22*cos(b)*cos(b))*cos(b)*cos(b) )*cos(b)*cos(b)) *sin(b)*cos(b)*1e-10;Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf))*cos( Bf)*cos(Bf);Z=y/Nf/cos(Bf);b2=(0.5+0.003369*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b3=0.333333-(0.166667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b5=0.2-(0.167-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);l=(1-(b3-b5*Z*Z)*Z*Z)*Z;L=l+L01; ///反算就出L Bl=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;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;a5=0.0083-(0.1667-(0.1968+0.0040*cB2)*cB2)*cB2;// x=X+N*sin(B)*cos(B)*l*l/2+N*sin(B)*pow(cos(B),3)*(5-t*t+9*n2+4*n2*n2)*pow(l, 4)/24+N*sin(B)*pow(cos(B),5)*(61-58*t*t+pow(t,4))*pow(l,6)/720;// y=N*cos(B)*l+N*pow(cos(B),3)*(1-t*t+n2)*pow(l,5)/6+N*pow(cos(B),5)*(5-18*t*t +pow(t,4)+14*n2-58*n2*t*t)*pow(l,5)/120;x=6367558.4969*B-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B)*cos(B);y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B);printf("输出高斯邻带坐标换算\n");printf("x=%f\ny=%f\n",x,y);}void main(){int i;printf("输入1 2 3 选择投影的计算方法\n"); scanf("%d",&i);switch (i){case 1: GaosZ_fun(); break;case 2: GaosF_fun(); break;case 3: GaosLT_fun(); break;default : break;}}。