[abc886y365]

R=6400;
theta=(-90:15:90)*pi/180; %纬度
fai=(-180:15:180)*pi/180; %经度
X=R*cos(theta)'*cos(fai);
Y=R*cos(theta)'*sin(fai);
Z=R*sin(theta)'*ones(size(fai)); %将原坐标转化为直角坐标
mesh(X,Y,Z)
colormap([0 0 0])
axis off
hold on
city=[40,116,31,122,36,140,37,-123,41,-76;31,122,36,140,37,-123,41,-76,40,116];
d=0;
s=0;
for k=1:2:9
theta=city(:,k)*pi/180;
fai=city(:,k+1)*pi/180;
x=R*cos(theta).*cos(fai);
y=R*cos(theta).*sin(fai);
z=R*sin(theta);
x1=x(1);y1=y(1);z1=z(1);
x2=x(2);y2=y(2);z2=z(2); %提取两个直角坐标点
Op1=[x1,y1,z1];Op2=[x2,y2,z2];
dk=R*acos(Op1*Op2'/R^2); %求球面距离
s=s+dk;
t=(0:20)/20;
X=(1-t)*x1+t*x2;
Y=(1-t)*y1+t*y2;
Z=(1-t)*z1+t*z2; %计算连线上各点的坐标
r=sqrt(X.*X+Y.*Y+Z.*Z);
X=R*X./r;Y=R*Y./r;Z=R*Z./r;
plot3(x,y,z,'r<',X,Y,Z,'b','LineWidth',2)
hold on
end
s
view(-93,68)