sos!!求个距离算法。

[yl_yz]

由两点的经、纬度和方向怎样计算出这两点的实际距离？3ks！

[letheanwater]

两个点的经纬度都知道了.
可以确定一个矩形,求这个矩形的对角线,不就可以了吗?

[yl_yz]

to fengyun925(深秋的落叶 卷起片片的凌乱) ：照你说的，更糊涂了。

[banner91]

你是用什么软件啊，普通算术是不可以用的
MAPXTREME用
bRC = CalcMapDistance(dblMapX1, dblMapY1, dblMapX2, dblMapY2, dblMapDistance) ''计算距离

就可以了

[banner91]

你用什么软件得？MAPX和MAPXTREME一样得啊，你再帮助文档搜索一下这个函数把，其他软件我就不知道了

[荒原独歌]

网上流行的三个方法，经实测，EarthDis误差较小。这是我改成的Delphi版本，VB的和VC的可以在CSDN上用关键字“距离”搜。
=====================================
unit Unit2;

interface

uses
Math;

function Distance(lon1: Double; lat1: Double; lon2: Double; lat2: Double): Double;
function GetDistance(lan1,lon1, lan2, lon2: Double) : Double;
function EarthDis(lan1, lon1, lan2, lon2: Double): Double;

implementation
const
R_A = 6378137;
R_B = 6356752;
//WGS-84 ³¤Öá°ë¾¶
EARTH_WGS84_A = 6378137.0000; //WGS-84Semi-major axis in meters
//WGS-84 ¶ÌÖá°ë¾¶
EARTH_WGS84_B = 6356752.3142; //WGS-84 Semi-minor axis in meters
//WGS-84 Æ«ÐÄÂʵÄƽ·½
EARTH_WGS84_E2 = 0.00669437999013; //WGS-84 Eccentricity squared
EARTH_WGS84_MPM = 1852.0; //meters per nmile
//WGS-84 µØÇò±âÂÊ
EARTH_WGS84_FLATTENING = 298.257223563; //WGS-84

function Distance(lon1: Double; lat1: Double; lon2: Double; lat2: Double): Double;
var
a_2d: Double;
e_2d: Double;
h_2d: Integer;
DEG_2_RAD: Double;
RAD_2_DEG: Double;

x_rads: Double;
y_rads: Double;
n_2ds: Double;
x_2d: Double;
y_2d: Double;
z_2d: Double;
x_radm: Double;
y_radm: Double;
n_2dm: Double;
x_2d_mark: Double;
y_2d_mark: Double;
z_2d_mark: Double;
curdistance: Double;
begin
a_2d := 6378137;
e_2d := 0.00669438;
h_2d := 15;
DEG_2_RAD := 0.01745329252;
RAD_2_DEG := 57.2957795129;

x_rads := Abs(lon1) * DEG_2_RAD;
y_rads := Abs(lat1) * DEG_2_RAD;


n_2ds := a_2d / Sqr(1 - e_2d * Sin(y_rads) * Sin(y_rads));

x_2d := (n_2ds + h_2d) * Cos(y_rads) * Cos(x_rads);
y_2d := (n_2ds + h_2d) * Cos(y_rads) * Sin(x_rads);
z_2d := (n_2ds * (1 - e_2d) + h_2d) * Sin(y_rads);

x_radm := Abs(lon2) * DEG_2_RAD;
y_radm := Abs(lat2) * DEG_2_RAD;

n_2dm := a_2d / Sqr(1 - e_2d * Sin(y_radm) * Sin(y_radm));

x_2d_mark := (n_2dm + h_2d) * Cos(y_radm) * Cos(x_radm);
y_2d_mark := (n_2dm + h_2d) * Cos(y_radm) * Sin(x_radm);
z_2d_mark := (n_2dm * (1 - e_2d) + h_2d) * Sin(y_radm);

curdistance := (x_2d_mark - x_2d) * (x_2d_mark - x_2d) + (y_2d_mark - y_2d) * (y_2d_mark - y_2d) + (z_2d_mark - z_2d) * (z_2d_mark - z_2d);
curdistance := Sqr(curdistance);

Result := curdistance;
end;


