format long;
x1 = input('请输入x1=');
y1= input('请输入y1=');
x2 = input('请输入x2=');
y2 = input('请输入y2=');
mx1 = input('请输入mx1=');
my1 = input('请输入my1=');
mx2 = input('请输入mx2=');
my2 = input('请输入my2=');
DX = x2 - x1;
DY = y2 - y1;
SAB = sqrt(DX*DX+DY*DY);
mDX = sqrt(mx1*mx1+mx2*mx2);
mDY = sqrt(my1*my1+my2*my2);
mSAB = sqrt((DX*mDX+DY*mDY)/sqrt(DX*DX+DY*DY));
m1= sqrt(abs((1/(1+(DY/DX)^2))*((mDY*DX-mDX*DY)/DX^2)))
m2 = sqrt(abs(((-1).*(1-DY*DY/(SAB*SAB)))*((DY*mDY*SAB*SAB-SAB*mSAB*DY*DY)/(SAB^4))))
if m1<m2
disp('方位角为：');
format long
if x1== x2 && y1 == y2
error('输入同一个点，请检查！');
elseif x1 == x2 && y1<y2
azimuth = 1/2*pi;
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
elseif x1 == x2 && y1>y2
azimuth=3/2*pi;
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
elseif y1 == y2 && x1<x2
azimuth=0;
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
elseif y1 == y2 && x1<x2
azimuth=pi;
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
elseif x1 > x2
azimuth=atan((y2-y1)/(x2-x1)) + pi;
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
elseif x1 < x2 && y1<y2
azimuth=atan((y2-y1)/(x2-x1));
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
elseif x1 < x2 && y1>y2
azimuth=atan((y2-y1)/(x2-x1)) + 2*pi;
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
end
D=sqrt(((x2-x1)^2)+((y2-y1)^2));
disp('两点间的距离为：');
fprintf('%04.4f米\n',SAB)
else
disp('方位角为：');
format long
if x1== x2 && y1 == y2
error('输入同一个点，请检查！');
elseif x1 == x2 && y1<y2
azimuth = 1/2*pi;
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
elseif x1 == x2 && y1>y2
azimuth=3/2*pi;
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
elseif y1 == y2 && x1<x2
azimuth=0;
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
elseif y1 == y2 && x1<x2
azimuth=pi;
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
elseif x1 > x2
azimuth=asin((y2-y1)/SAB) + pi;
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
elseif x1 < x2 && y1<y2
azimuth=asin((y2-y1)/SAB);
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
elseif x1 < x2 && y1>y2
azimuth=asin((y2-y1)/SAB) + 2*pi;
azimuth=azimuth/pi*180;
azimuth=degree2dms(azimuth);
end
disp('两点间的距离为：');
fprintf('%04.4f米\n',SAB)
end
X=[y1,y2];
Y=[x1,x2];
%绘图
plot(X,Y)
hold on
plot(X,Y,'x')
axis equal
title ('方位角计算')
xlabel('y轴')
ylabel('x轴')
text(y1,x1,['A','(',num2str(x1),',',num2str(y1),')'])
text(y2,x2,['B','(',num2str(x2),',',num2str(y2),')'])