高斯正反算程序代码--C语言
luyued 发布于 2011-06-13 18:48 浏览 N 次
#define PI 3.141592653589793
double DMS2RAD(double dmsAngle){int degAngle,minAngle,intSignOfDms;double radAngle,secAngle;
intSignOfDms = 1;if(dmsAngle<0) intSignOfDms=-1; dmsAngle=fabs(dmsAngle);
degAngle=(int)(dmsAngle+0.0001);minAngle=(int)((dmsAngle-degAngle)*100+0.0001);secAngle=(dmsAngle-degAngle-minAngle/100.0)*10000.0; radAngle=(degAngle+minAngle/60.0+secAngle/3600.0)*PI/180.0; radAngle=radAngle*intSignOfDms; return radAngle;}
double RAD2DMS(double radAngle){ int degAngle,minAngle,intSignOfRad;double secAngle,dmsAngle;
intSignOfRad = 1;if(radAngle<0) intSignOfRad=-1; radAngle=fabs(radAngle);
secAngle=radAngle*180.0/PI*3600.0; degAngle=(int)(secAngle/3600+0.0001); minAngle=(int)((secAngle-degAngle*3600.0)/60.0+0.0001); secAngle=secAngle-degAngle*3600.0-minAngle*60.0;dmsAngle=degAngle+minAngle/100.0+secAngle/10000.0; dmsAngle=dmsAngle*intSignOfRad;return dmsAngle;}
void a0a2a4a6a8(double a,double e,double *Coeficient_a0){ double m0,m2,m4,m6,m8; 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;
Coeficient_a0[0]=m0+m2/2+3*m4/8+5*m6/16+35*m8/128; Coeficient_a0[1]=m2/2+m4/2+15*m6/32+7*m8/16; Coeficient_a0[2]=m4/8+3*m6/16+7*m8/32; Coeficient_a0[3]=m6/32+m8/16; Coeficient_a0[4]=m8/128;}void main(){ int h,k; double a,Alfa; double dmslat,dmslon,dmsl0; double a0 ,a2 ,a4, a6, a8; double radlat,radlon,radl0,l; double b,t,sb,cb,ita,e1,e; double X,l0; double N,c,v; double coor_x,coor_y; double Bf0,Bf1,dB,FBf,bf; double itaf, tf; double Nf,Mf; double B,L,dietaB,dietal; double sinBf,cosBf; double *Coeficient_a0; Coeficient_a0=(double *)malloc(5*sizeof(double)); printf("正算请选择1， 反算请选择2：\n"); scanf("%d",&k); if(k==1) //正算 { printf("请选择坐标系:\n"); printf(" 选择WGS-84坐标系,请按1\n"); printf(" 选择BJ-54坐标系,请按2\n"); printf(" 选择GDZ-80坐标系,请按3\n"); printf(" 其他坐标系,请按4\n"); scanf("%d",&h); if(h==1) a=6378137,Alfa=1.0/298.257223563; if(h==2) a=6378245,Alfa=1.0/298.3; if(h==3) a=6378140,Alfa=1.0/298.257;if(h==4) {printf("输入椭球长轴：");scanf("%lf",&a);printf("输入椭球扁率分母："); scanf("%lf",&Alfa);Alfa=1.0/Alfa;} printf("请输入已知点纬度:\n"); scanf("%lf",&dmslat); printf("请输入已知点经度:\n"); scanf("%lf",&dmslon); printf("请输入中央子午线经度 :\n"); scanf("%lf",&dmsl0); radlat=DMS2RAD(dmslat); radlon=DMS2RAD(dmslon); radl0=DMS2RAD(dmsl0); l=radlon-radl0;
b=a*(1-Alfa); sb=sin(radlat); cb=cos(radlat); t=sb/cb; e1=sqrt((a/b)*(a/b)-1); e=sqrt(1-(b/a)*(b/a)); ita=e1*cb; a0a2a4a6a8(a,e,Coeficient_a0); a0=Coeficient_a0[0]; a2=Coeficient_a0[1]; a4=Coeficient_a0[2]; a6=Coeficient_a0[3]; a8=Coeficient_a0[4];
X=a0*radlat-sb*cb*((a2-a4+a6)+(2*a4-16*a6/3)*sb*sb+16*a6*sb*sb*sb*sb/3.0); c=a*a/b; v=sqrt(1+e1*e1*cb*cb); N=c/v;
coor_x=X+N*sb*cb*l*l/2+N*sb*cb*cb*cb*(5-t*t+9*ita*ita+4*ita*ita*ita*ita)*l*l*l*l/24+N*sb*cb*cb*cb*cb*cb*(61-58*t*t+t*t*t*t)*l*l*l*l*l*l/720; coor_y=N*cb*l+N*cb*cb*cb*(1-t*t+ita*ita)*l*l*l/6+N*cb*cb*cb*cb*cb*(5-18*t*t+t*t*t*t+14*ita*ita-58*ita*ita*t*t)*l*l*l*l*l/120;
printf("coor_x=%.4lf\n",coor_x); printf("coor_y=%.4lf\n",coor_y);
}
if(k==2)//反算 { printf("请选择坐标系:\n"); printf(" 选择WGS-84坐标系,请按1\n"); printf(" 选择BJ-54坐标系,请按2\n"); printf(" 选择GDZ-80坐标系,请按3\n"); printf(" 其他坐标系,请按4\n"); scanf("%d",&h); if(h==1) a=6378137,Alfa=1.0/298.257223563; if(h==2) a=6378245,Alfa=1.0/298.3; if(h==3) a=6378140,Alfa=1.0/298.257;if(h==4) {scanf("输入椭球长轴：%lf",&a); scanf("输入椭球扁率分母：%lf",&Alfa);Alfa=1.0/Alfa;} printf("请输入l0:\n"); scanf("%lf",&l0); l0=l0*PI/180.0; printf("请输入x坐标:\n"); scanf("%lf",&coor_x); printf("请输入y坐标:\n"); scanf("%lf",&coor_y);
b=a*(1-Alfa); e1=sqrt((a/b)*(a/b)-1); e=sqrt(1-(b/a)*(b/a));
a0a2a4a6a8(a,e,Coeficient_a0); a0=Coeficient_a0[0]; a2=Coeficient_a0[1]; a4=Coeficient_a0[2]; a6=Coeficient_a0[3]; a8=Coeficient_a0[4];
X=coor_x; Bf0=X/a0; do{ sinBf=sin(Bf0);cosBf=cos(Bf0); FBf=-sinBf*cosBf*((a2-a4+a6)+(2.0*a4-16.0*a6/3.0)*sinBf*sinBf+(16.0/3.0)*a6*(sinBf*sinBf*sinBf*sinBf)); Bf1=(X-FBf)/a0; dB=Bf1-Bf0; Bf0=Bf1; } while(fabs(dB*180.0/PI*3600.0)>0.00001); bf=Bf1; c=a*a/b; v=sqrt(1+e1*e1*cos(bf)*cos(bf)); Nf=c/v; Mf=c/(v*v*v); tf=sin(bf)/cos(bf); itaf=e1*cos(bf);
dietaB=tf*coor_y*coor_y/(2*Mf*Nf)-tf*(5+3*tf*tf+itaf*itaf-9*tf*tf*itaf*itaf)*coor_y*coor_y*coor_y*coor_y/(24*Mf*Nf*Nf*Nf)+(61+90*tf*tf+45*tf*tf*tf*tf)*coor_y*coor_y*coor_y*coor_y*coor_y*coor_y/(720*Mf*Nf*Nf*Nf*Nf*Nf); dietal=coor_y/(Nf*cos(bf)+(1+2*tf*tf+itaf*itaf)*cos(bf)*coor_y*coor_y/(6*Nf))+(5+44*tf*tf+32*tf*tf*tf*tf-2*itaf*itaf-16*itaf*itaf*tf*tf)/(360*Nf*Nf*Nf*Mf*Mf*cos(bf));
B=bf-dietaB; L=l0+dietal; dmslat=RAD2DMS(B); dmslon=RAD2DMS(L); printf("已知点的纬度(ddd.mmss)为：.9lf\n",dmslat); printf("已知点的经度(ddd.mmss)为：.9lf\n",dmslon);
}}
上一篇：4月24日 ALFA 下一篇：朗文斯汀
相关资讯
- 07-01· 禁教唐诗算术能还幼儿快
- 07-01· 2011年06月17日
- 07-01· 唐诗宋词英译：李商隐 筹
- 07-01· 仿评《唐诗1000首》第186首
- 07-01· 没事干的时候背背唐诗吧
- 07-01· [转载]唐诗中“斜”字该读
- 07-01· 湖南醴陵瓷业转型升级
- 07-01· 奇瑞风云2两厢黑色|2010款
- 07-01· 摩根士丹利华鑫摩根士丹
- 07-01· 摩根士丹利华鑫近期优选
图文资讯
最新资讯
- 07-01· 中金投行部大摩出售中金
- 07-01· 摩根士丹利招聘6月2日【实
- 07-01· 营养防病圣典
- 07-01· 《博伽梵歌原意》之第十
- 07-01· [不错]斑斓圣典---减肥中常
- 07-01· 武乐圣典《太极武当》:武
- 07-01· 铁血英雄-现阶段战功牌兑
- 07-01· 2011年06月10日【原创】南歌
- 07-01· 【淘宝网信息】- 2010年的
- 07-01· 深圳品牌女装有哪些?
广告赞助商
最多关注
- 今日
- 本周
- 年度
最多评论
品牌排行榜