《matlab高斯正反算及换带》由会员分享，可在线阅读，更多相关《matlab高斯正反算及换带（9页珍藏版）》请在金锄头文库上搜索。
1、本程序采用 matlab 编程， 用于测绘工程同学的学习。 参考网上资料， 加以学习。 内容包括： 高斯正算 高斯反算 换带计算 Bf 的计算 X 的计算 涉及度化弧度 涉及弧度化度 By cowboy function main disp(欢迎使用高斯投影正反算及相邻带的坐标换算程序); disp(1：高斯正算 2：高斯反算 3：换带计算); K=0; while (K3) K=input(请根据上列选择计算类型 K=); switch K case 1 GSZS; case 2 GSFS; case 3 HDJS; otherwise disp(K 值无效 （1-3）); end disp
2、(程序作者：屌丝); disp(指导老师：高富帅); end function GSZS %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 (T3) T=input(请根据上列选择椭球模
3、型 T=); switch T case 1 a=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); case 2 a=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); case 3 a=6378137.0000000000; f=1
4、/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*c os(B)*(sin(B)5 otherwise disp(T 值无效 （1-3）); end end e=(sqrt(a2-b2)/a; e1=(sqrt(a2-b2)/b; V=sqrt(1+(e12)*(cos(B)2); c=(a2)/b; M=c/(V3); N=c/V; t=tan(B); n=sqrt(e12)*(cos(B)2); l=L-L0; xp1=X; xp2=(N*
5、sin(B)*cos(B)*l2)/2; xp3=(N*sin(B)*(cos(B)3)*(5-t2+9*n2+4*n4)*l4)/24; xp4=(N*sin(B)*(cos(B)5)*(61-58*t2+t4)*l6)/720; x=xp1+xp2+xp3+xp4; yp1=N*cos(B)*l; yp2=N*(cos(B)3*(1-t2+n2)*l3/6; yp3=N*(cos(B)5*(5-18*t2+t4+14*n2-58*(n2)*(t2)*l5/120; y=yp1+yp2+yp3; r1=l*sin(B); r2=(1/3)*sin(B)*(cos(B)2*(l3)*(1+3*
6、n2+2*n4); r3=(1/15)*sin(B)*(cos(B)2*(l5)*(2-t2); r=r1+r2+r3; format long g x y R=HHD(r) End function GSFS %GSZS 是将高斯坐标换算为大地坐标的子函数 %此函数要调用 DHH 和 HHD 两个子函数 %此函数包含子午线收敛角的计算 disp(你选择的是高斯反算); x=input(输入高斯坐标 x=); y=input(输入高斯坐标 y=); L0=input(输入所用中央子午线 L0=); disp(1:克拉索夫斯基椭球 2:1975 年国际椭球 3:WGS-84 椭球); T=0;
7、while (T3) T=input(请根据上列选择椭球模型 T=); switch T case 1 a=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)*si n(B)*cos(B); case 2 a=6378140.0000000000; b=6356755.2881575287; B=x/6367452.1328; B=B+(50228976+(293697+(2383+22*(
8、cos(B)2)*(cos(B)2)*(cos(B)2)*10(-10)*si n(B)*cos(B); case 3 a=6378137.0000000000; f=1/298.257223563; b=a*(1-f); e=(sqrt(a2-b2)/a; m0=a*(1-e2); m2=(3/2)*e2*m0; m4=(5/4)*e2*m2; m6=(7/6)*e2*m4; m8=(9/8)*e2*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=
9、(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; %make B0 != B1 so as to pass the while test at first time B1 = 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 %main loop B
10、=B1; %Bf otherwise disp(T 值无效 （1-3）); end end e=(sqrt(a2-b2)/a; e1=(sqrt(a2-b2)/b; V=sqrt(1+(e12)*(cos(B)2); c=(a2)/b; M=c/(V3); N=c/V; t=tan(B); n=sqrt(e12)*(cos(B)2); lp1=y/(N*cos(B); lp2=(1+2*t2+n2)*(y3)/(6*cos(B)*N3); lp3=(5+28*t2+24*t4+6*n2+8*(n2)*(t2)*(y5)/(120*cos(B)*N5); l=lp1-lp2+lp3; Bp1=B
11、; Bp2=(t*y2)/(2*M*N); Bp3=(t/(24*M*N3)*(5+3*t2+n2-9*(n2)*(t2)*y4; Bp4=(t/(720*M*N5)*(61+90*t2+45*(t4)*y6; B=Bp1-Bp2+Bp3-Bp4; r1=l*sin(B); r2=(1/3)*sin(B)*(cos(B)2*(l3)*(1+3*n2+2*n4); r3=(1/15)*sin(B)*(cos(B)2*(l5)*(2-t2); r=r1+r2+r3; format long g L=HHD(l)+L0 B=HHD(B) R=HHD(r) function HDJS disp(你选择
12、的是换带计算) x=input(输入高斯坐标 x=); y=input(输入高斯坐标 y=); L0=input(输入原中央子午线 L0=); disp(1:克拉索夫斯基椭球 2:1975 年国际椭球 3:WGS-84 椭球); T=0; while (T3) T=input(请根据上列选择原椭球模型 T=); switch T case 1 a=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(-
13、10)*si n(B)*cos(B); case 2 a=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)*si n(B)*cos(B); case 3 a=6378137.0000000000; f=1/298.257223563; b=a*(1-f); e=(sqrt(a2-b2)/a; m0=a*(1-e2); m2=(3/2)*e2*m0; m4=(5/4)*e2*m2; m6=(
14、7/6)*e2*m4; m8=(9/8)*e2*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; %make B0 != B1 so as to pass the while test at first time B1 = 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 %main loop B=B1; %Bf otherwise disp(T 值无效 （
