使用MATLAB计算北京到纽约的新旧航线里程

数学实验的PPT课件参见百度文库
http://wenku.baidu.com/view/2eda008884868762caaed5a6.html
本文给出了实验的示例程序以及输出结果。

以下是mlab21.m的代码，用于计算北京到纽约的航线例程：

	Beijing=[40 116];
	Shanghai=[31 122];
	Tokyo=[36 140];
	SanFrancisco=[37 -123];
	NewYork=[41 -76];
	city=[Beijing;Shanghai;Tokyo;SanFrancisco;NewYork];
	R=6400+10;
	theta=city(:,1)*pi/180; %第一列，即纬度
	fai=city(:,2)*pi/180; %第二列，即经度
	x=R*cos(theta).*cos(fai);
	y=R*cos(theta).*sin(fai);
	z=R*sin(theta);
	op=[x,y,z];
	format bank
	Dmatrix=R*acos(op*op'/R.^2)

view raw 2011040601.m hosted with ❤ by GitHub

执行mlab21后得到以下结果

Dmatrix =
             0       1144.90       2155.64       9608.53      10993.77
       1144.90             0       1766.78       9936.44      11870.27
       2155.64       1766.78             0       8283.06      10764.64
       9608.53       9936.44       8283.06             0       4061.47
      10993.77      11870.27      10764.64       4061.47             0

可以清楚的看到，北京直飞纽约的航线例程为10993.77km，
而旧航线里程为1144.90+1766.78+8283.06+4061.47= 15256.21 (km).
假设飞机按照统计的平均时速980km/h匀速前进，那么新航线比旧航线至少节约了4.35小时，这尚未包含在各个中转站花费的候机时间。

以下是mlab22.m的代码，用于在球面上绘制新旧航线：

	Beijing=[40 116];
	Shanghai=[31 122];
	Tokyo=[36 140];
	SanFrancisco=[37 -123];
	NewYork=[41 -76];
	city=[Beijing;Shanghai;Tokyo;SanFrancisco;NewYork];
	R=6400+10;
	theta=city(:,1)*pi/180; %第一列，即纬度
	fai=city(:,2)*pi/180; %第二列，即经度
	x=R*cos(theta).*cos(fai);
	y=R*cos(theta).*sin(fai);
	z=R*sin(theta);
	op=[x,y,z];
	format bank
	Dmatrix=R*acos(op*op'/R.^2)

view raw 2011040602.m hosted with ❤ by GitHub

这一段程序调用了skyway(p1,p2,color)函数，代码如下：

	function skyway(p1,p2,color)
	city=[p1;p2];
	theta=city(:,1)*pi/180;fai=city(:,2)*pi/180;
	x=cos(theta).*cos(fai);
	y=cos(theta).*sin(fai);
	z=sin(theta);
	t=linspace(0,1,50);
	xt=(1-t)*x(1)+t*x(2);
	yt=(1-t)*y(1)+t*y(2);
	zt=(1-t)*z(1)+t*z(2);
	r=sqrt(xt.^2+yt.^2+zt.^2);
	xt=xt./r;yt=yt./r;zt=zt./r;
	figure(1),
	plot3(x,y,z,'ro',xt,yt,zt,color,'linewidth',2);

view raw 2011040603.m hosted with ❤ by GitHub

执行mlab.22之后，可以在Figure(1)窗口中看到绘制的新旧航线图：

视角一(示太平洋)

视角二(示北极)