{'------------------------------------------
Private Const PI = 3.1415926 'Pi

Private m_R As Double

Public Enum Distance_Mode '
Miles = 0
Feet = 1
End Enum

Private Type PointType
lat As Double
lon As Double
End Type

Private m_PointA As PointType
Private m_PointB As PointType
Private m_DistMode As Distance_Mode
'to get the parameters of DistanceMode
Public Property Get Dist_Mode() As Distance_Mode
Dist_Mode = m_DistMode
End Property

Public Property Let Dist_Mode(ByVal vNewValue As Distance_Mode)
m_DistMode = vNewValue
End Property

Private Sub Class_Initialize()
m_R = (R_A + R_B) / 2
End Sub
'get the first point
Public Function SetPointA(lat As Double, lon As Double) As Boolean
If (lat <= 90 And lat >= -90) And (lon >= -180 And lon <= 180) Then
m_PointA.lat = lat
m_PointA.lon = lon
SetPointA = True
Else
SetPointA = False
End If
End Function
'get the second point
Public Function SetPointB(lat As Double, lon As Double) As Boolean
If (lat <= 90 And lat >= -90) And (lon >= -180 And lon <= 180) Then
m_PointB.lat = lat
m_PointB.lon = lon
SetPointB = True
Else
SetPointB = False
End If
End Function }
//calculate the distance
function GetDistance(lan1,lon1, lan2, lon2: Double) : Double;
var
angAB: Double;
tempV: Double;
tempD: Double ;

angLatA: Double;
angLatB: Double;
angLonA: Double;
angLonB: Double;
m_R : Double;
itempV : Integer;
begin

m_R := (R_A + R_B) / 2;

angLatA := lan1 / 180 * PI;
angLonA := lon1 / 180 * PI;
angLatB := lan2 / 180 * PI;
angLonB := lon2 / 180 * PI;

tempV := Cos(angLatA) * Cos(angLatB) * Cos(angLonB - angLonA) + Sin(angLatA) * Sin(angLatB);
if (tempV <> 1) and (tempV <> -1) then
begin
angAB := ArcTan(-tempV /Sqr(-tempV * tempV + 1)) + 2 * ArcTan(1) ;
tempD := m_R * angAB;
end;
//angAB := 0
if tempV = 1 then tempD := 0;

//angAB := 180
if tempV = -1 then tempD := m_R * PI ;

{if m_DistMode = Miles then
tempD := tempD / 1609
else
tempD := tempD / 1609 * 5280; }

Result := tempD;

end;

function EarthDis(lan1, lon1, lan2, lon2: Double): Double;
var
D: double ;
tmpVal: double ;
fi1,fi2: double ;
drda: double;
tmp1,tmp2: double ;
begin
drda := lon2 - lon1;
//drda := drda/60.0;// in degrees
drda := drda*PI/180.0;// in radians
fi1 := lan1;
fi2 := lan2;

fi1 := fi1 * PI / 180.0;//in in radians
fi2 := fi2 * PI / 180.0;//in in radians

tmpVal:=sin(fi1)*sin(fi2)+cos(fi1)*cos(fi2)*cos(drda);
{if(fabs(tmpVal)>1.0) then
begin
AfxMessageBox("Invalidate input valve of arccos!",MB_ICONSTOP);
return;
end; }
D:=arccos(tmpVal);// in radians
//dis:=(dis*180.0/PI)*60.0; // in minutes

tmp1:=(sin(fi1)+sin(fi2));
tmp2:=(sin(fi1)-sin(fi2));
tmpVal := ( (3*sin(D) - D)*tmp1*tmp1 ) / ( 1+cos(D) );
tmpVal := tmpVal - ( (3*sin(D) + D)*tmp2*tmp2 ) / ( 1-cos(D) );

Result :=EARTH_WGS84_A*D + (EARTH_WGS84_A/(4*EARTH_WGS84_FLATTENING)) * tmpVal;//in meters
//Result :=Result/EARTH_WGS84_MPM ;//in nmiles

//double C=(EARTH_WGS84_A-EARTH_WGS84_B)/EARTH_WGS84_A;
//double C1=1.0/C;
end;



end.