当前位置:文档之家› 山东农业大学大地测量学课程设计C++

山东农业大学大地测量学课程设计C++

一、通过四参数进行平面坐标转换#include <iostream>#include <math.h>using namespace std;int main(){ double x1,y1,x2,y2,x3,y3,x4,y4,z,k,x0,y0,x,y,X,Y,Z,c,a,b;cout<<"请输入原坐标系两点坐标:"<<endl;cout<<"A ";cin>>x1>>y1;cout<<"B "; cin>>x2>>y2;cout<<"请输入转换后两点坐标:"<<endl;cout<<"A' "; cin>>x3>>y3;cout<<"B' ";cin>>x4>>y4;z=atan((y2-y1)/(x2-x1))-atan((y4-y3)/(x4-x3));if(z/3.1415926*180>90)z=3.1415926-z;k=(sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)))/(sqrt((x4-x3)*(x4-x3)+(y4-y3)*(y4-y3)));x0=x3-k*cos(z)*x1+k*sin(z)*y1;y0=y3-k*sin(z)*x1-k*cos(z)*y1;a=int(z/3.1415926*180);b=int((z/3.1415926*180-a)*60);c=int(((z/3.1415926*180-a)*60-b)*60);cout<<"转换角:"<<a<<"°"<<b<<"′"<<c<<"″"<<"缩放比例:"<<k<<" 平移参数X:"<<x0<<" 平移参数Y:"<<y0<<endl;cout<<"请输入待求点原坐标:";cin>>x>>y;X=x0+k*cos(z)*x-k*sin(z)*y;Y=y0+k*cos(z)*y+k*sin(z)*x;cout<<"转换后坐标:"<<"X= "<<X<<" 、"<<"Y= "<<Y<<endl;return 0;}图1、四参数坐标变换此例子是老师上课讲的例题,程序中没有控制小数位数,都是两位。

二、大地坐标与空间直角坐标的转换#include<iostream>#include<cmath>#include<iomanip>#include <math.h>#define E 0.006693422#define A 6378245#define PAI 3.141592654using namespace std;int main(){ void kd();void dk(); int m;cout<<"空间直角坐标换算到大地坐标输入 1 ;大地坐标换算到空间直角坐标输入2 "<<endl;cin>>m;if(m==1) kd();else if(m==2) dk();else cout<<"输入错误,重新输入"<<endl;return 0;}void kd(){ double x,y,z,b1,b2,b,l,h,d,f,m,D,F,M;cout<<"输入空间直角坐标X,Y,Z"<<endl;cin>>x; cin>>y; cin>>z;b1=atan(z/sqrt(x*x+y*y));b2=atan(((A*E*sin(b1))/sqrt(1-E*pow(sin(b1),2))+z)/(sqrt(x*x+y*y)));while(fabs(b2-b1)>0.000001){ b1=b2; b2=atan(((A*E*sin(b1))/sqrt(1-E*pow(sin(b1),2))+z)/(sqrt(x*x+y*y))); }b=b2; l=atan(y/x); h=sqrt(x*x+y*y)/cos(b)-A/sqrt(1-E*pow(sin(b),2));b=b*180/PAI; l=l*180/PAI;d=int(b);f=int((b-d)*60);m=((b-d)*60-f)*60;D=int(l);F=int((l-D)*60);M=((l-D)*60-F)*60;cout<<"大地纬度B="<<d<<"°"<<f<<"′"<<fixed<<setprecision(4)<<m<<"″"<<endl;cout<<"大地经度L="<<D<<"°"<<F<<"′"<<fixed<<setprecision(4)<<M<<"″"<<endl;cout<<"大地高H="<<fixed<<setprecision(3)<<h<<endl; }void dk(){ double B,L,H,x,y,z,d,f,m,D,F,M;cout<<"输入大地坐标B,L,H"<<endl;cout<<" °′″"<<endl;cin>>d>>f>>m;cin>>D>>F>>M;cin>>H;B=(d+f/60+m/3600)*PAI/180;L=(D+F/60+M/3600)*PAI/180;cout<<"弧度制B="<<B<<endl;cout<<"弧度制L="<<L<<endl;x=(A/sqrt(1-E*pow(sin(B),2))+H)*cos(B)*cos(L);y=(A/sqrt(1-E*pow(sin(B),2))+H)*cos(B)*sin(L);z=(A*(1-E)/sqrt(1-E*pow(sin(B),2))+H)*sin(B);cout<<"空间坐标系X="<<fixed<<setprecision(3)<<x<<endl;cout<<"空间坐标系Y="<<y<<endl;cout<<"空间坐标系Z="<<z<<endl;}图2、XYZ-BHL图3、BHL-XYZ此程序是依照《大地测量学基础》补充作业和辅助材料中作业四来完成,实例是作业四中的实例,此程序中没有提供椭球参数的选择,固定的给出了a=6378245,e方=0.006693422.程序不够完善。

