高斯正反算及换带计算matlab源代码-附截图

高斯投影坐标正、反算及相邻带的坐标换算MATLAB源代码主程序截图换带计算截图正算截图反算截图以上是课本算例截图精度0.001m和0.001”主程序程序引用代码MAIN程序名称主程序高斯投影坐标正、反算及相邻带的坐标换算MATLAB源代码子程序1disp('欢迎使用高斯投影正反算及相邻带的坐标换算程序');disp('1：高斯正算2：高斯反算3：换带计算');K=0;while(K1||K3)K=input('请根据上列选择计算类型K=');switchKcase1GSZS;case2GSFS;case3HDJS;otherwisedisp('K值无效（1-3）');enddisp('程序作者：桂林理工大学刘敬涛');disp('指导老师：王新桥刘斌卢献健');disp('引用请注明出处');end程序引用代码GSZS程序名称高斯正算functionGSZS%GSZS是将大地坐标换算为高斯坐标的子函数%此函数要调用DHH和HHD两个子函数%此函数包含子午线收敛角的计算disp('你选择的是高斯正算');B=input('输入大地坐标B=');L=input('输入大地坐标L=');L0=input('输入所用中央子午线L0=');B=DHH(B);L=DHH(L);L0=DHH(L0);disp('1:克拉索夫斯基椭球2:1975年国际椭球3:WGS-84椭球');T=0;while(T1||T2)T=input('请根据上列选择椭球模型T=');switchTcase1a=6378245.0000000000;b=6356863.0187730473;X=111134.861*(B*180/pi)-16036.480*sin(2*B)+16.828*sin(4*B)-0.022*sin(6*B);case2高斯投影坐标正、反算及相邻带的坐标换算MATLAB源代码子程序2a=6378140.0000000000;b=6356755.2881575287;X=111133.005*(B*180/pi)-16038.528*sin(2*B)+16.833*sin(4*B)-0.022*sin(6*B);otherwisedisp('T值无效（1-2）');endende=(sqrt(a^2-b^2))/a;e1=(sqrt(a^2-b^2))/b;V=sqrt(1+(e1^2)*(cos(B))^2);c=(a^2)/b;M=c/(V^3);N=c/V;t=tan(B);n=sqrt((e1^2)*(cos(B))^2);l=L-L0;xp1=X;xp2=(N*sin(B)*cos(B)*l^2)/2;xp3=(N*sin(B)*((cos(B))^3)*(5-t^2+9*n^2+4*n^4)*l^4)/24;xp4=(N*sin(B)*((cos(B))^5)*(61-58*t^2+t^4)*l^6)/720;x=xp1+xp2+xp3+xp4;yp1=N*cos(B)*l;yp2=N*(cos(B))^3*(1-t^2+n^2)*l^3/6;yp3=N*(cos(B))^5*(5-18*t^2+t^4+14*n^2-58*(n^2)*(t^2))*l^5/120;y=yp1+yp2+yp3;r1=l*sin(B);r2=(1/3)*sin(B)*(cos(B))^2*(l^3)*(1+3*n^2+2*n^4);r3=(1/15)*sin(B)*(cos(B))^2*(l^5)*(2-t^2);r=r1+r2+r3;formatlonggxyR=HHD(r)程序引用代码GSFS程序名称高斯反算functionGSFS%GSZS是将高斯坐标换算为大地坐标的子函数%此函数要调用DHH和HHD两个子函数%此函数包含子午线收敛角的计算disp('你选择的是高斯反算');x=input('输入高斯坐标x=');y=input('输入高斯坐标y=');高斯投影坐标正、反算及相邻带的坐标换算MATLAB源代码L0=input('输入所用中央子午线L0=');disp('1:克拉索夫斯基椭球2:1975年国际椭球3:WGS-84椭球');T=0;while(T1||T2)T=input('请根据上列选择椭球模型T=');switchTcase1a=6378245.0000000000;b=6356863.0187730473;B=x/6367558.4969;B=B+(50221746+(293622+(2350+22*(cos(B))^2)*(cos(B))^2)*(cos(B))^2)*10^(-10)*sin(B)*cos(B);case2a=6378140.0000000000;b=6356755.2881575287;B=x/6367452.1328;B=B+(50228976+(293697+(2383+22*(cos(B))^2)*(cos(B))^2)*(cos(B))^2)*10^(-10)*sin(B)*cos(B);otherwisedisp('T值无效（1-2）');endende=(sqrt(a^2-b^2))/a;e1=(sqrt(a^2-b^2))/b;V=sqrt(1+(e1^2)*(cos(B))^2);c=(a^2)/b;M=c/(V^3);N=c/V;t=tan(B);n=sqrt((e1^2)*(cos(B))^2);lp1=y/(N*cos(B));lp2=(1+2*t^2+n^2)*(y^3)/(6*cos(B)*N^3);lp3=(5+28*t^2+24*t^4+6*n^2+8*(n^2)*(t^2))*(y^5)/(120*cos(B)*N^5);l=lp1-lp2+lp3;Bp1=B;Bp2=(t*y^2)/(2*M*N);Bp3=(t/(24*M*N^3))*(5+3*t^2+n^2-9*(n^2)*(t^2))*y^4;Bp4=(t/(720*M*N^5))*(61+90*t^2+45*(t^4))*y^6;B=Bp1-Bp2+Bp3-Bp4;r1=l*sin(B);r2=(1/3)*sin(B)*(cos(B))^2*(l^3)*(1+3*n^2+2*n^4);r3=(1/15)*sin(B)*(cos(B))^2*(l^5)*(2-t^2);r=r1+r2+r3;formatlonggL=HHD(l)+L0B=HHD(B)高斯投影坐标正、反算及相邻带的坐标换算MATLAB源代码子程序3R=HHD(r)程序引用代码HDJS程序名称换带计算functionHDJSdisp('你选择的是换带计算')x=input('输入高斯坐标x=');y=input('输入高斯坐标y=');L0=input('输入原中央子午线L0=');disp('1:克拉索夫斯基椭球2:1975年国际椭球3:WGS-84椭球');T=0;while(T1||T2)T=input('请根据上列选择原椭球模型T=');switchTcase1a=6378245.0000000000;b=6356863.0187730473;B=x/6367558.4969;B=B+(50221746+(293622+(2350+22*(cos(B))^2)*(cos(B))^2)*(cos(B))^2)*10^(-10)*sin(B)*cos(B);case2a=6378140.0000000000;b=6356755.2881575287;B=x/6367452.1328;B=B+(50228976+(293697+(2383+22*(cos(B))^2)*(cos(B))^2)*(cos(B))^2)*10^(-10)*sin(B)*cos(B);otherwisedisp('T值无效（1-2）');endende=(sqrt(a^2-b^2))/a;e1=(sqrt(a^2-b^2))/b;V=sqrt(1+(e1^2)*(cos(B))^2);c=(a^2)/b;M=c/(V^3);N=c/V;t=tan(B);n=sqrt((e1^2)*(cos(B))^2);lp1=y/(N*cos(B));lp2=(1+2*t^2+n^2)*(y^3)/(6*cos(B)*N^3);lp3=(5+28*t^2+24*t^4+6*n^2+8*(n^2)*(t^2))*(y^5)/(120*cos(B)*N^5);l=lp1-lp2+lp3;Bp1=B;Bp2=(t*y^2)/(2*M*N);Bp3=(t/(24*M*N^3))*(5+3*t^2+n^2-9*(n^2)*(t^2))*y^4;高斯投影坐标正、反算及相邻带的坐标换算MATLAB源代码Bp4=(t/(720*M*N^5))*(61+90*t^2+45*(t^4))*y^6;B=Bp1-Bp2+Bp3-Bp4;formatlonggL=HHD(l)+L0B=HHD(B)L00=input('输入新中央子午线L00=');B=DHH(B);L=DHH(L);L00=DHH(L00);disp('1:克拉索夫斯基椭球2:1975年国际椭球3:WGS-84椭球');T=0;while(T1||T2)T=input('请根据上列选择新椭球模型T=');switchTcase1a=6378245.0000000000;b=6356863.0187730473;X=111134.861*(B*180/pi)-16036.480*sin(2*B)+16.828*sin(4*B)-0.022*sin(6*B);case2a=6378140.0000000000;b=6356755.2881575287;X=111133.005*(B*180/pi)-16038.528*sin(2*B)+16.833*sin(4*B)-0.022*sin(6*B);otherwisedisp('T值无效（1-2）');endende=(sqrt(a^2-b^2))/a;e1=(sqrt(a^2-b^2))/b;V=sqrt(1+(e1^2)*(cos(B))^2);c=(a^2)/b;M=c/(V^3);N=c/V;t=tan(B);n=sqrt((e1^2)*(cos(B))^2);l=L-L00;xp1=X;xp2=(N*sin(B)*cos(B)*l^2)/2;xp3=(N*sin(B)*((cos(B))^3)*(5-t^2+9*n^2+4*n^4)*l^4)/24;xp4=(N*sin(B)*((cos(B))^5)*(61-58*t^2+t^4)*l^6)/720;x=xp1+xp2+xp3+xp4;yp1=N*cos(B)*l;yp2=N*(cos(B))^3*(1-t^2+n^2)*l^3/6;yp3=N*(cos(B))^5*(5-18*t^2+t^4+14*n^2-58*(n^2)*(t^2))*l^5/120;高斯投影坐标正、反算及相邻带的坐标换算MATLAB源代码子程序4子程序5y=yp1+yp2+yp3;r1=l*sin(B);r2=(1/3)*sin(B)*(cos(B))^2*(l^3)*(1+3*n^2+2*n^4);r3=(1/15)*sin(B)*(cos(B))^2*(l^5)*(2-t^2);r=r1+r2+r3;formatlonggx2=xy2=yR=HHD(r)程序引用代码DHH程序名称度化弧度functionHD=DHH(BD)%DHH是将度化算为弧度的函数%度的表达形式形如（三十度四十五分二十三秒：30.4523）D=fix(BD);F=fix((BD-D)*100);M=((BD-D)*100-F)*100;HD=((D+F/60+M/3600)/180)*pi;程序引用代码HHD程序名称弧度化度functionBD=HHD(HD)%HHD是将弧度化算为度的函数%度的表达形式形如（三十度四十五分二十三秒：30.4523）SD=HD