公式省略了相邻单元格高斯投影正反算编程

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

（1）高斯正算基本公式
（2）高斯反算基本公式


* *





以上主要通过大地测量学基础课程得到，这不进行详细的推
导，只是列出基本公式指导编程的进 行。
二．编程的基本方法和流程图
（1）编程的基本方法
高斯投影正反算基本上 运用了所有的编程基本语句，本
文中是利用C++语言进行基本的设计。高斯正算中对椭球参
数 和带宽的选择主要运用了选择语句。而高斯反算中除了选
择语句的应用，在利用迭代算法求底点纬度还应 用了循环语
句。编程中还应特别注意相关的度分秒和弧度之间的相互转
* *
换，这是极其重要的。
（2）相关流程图
1）正算

* *


开始


输入大地坐标B，L
和经差L0
选择椭球参数
选择带宽
36度带
3度带 6度带
计算带号

计算带号


计算弧长

计算弧长

计算平面坐标x,y

计算平面坐标x,y

打印x,y
打印x,y

* *
2）反算

开始



输入自然值坐标x,y
和经差L0
选择椭球参数
利用迭代算法
求解底点纬度
利用公式计算B
和L
打印B和L
* *
三．编程的相关代码
（1）正算
# include
# include
# include
# include
#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(
assert(r!=NULL);
FILE *w=fopen(
assert(r!=NULL);
* *


int i=0;
while(fscanf(r,
!=EOF)
{
double a,b;
int zuobiao;
printf(请输入坐标系：北京54=1，西安80=2，
WGS84=3：





scanf(
getchar();
if(zuobiao==1)
{
a=6378245;
b=6356863.;


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((ab)*(ab)-1);求椭球的第一，第二曲率

double m0,m2,m4,m6,m8;


double a0,a2,a4,a6,a8;
m0=a*(1-e*e);
m2=3*e*e*m02;








m4=5*e*e*m24;
m6=7*e*e*m46;
m8=9*e*e*m68;
a0=m0+m22+3*m48+5*m616+35*m8128;
a2=m22+m42+15*m632+7*m816;
a4=m48+3*m616+7*m832;
a6=m632+m816;
a8=m8128;
* *




double Bmiao,Lmiao, L0miao;

Bmiao=(int)(g[i].B)* 3600.0+(int)((g[i].B-(int)(g[i].B))*1
00.0)*60. 0+(g[i].B*100-(int)(g[i].B*100))*100.0;
< br>Lmiao=(int)(g[i].L)*3600.0+(int)((g[i].L-(int)( g[i].L))*10
0.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=pi180.03600.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=asqrt(1-e*e*sin(B1)*sin(B1));
int Zonewide;
int Zonenumber;
printf(请输入带宽：3度带或6度带Zonewide=






scanf(
getchar();

if(Zonewide==3)
{

Zonenumber=(int )((g[i].L-Zonewide2)Zonewide+1);









}
else if(Zonewide==6)
{
Zonenumber=(int)g[i].LZonewide+1;
}
else
{
printf(错误
exit(0);
* *



}选择带宽
double
FE=Zonenumber*1000000+500000;改写为国家
通用坐标
< br>y=FE+N*A+A*A*A*N*(1-T*T+n*n)6+A*A*A*A*A*N*(5-1< br>8*T*T+T*T*T*T+14*n*n-58*n*n*T*T)120;

x=X+tan(B1)*N*A*A2+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< br>+T*T*T*T)720;




printf(所选坐标系的转换结果：x=%lf y=%lfn


}
fprintf(w,输出结果到文本文件


fclose(r);
fclose(w);



}
system(
return 0;
* *

（2）反算
# include
# include
# include
# include
#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(
assert(r!=NULL);
FILE *w=fopen(
assert(r!=NULL);


int i=0;
while(fscanf(r,%lf %lf
!=EOF)文件为空，无法打开
{
double a1=6378245.;克拉索夫斯基椭
球参数
double b1=6356863.;
double a75=6378140.;1975国际椭
球参数


double b75=6356755.2881575287;
double a84=6378137.;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(选择参考椭球：1=克拉索夫斯基椭球，2=1975
国际椭球，3=WGS-84系椭 球:
scanf(




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=aq;
double H,I;
A=g[i].yN;






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*pi180.0+AG-H(6*G)*(1.0+2*B+C)+I(12
* *
0*G)*(5.0+28*B+24*B*B+6*C+8*B*C);利用公式求解经
纬度








int Bdu,Bfen,Ldu,Lfen;
double Bmiao,Lmiao;
Ldu=int(LL0pi*180);
Lfen=int((LL0pi*180)*60-Ldu*60);
Lmiao=LL0pi*180*3600-Ldu*3600-Lfen*60;
Bdu=int(BB0pi*180);
Bfen=int((BB0pi*180)*60-Bdu*60);
Bmiao=BB0pi*180*3600-Bdu*3600-Bfen*60;
将弧度转化为角度
printf(所选坐标系的转换结果：%d度%d分%lf
秒 %d度%d分%lf秒
n
fprintf(w,°%d’%lf”%d°%d’%lf”
n将结果输出到文本
文件



}


fclose(r);
fclose(w);
* *


system(
return 0;
}
double sort(double a,double b)
{
double e,e1;
e=sqrt(1-(ba)*(ba));
e1=sqrt((ab)*(ab)-1);



double m0,m2,m4,m6,m8;
double a0,a2,a4,a6,a8;
m0=a*(1-e*e);
m2=3*e*e*m02;








m4=5*e*e*m24;
m6=7*e*e*m46;
m8=9*e*e*m68;
a0=m0+m22+3*m48+5*m616+35*m8128;
a2=m22+m42+15*m632+7*m816;
a4=m48+3*m616+7*m832;
a6=m632+m816;
a8=m8128;
B1=g[i].xa0;
do
* *
{

F=-a2*sin(2*B1)2+a4*sin(4* B1)4-a6*sin(6*B1)6+
a8*sin(8*B1)8;




} while(fabs(B3-B2)>10e-10);利用迭代算法求解底点纬度

return B2;

}
B2=(g[i].x-F)a0;
B3=B1;
B1=B2;