三、高斯投影正反算#include<iostream>#include<cmath>#include<iomanip>#include <math.h>#include <stdio.h>#include <stdlib.h>#include <malloc.h>#define PAI 3.141592654using namespace std;int main(){ void zs();void fs(); int m,h;cout<<"高斯投影正算,请按1 ;高斯投影反算,请按2 "<<endl;cin>>m;if(m==1)zs();else if(m==2)fs();else cout<<"输入错误,重新输入"<<endl;return 0;}void zs(){cout<<"请选择坐标系:"<<endl;cout<<"选择北京-54坐标系,请按1"<<endl;cout<<"选择西安-80坐标系,请按2"<<endl;int h;double d,f,m,D,F,M,N,a0,a4,a6,a3,a5,x,y,l,bm,rm,B,L,L0;cin>>h;if(h==1){cout<<"B= °′″"<<endl;cin>>d>>f>>m;cout<<"L= °′″"<<endl;cin>>D>>F>>M;B=(d+f/60+m/3600)*PAI/180;L=(D+F/60+M/3600);bm=d*3600+f*60+m;rm=206264.8063;if(((int(L/6)+1)*6-3)>L)L0=((int(L/6))*6-3);else L0=((int(L/6)+1)*6-3);l=(L-L0)*3600/rm;N=6399698.902-(21562.267-(108.973-0.612*pow(cos(B),2))*pow(cos(B),2))*pow(cos(B),2);a0=32140.404-(135.3302-(0.7092-0.0040*pow(cos(B),2))*pow(cos(B),2))*pow(cos(B),2);a4=(0.25+0.00252*pow(cos(B),2))*pow(cos(B),2)-0.04166;a6=(0.166*pow(cos(B),2)-0.084)*pow(cos(B),2);a3=(0.3333333+0.001123*pow(cos(B),2))*pow(cos(B),2)-0.1666667;a5=0.0083-(0.1667-(0.1968+0.0040*pow(cos(B),2))*pow(cos(B),2))*pow(cos(B),2);x=6367558.4969*bm/rm-(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);cout<<"采用北京-54坐标系、"<<endl;cout<<"高斯投影正算后x=:"<<fixed<<setprecision(4)<<x<<"高斯投影正算后y=:"<<fixed<<setprecision(4)<<y<<endl;}if(h==2){cout<<"B= °′″"<<endl;cin>>d>>f>>m;cout<<"L= °′″"<<endl;cin>>D>>F>>M;B=(d+f/60+m/3600)*PAI/180;L=(D+F/60+M/3600);bm=d*3600+f*60+m;rm=206264.8063;L0=((int(L/6)+1)*6-3);if(((int(L/6)+1)*6-3)>L)L0=((int(L/6))*6-3);else L0=((int(L/6)+1)*6-3);l=(L-L0)*3600/rm;N=6399596.652-(21565.045-(108.9996-0.603*pow(cos(B),2))*pow(cos(B),2))*pow(cos(B),2);a0=32144.5189-(135.3646-(0.7034-0.0041*pow(cos(B),2))*pow(cos(B),2))*pow(cos(B),2);a4=(0.25+0.00253*pow(cos(B),2))*pow(cos(B),2)-0.04167;a6=(0.167*pow(cos(B),2)-0.083)*pow(cos(B),2);a3=(0.3333333+0.001123*pow(cos(B),2))*pow(cos(B),2)-0.1666667;a5=0.00878-(0.1702-0.20382*pow(cos(B),2))*pow(cos(B),2);x=6367452.1328*bm/rm-(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);cout<<"采用西安-80坐标系、"<<endl;cout<<"高斯投影正算后x=:"<<fixed<<setprecision(4)<<x<<"高斯投影正算后y=:"<<fixed<<setprecision(4)<<y<<endl;}}void fs(){cout<<"请选择坐标系:"<<endl;cout<<"选择北京-54坐标系,请按1"<<endl;cout<<"选择西安-80坐标系,请按2"<<endl;int h; double x,y,rm,b,bf,B,l,nf,b2,b3,b4,b5,z,d,f,m,L,L0,D,M;cin>>h;if(h==1){cout<<"请输入X:"<<endl;cin>>x;cout<<"请输入Y: "<<endl;cin>>y;cout<<"请输入LO(中央经度): "<<endl;cin>>L0;rm=206264.8063;b =x/6367558.4969;bf=b+(50221746+(293622+(2350+22*pow(cos(b),2))*pow(cos(b),2))*pow(cos(b),2))*0.00000000 01*sin(b)*cos(b);nf=6399698.902-(21562.267-(108.973-0.612*pow(cos(bf),2))*pow(cos(bf),2))*pow(cos(bf),2);z=y/(nf*cos(bf));b2=(0.5+0.003369*pow(cos(bf),2))*sin(bf)*cos(bf);b3=0.333333-(0.166667-0.001123*pow(cos(bf),2))*pow(cos(bf),2);b4=0.25+(0.16161+0.00562*pow(cos(bf),2))*pow(cos(bf),2);b5=0.2-(0.1667-0.0088*pow(cos(bf),2))*pow(cos(bf),2);B=bf*rm-(1-(b4-0.12*z*z)*z*z)*z*z*b2*rm;l=(1-(b3-b5*z*z)*z*z)*z*rm;cout<<"采用北京-54坐标系、"<<endl;cout<<"B=(弧度制)"<<fixed<<setprecision(4)<<B<<" l=:"<<fixed<<setprecision(4)<<l<<endl;d=B/3600;f=int((d-int(B/3600))*60);m=(d-int(d)-f/60)*3600;cout<<"B="<<int(d)<<"°"<<fixed<<setprecision(0)<<f<<"′"<<fixed<<setprecision(4)<<m<<"″"<<endl;L=L0+int(l/3600);D=int((l/3600-int(l/3600))*60);M=(((l/3600-int(l/3600))*60)-D)*60;cout<<"L="<<fixed<<setprecision(0)<<L<<"°"<<fixed<<setprecision(0)<<D<<"′"<<fixed<<setprecision(4)<<M<<"″"<<endl;}if(h==2){cout<<"请输入X:"<<endl;cin>>x;cout<<"请输入Y: "<<endl;cin>>y;cout<<"请输入LO(中央经度): "<<endl;cin>>L0;rm=206264.8063;b =x/6367452.133;bf=b+(50228976+(293697+(2383+22*pow(cos(b),2))*pow(cos(b),2))*pow(cos(b),2))*0.00000000 01*sin(b)*cos(b);nf=6399698.902-(21562.267-(108.973-0.612*pow(cos(bf),2))*pow(cos(bf),2))*pow(cos(bf),2);z=y/(nf*cos(bf));b2=(0.5+0.00336975*pow(cos(bf),2))*sin(bf)*cos(bf);b3=0.333333-(0.166667-0.001123*pow(cos(bf),2))*pow(cos(bf),2);b4=0.25+(0.161612+0.005617*pow(cos(bf),2))*pow(cos(bf),2);b5=0.2-(0.1667-0.00878*pow(cos(bf),2))*pow(cos(bf),2);B=bf*rm-(1-(b4-0.12*z*z)*z*z)*z*z*b2*rm;l=(1-(b3-b5*z*z)*z*z)*z*rm;cout<<"采用西安-80坐标系、"<<endl;cout<<"B=(弧度制)"<<fixed<<setprecision(4)<<B<<" l=:"<<fixed<<setprecision(4)<<l<<endl;d=B/3600;f=int((d-int(B/3600))*60);m=(d-int(d)-f/60)*3600;cout<<"B="<<int(d)<<"°"<<fixed<<setprecision(0)<<f<<"′"<<fixed<<setprecision(4)<<m<<"″"<<endl;L=L0+int(l/3600);D=int((l/3600-int(l/3600))*60);M=(((l/3600-int(l/3600))*60)-D)*60;cout<<"L="<<fixed<<setprecision(0)<<L<<"°"<<fixed<<setprecision(0)<<D<<"′"<<fixed<<setprecision(4)<<M<<"″"<<endl;}}高斯投影正反算实例:图4、高斯投影正算-北京54图5、高斯投影正算-西安80图6、高斯投影正算-北京80图7、高斯投影正算-西安80图8、高斯投影反算-北京54图9、高斯投影反算-西安80此部分西安-80坐标系下的高斯投影反算结果不对,是因为在此部分Nf的计算是用了北京54坐标系下Nf的计算公式,因此出现小数部分的错误。

相关主题