matlab-高斯正反算以及换带

本程序采用matlab编程，用于测绘工程同学的学习。参考网上资料，加以学习。内容包括：高斯正算高斯反算换带计算Bf的计算X的计算涉及度化弧度涉及弧度化度Bycowboyfunctionmaindisp('欢迎使用高斯投影正反算及相邻带的坐标换算程序');disp('1：高斯正算2：高斯反算3：换带计算');K=0;while(K1||K3)K=input('请根据上列选择计算类型K=');switchKcase1GSZS;case2GSFS;case3HDJS;otherwisedisp('K值无效（1-3）');enddisp('程序作者：屌丝');disp('指导老师：高富帅');endfunctionGSZS%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||T3)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);case3a=6378137.0000000000;f=1/298.257223563;b=a*(1-f);X=6367449.1458*B-32009.8185*cos(B)*sin(B)-133.9975*cos(B)*(sin(B))^3-0.6975*cos(B)*(sin(B))^5otherwisedisp('T值无效（1-3）');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)EndfunctionGSFS%GSZS是将高斯坐标换算为大地坐标的子函数%此函数要调用DHH和HHD两个子函数%此函数包含子午线收敛角的计算disp('你选择的是高斯反算');x=input('输入高斯坐标x=');y=input('输入高斯坐标y=');L0=input('输入所用中央子午线L0=');disp('1:克拉索夫斯基椭球2:1975年国际椭球3:WGS-84椭球');T=0;while(T1||T3)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);case3a=6378137.0000000000;f=1/298.257223563;b=a*(1-f);e=(sqrt(a^2-b^2))/a;m0=a*(1-e^2);m2=(3/2)*e^2*m0;m4=(5/4)*e^2*m2;m6=(7/6)*e^2*m4;m8=(9/8)*e^2*m6;a0=m0+(1/2)*m2+(3/8)*m4+(5/16)*m6+(35/128)*m8;a2=(1/2)*m2+(1/2)*m4+(15/32)*m6+(7/16)*m8;a4=(1/8)*m4+(3/16)*m6+(7/32)*m8;a6=(1/32)*m6+(1/16)*m8;a8=(1/128)*m8;B0=x/a0+10*eps;%makeB0!=B1soastopassthewhiletestatfirsttimeB1=x/a0;while(abs(B1-B0)eps)B0=B1;B1=(x-(-a2*sin(B0*2)*0.5+a4*sin(B0*4)*0.25-a6*sin(B0*6)/6.0+a8*sin(B0*8)*0.125))/a0;end%mainloopB=B1;%Bfotherwisedisp('T值无效（1-3）');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)R=HHD(r)functionHDJSdisp('你选择的是换带计算')x=input('输入高斯坐标x=');y=input('输入高斯坐标y=');L0=input('输入原中央子午线L0=');disp('1:克拉索夫斯基椭球2:1975年国际椭球3:WGS-84椭球');T=0;while(T1||T3)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);case3a=6378137.0000000000;f=1/298.257223563;b=a*(1-f);e=(sqrt(a^2-b^2))/a;m0=a*(1-e^2);m2=(3/2)*e^2*m0;m4=(5/4)*e^2*m2;m6=(7/6)*e^2*m4;m8=(9/8)*e^2*m6;a0=m0+(1/2)*m2+(3/8)*m4+(5/16)*m6+(35/128)*m8;a2=(1/2)*m2+(1/2)*m4+(15/32)*m6+(7/16)*m8;a4=(1/8)*m4+(3/16)*m6+(7/32)*m8;a6=(1/32)*m6+(1/16)*m8;a8=(1/128)*m8;B0=x/a0+10*eps;%makeB0!=B1soastopassthewhiletestatfirsttimeB1=x/a0;while(abs(B1-B0)eps)B0=B1;B1=(x-(-a2*sin(B0*2)*0.5+a4*sin(B0*4)*0.25-a6*sin(B0*6)/6.0+a8*sin(B0*8)*0.125))/a0;end%mainloopB=B1;%Bfotherwisedisp('T值无效（1-3）');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;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||T3)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*