Calculation of spherical distance between two points of spherical surface (2)
1. Calculate
The theoretical basis of the calculation of two spherical distances is mentioned in the above, in which this is calculated.
Similarly, hypotheses of point A is east, latitude is north latitude B1 degrees; the length of point B is east, latitude is the South latitude B2 degrees.
Therefore, α = B1, β = B2, the value of the ∠ACE is θ, then θ = | A1 - A2 |.
Can be calculated
| OC | = r * sin (α) | α | = r * sin (β)
| AC | = r * COS (α) | BD | = r * COS (β)
Due to | CE | = | BD |, | CD | = | OD | | OC |
From the cosine theorem, you can know
| AE | 2 = | AC | 2 | CE | 2 - 2 * | AC | * | CE | * COS (θ)
ð | AE | 2 = R2 * COS2 (α) R2 * COS2 (β) - 2 * R2 * COS (α) COS (β) * COS (θ)
| Be | = | CD |
From | AB | 2 = | AE | 2 | BE | 2
ð | AB | 2 = R2 * COS2 (α) R2 * COS2 (β) - 2 * R2 * COS (α) COS (β) * COS (θ) R2 * SIN2 (α) R2 * SIN2 (β) ) 2 * r2 * sin (α) sin (β)
ð | AB | 2 = 2R2 (1 sin (α) sin (β) - COS2 (α) Cos2 (β) COS (θ))
In triangular AOB, you can also use the cosine to get
ð | OA | 2 | OB | 2- 2 | OA | * | OB | * COS (AOB) = | OB | 2
ð COS (AOB) = COS (α) * COS (β) * COS (θ) - sin (α) sin (β)
This makes it possible to obtain the value of AOB (in an arc), and the last spherical distance is R * (AOB).
Ok, huh, huh.
Note that
(1) When A, b is the same as the same east, θ = | A1 - A2 |, then when A, B is not the same as the east, θ = | A1 A2 |, if θ> 180 degrees, θ = 360 -θ (assuming unit is angle, not curvature).
(2) When A, B is the Southern Hemisphere or the same as the northern hemisphere, | CD | = | | OC | - | OD | | | OD |.
2. Program implementation
(1) Function interface definition
INT CALGLOBEDISTANCE (const st_globepnt a_stpnta,
Const ST_GLOBEPNT A_STPNTB,
Const Double A_Dradius,
Double & a_ddistance)
return value:
0: Perform OK, return the value correctly;
1: A_stpnta or A_STPNTB East, Western Value is incorrect;
2: A_stpnta or the South latitude of A_STPNTB, the north latitude is incorrect;
3: A_STPNTA or A_STPNTB's longitude range;
4: A_stpnta or a_stpntb latitude range offline; 10: The radius value is incorrect (not less than 0).
parameter:
A_STPNTA
The coordinates of a given first point, please refer to the following ST_GLOBEPNT definitions for structural ST_GLOBEPNT;
A_STPNTB
The coordinates of a given second point, please refer to the following ST_GLOBEPNT definitions for structural ST_GLOBEPNT;
A_DRADIUS
The radius of a given sphere cannot be less than 0;
a_ddistance
The last obtained spherical distance, if an error occurs during the calculation process, its value is negative.
Structural ST_GLOBEPNT
Typedef struct
{
BYTE BTEORW; // Take the value of 1, East Economic; 0, West;
Byte btnors; // Take the value of 1, north latitude; 0, south latitude;
Double Dlong; // longitude, value range [0.0, 180.0)
Double DLAT; // latitude, value range [0.0, 90.0]
} ST_GLOBEPNT;
(2) program code
Constant definition
Typedef struct
{
Char bteorw; // Take the value of 1, East (0, West;
Char btnors; // Take the value of 1, north latitude; 0, south latitude;
Double Dlong; // longitude, value range [0.0, 180.0)
Double DLAT; // latitude, value range [0.0, 90.0]
} ST_GLOBEPNT;
Const char east = 1;
Const char West = 0;
Const char north = 1;
Const char number = 0;
Const int exc_ok = 0; // 0: Perform OK, return value correctly;
Const int eorw_error = 1; // 1: East, Western Value;
Const int NORS_ERROR = 2; // 2: South latitude, north latitude.
Const int long_error = 3; // 3: longitude range offline;
Const int lat_error = 4; // 4: latitude range offline;
Const int RADIUS_ERROR = 10; // 10: The radius value (cannot be less than 0).
Const Double Accuracy = 0.00001; // Accuracy
Const Double Pi = 3.14159265359;
Function implementation
Auxiliary function int checkgpnt (const st_globepnt & a_stpnt)
Whether the function judgment point meets the requirements, the implementation is as follows:
INT Checkgpnt (const st_globepnt & a_stpnt)
{
INT IRET = EXC_OK;
// Check if East or West
IF (a_stpnt.bteorw! = EAST && A_STPNT.BTEORW! = WEST)
{
IRET = EORW_ERROR;
Return IRet;
}
// Check if Nort or South
IF (a_stpnt.btnors! = north && a_stpnt.btnors! = South) {
IRET = EORW_ERROR;
Return IRet;
}
// Check the Range of Long.
IF (a_stpnt.dlong - 180.0> Accuracy || A_STPNT.DLONG <= - Accuracy)
{
IRET = long_error;
Return IRet;
}
// Check the Range of Lat
IF (a_stpnt.dlat - 90.0> Accuracy || A_STPNT.DLAT <= -AccuracY)
{
Iret = lat_error;
Return IRet;
}
Return IRet;
}
The interface functions are as follows:
INT CALGLOBEDISTANCE (const st_globepnt & a_stpnta,
Const ST_GLOBEPNT & A_STPNTB,
Const Double & A_Dradius,
Double & a_ddistance)
{
INT IRET = EXC_OK;
A_DDISTANCE = -1.0;
Double dtheta; // the Angle of the Two longs
Double dalpha; // the Angle of a_stpnta's lat
Double dbeta; // the Angle of a_stpntb's lat
Double Dresult; // the Angle of the Two Lines
// Which One is through the given point and the center
Double dcosResult; // = COS (DRESULT);
// Check if The Given Points IS OK - Begin -------
IRET = CheckGPNT (A_STPNTA);
IF (IRET! = EXC_OK)
{
Return IRet;
}
IRET = CheckGPNT (A_STPNTB);
IF (IRET! = EXC_OK)
{
Return IRet;
}
// Check if the given points is ok - Nend ---------
// Check the Range of Radius
IF (A_DRADIUS <-Accurac "
{
IRET = RADIUS_ERROR;
Return IRet;
}
Dalpha = a_stpnta.dlat;
Dbeta = a_stpntb.dlat;
// Calculate the the Angle of the two longs --- begin --------
IF (A_STPNTA.BTEORW == a_stpntb.bteorw)
{
DTheta = FABS (a_stpnta.dlong - a_stpntb.dlong);
} else //! =
{
DTheta = a_stpnta.dlong a_stpntb.dlong;
}
IF (dtheta - 180.00> -accurac ")
{
DTheta = 360.0 - DTHETA;
}
// Calculate the the Angle of the two longs --- End --------
// Change Degree to Radian --- Begin ------ Dalpha = DALPHA / 180.0 * PI;
Dbeta = dbeta / 180.0 * pi;
DTheta = dtheta / 180.0 * pi;
// Change Degree to Radian --- End ------------
// Calculate the Angle of the Two Lines --- Begin ---------
IF (a_stpnta.btnors! = a_stpntb.btnors)
{
DCOSRESULT = COS (Dalpha) * COS (DBETA) * COS (DTHETA) - SIN (Dalpha) * sin (Dbeta);
} else //! =
{
DcosResult = COS (Dalpha) * COS (Dbeta) * COS (DTHETA) SIN (Dalpha) * sin (Dbeta);
}
Dresult = acos (dcosResult); // Dresult Ranges [0, Pi], Needn't Change
// Calculate the angle of the two --- ----------
// Get the distance arround the globe
A_DDISTANCE = a_dradius * Dresult;
Return IRet;
}
(3) Supplementary description
Using the program calculation, there may be certain errors, and the method of reducing the error may be used to calculate the algorithm mentioned in the calculation method, call, first, and not consider this.
3. About this article
This article is
Version 1.2
,in
2007
year
5
month
8
Day, in
1.1
Based on the basis of modification.
Version 1.1
Delivery
2005
year
1
month
12
day.
You can find the original text and related source code in this article "http://zijinshi.cn (website server can be needed to access) to find the original text and related source code, of course, you can also from my blog: http: // This article found on blog.9cbs.net/zijinshi.
My email: frank@zijinshi.cn, welcome to enlighten me. Bamboo
(Finish)
