高斯正反算及换带计算matlab源代码_附截图
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
高斯投影坐标正、反算及相邻带的坐标换算 MATLA源B 代码
主程序截图
正算截图
换带计算截图
反算截图
以上是课本算例截图 精度 0.001m 和 0.001” 主程序
程序引用代码
MAIN
程序名称
主程序
高斯投影坐标正、反算及相邻带的坐标换算
disp(' 欢迎使用高斯投影正反算及相邻带的坐标换算程序 ');
disp('T 值无效 （1-2）'); end
end
e=(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);
otherwise
disp('T 值无效 （1-2）'); end end
e=(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);
r3=(1/15)*sin(B)*(cos(B))^2*(l^5)*(2-t^2);
r=r1+r2+r3;
format long g
L=HHD(l)+L0
B=HHD(B)
高斯投影坐标正、反算及相邻带的坐标换算
R=HHD(r)
子程序 3
MATLA源B 代码
程序引用代码
HDJS
程序名称
换带计算
function HDJS
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;
disp(' 你选择的是换带计算 ')
x=input(' 输入高斯坐标 x=');
y=input(' 输入高斯坐标 y=');
L0=input(' 输入原中央子午线 L0=');
disp('1: 克拉索夫斯基椭球
2:1975 年国际椭球
T=0;
while (T<1||T>2)
T=input(' 请根据上列选择原椭球模型 T=');
case 2
高斯投影坐标正、反算及相邻带的坐标换算 MATLA源B 代码
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); otherwise disp('T 值无效 （1-2）'); end end e=(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; format long g x y R=HHD(r)
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;
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);
刘敬涛 ');
disp(' 指导老师：王新桥 刘斌 卢献健 ');
disp(' 引用请注明出处 '); end
子程序 1
MATLA源B 代码
程序引用代码
GSZS
程序名称
高斯正算
function GSZS %GSZS是将大地坐标换算为高斯坐标的子函数
%此函数要调用 DHH 和 HHD 两个子函数
%此函数包含子午线收敛角的计算
子程序 4
MATLA源B 代码
程序引用代码
DHH
function HD=DHH(BD) %DHH 是将度化算为弧度的函数 %度的表达形式形如（三十度四十五分二十三秒： D=fix(BD);
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);
b=6356755.2881575287;
X=111133.005*(B*180/pi)-16038.528*sin(2*B)+16.833*sin(4*B)-0.022*sin(6*B);
otherwise
disp('T 值无效 （1-2）'); end
end
e=(sqrt(a^2-b^2))/a;
3:WGS-84 椭球 ');
T=0;
while (T<1||T>2)
T=input(' 请根据上列选择椭球模型 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);
程序名称
高斯反算
高斯投影坐标正、反算及相邻带的坐标换算 MATLA源B 代码
L0=input(' 输入所用中央子午线 L0=');
disp('1: 克拉索夫斯基椭球
2:1975 年国际椭球
3:WGS-84 椭球 ');
T=0;
while (T<1||T>2)
T=input(' 请根据上列选择椭球模型 T=');
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; format long g x2=x y2=y R=HHD(r)
MATLA源B 代码
L00=input(' 输入新中央子午线 L00=');
B=DHH(B);
L=DHH(L);
L00=DHH(L00);
disp('1: 克拉索夫斯基椭球 T=0;
2:1975 年国际椭球
3:WGS-84 椭球 ');
while (T<1||T>2)
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)*sin(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)*sin(B)*cos(B);
otherwise
disp('1 ：高斯正算
2：高斯反算
3：换带计算 ');
K=0;
while (K<1||K>3)
K=input(' 请根据上列选择计算类型 K=');
switch K
case 1
GSZS;
case 2
GSFS;
case 3 HDJS;
otherwise
disp('K 值无效 （ 1-3）');
end disp(' 程序作者：桂林理工大学
子程序 2
程序引用代码
GSFS
function GSFS %GSZS是将高斯坐标换算为大地坐标的子函数 %此函数要调用 DHH 和 HHD 两个子函数 %此函数包含子午线收敛角的计算 disp(' 你选择的是高斯反算 '); x=input(' 输入高斯坐标 x='); y=input(' 输入高斯坐标 y=');
switch T
3:WGS-84 椭球 ');
case 1 a=6378245.0000000000;
Hale Waihona Puke Baidu
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);
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);
disp(' 你选择的是高斯正算 '); B=input(' 输入大地坐标 B=');
L=input(' 输入大地坐标 L=');
L0=input(' 输入所用中央子午线 L0=');
B=DHH(B);
L=DHH(L);
L0=DHH(L0);
disp('1: 克拉索夫斯基椭球
2:1975 年国际椭球
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;
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)*sin(B)*cos(B);
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; format long g L=HHD(l)+L0 B=HHD(B)
主程序截图
正算截图
换带计算截图
反算截图
以上是课本算例截图 精度 0.001m 和 0.001” 主程序
程序引用代码
MAIN
程序名称
主程序
高斯投影坐标正、反算及相邻带的坐标换算
disp(' 欢迎使用高斯投影正反算及相邻带的坐标换算程序 ');
disp('T 值无效 （1-2）'); end
end
e=(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);
otherwise
disp('T 值无效 （1-2）'); end end
e=(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);
r3=(1/15)*sin(B)*(cos(B))^2*(l^5)*(2-t^2);
r=r1+r2+r3;
format long g
L=HHD(l)+L0
B=HHD(B)
高斯投影坐标正、反算及相邻带的坐标换算
R=HHD(r)
子程序 3
MATLA源B 代码
程序引用代码
HDJS
程序名称
换带计算
function HDJS
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;
disp(' 你选择的是换带计算 ')
x=input(' 输入高斯坐标 x=');
y=input(' 输入高斯坐标 y=');
L0=input(' 输入原中央子午线 L0=');
disp('1: 克拉索夫斯基椭球
2:1975 年国际椭球
T=0;
while (T<1||T>2)
T=input(' 请根据上列选择原椭球模型 T=');
case 2
高斯投影坐标正、反算及相邻带的坐标换算 MATLA源B 代码
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); otherwise disp('T 值无效 （1-2）'); end end e=(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; format long g x y R=HHD(r)
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;
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);
刘敬涛 ');
disp(' 指导老师：王新桥 刘斌 卢献健 ');
disp(' 引用请注明出处 '); end
子程序 1
MATLA源B 代码
程序引用代码
GSZS
程序名称
高斯正算
function GSZS %GSZS是将大地坐标换算为高斯坐标的子函数
%此函数要调用 DHH 和 HHD 两个子函数
%此函数包含子午线收敛角的计算
子程序 4
MATLA源B 代码
程序引用代码
DHH
function HD=DHH(BD) %DHH 是将度化算为弧度的函数 %度的表达形式形如（三十度四十五分二十三秒： D=fix(BD);
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);
b=6356755.2881575287;
X=111133.005*(B*180/pi)-16038.528*sin(2*B)+16.833*sin(4*B)-0.022*sin(6*B);
otherwise
disp('T 值无效 （1-2）'); end
end
e=(sqrt(a^2-b^2))/a;
3:WGS-84 椭球 ');
T=0;
while (T<1||T>2)
T=input(' 请根据上列选择椭球模型 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);
程序名称
高斯反算
高斯投影坐标正、反算及相邻带的坐标换算 MATLA源B 代码
L0=input(' 输入所用中央子午线 L0=');
disp('1: 克拉索夫斯基椭球
2:1975 年国际椭球
3:WGS-84 椭球 ');
T=0;
while (T<1||T>2)
T=input(' 请根据上列选择椭球模型 T=');
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; format long g x2=x y2=y R=HHD(r)
MATLA源B 代码
L00=input(' 输入新中央子午线 L00=');
B=DHH(B);
L=DHH(L);
L00=DHH(L00);
disp('1: 克拉索夫斯基椭球 T=0;
2:1975 年国际椭球
3:WGS-84 椭球 ');
while (T<1||T>2)
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)*sin(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)*sin(B)*cos(B);
otherwise
disp('1 ：高斯正算
2：高斯反算
3：换带计算 ');
K=0;
while (K<1||K>3)
K=input(' 请根据上列选择计算类型 K=');
switch K
case 1
GSZS;
case 2
GSFS;
case 3 HDJS;
otherwise
disp('K 值无效 （ 1-3）');
end disp(' 程序作者：桂林理工大学
子程序 2
程序引用代码
GSFS
function GSFS %GSZS是将高斯坐标换算为大地坐标的子函数 %此函数要调用 DHH 和 HHD 两个子函数 %此函数包含子午线收敛角的计算 disp(' 你选择的是高斯反算 '); x=input(' 输入高斯坐标 x='); y=input(' 输入高斯坐标 y=');
switch T
3:WGS-84 椭球 ');
case 1 a=6378245.0000000000;
Hale Waihona Puke Baidu
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);
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);
disp(' 你选择的是高斯正算 '); B=input(' 输入大地坐标 B=');
L=input(' 输入大地坐标 L=');
L0=input(' 输入所用中央子午线 L0=');
B=DHH(B);
L=DHH(L);
L0=DHH(L0);
disp('1: 克拉索夫斯基椭球
2:1975 年国际椭球
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;
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)*sin(B)*cos(B);
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; format long g L=HHD(l)+L0 B=HHD(B)