高斯正反算的程序。vc编辑

//大地测量学基础作业 //任意椭球坐标系的高斯投影正反算 //2003年11月 #include <math.h> #include <stdio.h> #include <stdlib.h> #include <malloc.h> #define PI 3.141592653589793 double DMS2RAD(double dmsAngle) { int degAngle,minAngle; double radAngle,secAngle; 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; return radAngle; } double RAD2DMS(double radAngle) { int degAngle,minAngle; double secAngle,dmsAngle; 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; 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; /*计算a0 a2 a4 a6 a8*/ 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; /*计算a0 a2 a4 a6 a8*/ 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*/ X=a0*radlat-sb*cb*((a2-a4+a6)+(2*a4-16*a6/3)*sb*sb+16*a6*sb*sb*sb*sb/3.0); /*计算卯酉圈半径N*/ 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; } /*输入l0，已知点坐标*/ printf("请输入l0:\n"); scanf("%lf",&l0); l0=l0*3.1415926535897932384626433832795/180.0; printf("请输入x坐标:\n"); scanf("%lf",&coor_x); printf("请输入y坐标:\n"); scanf("%lf",&coor_y); /*计算b,e1,e*/ b=a*(1-Alfa); e1=sqrt((a/b)*(a/b)-1); e=sqrt(1-(b/a)*(b/a)); /*计算a0 a2 a4 a6 a8*/ 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,v,N,M*/ 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,dietal*/ 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)为：%15.9lf\n",dmslat); printf("已知点的经度(ddd.mmss)为：%15.9lf\n",dmslon); } }