matlab高斯正反算程序6,高斯正反算MATLAB程序

function [x,y]=zgauss(B,L,B0)

%B是弧度为单位

%L B0都是一度为单位

e2=0.006693421622966; %第一偏心率

ee2=0.006738525414683; %第二偏心率

a=6378245; %长半轴

X=111134.8611*B0-16036.4803*sin(2*B)+16.8281*sin(4*B)-0.022*sin(6*B); %子午弧长

p=3600*180/pi; %ρ"

W=sqrt(1-e2*(sin(B)^2)); %W

N=a/W; %卯酉圈半径

t=tan(B);

yita=ee2*cos(B)*cos(B); %η的平方

l=L-117; %L

ll=l*3600;

x=X+(N*sin(B)*cos(B)*ll^2)/(2*p^2)+N*sin(B)*(cos(B)^3)*(5-t^2-9*yita) *(ll^4)/(24*p^4);%x坐标

y=N*cos(B)*ll/p+N*(cos(B)^3)*(1-t^2+yita)*ll^3/(6*p^3)+N*(cos(B)^5)*( 5-18*t^2+t^4)*ll^5/(120*p^5);

function [BB,l]=fgauss(x,y)

%在克拉索夫椭球上

e2=0.006693421622966; %第一偏心率

ee2=0.006738525414683; %第二偏心率

a=6378245; %长半轴

a0=111134.8611;

Bf(1)=x/a0;

Bf(1)=Bf(1)*pi/180; %把度换算成弧度

%以下用迭代法解出 Bf

for i=1:5

Bf(i+1)=(x-(-16036.4803*sin(2*Bf(i))+16.8281*sin(4*Bf(i))-0.022*sin(6 *Bf(i))))/a0*pi/180;

end

%代入公式计算

BBf=Bf(5);

tf=tan(BBf);

yitaf2=ee2*(cos(BBf)^2);

Wf=sqrt(1-e2*(sin(BBf)^2));

Mf=a*(1-e2)/Wf;

Nf=a/Wf;

BB=BBf-tf*y^2/(2*Mf*Nf)+tf^3*(5+3*tf^2+yitaf2-9*yitaf2*tf^2)*y^4/(24* Mf*Nf^3);