Theory
We can think of each polyline as a line segment that represents a plane defined by the end points of the line segment and the center of the Earth. The point of intersection would be the point where two planes intersects which would create a line that passes through the center of the Earth, thus representing a point on the surface of the Earth.
Figure 1 Intersection of two planes defined by line segments
We will look at the solution when Cartesian coordinates are used. A line segment can be defined by two points P1 and P2. Each point having an x, y and z component. We can define a plane that contains P1, P2 and the center of the Earth (P0) taking the dot product of P0 and the cross product of P1 and P2 which will equal zero.
We can evaluate this as follows:
We can represent a second line segment the same way which consists of points P3, and P4.
We can then solve for x and Y in terms of Z as follows:
The point of intersection with this line and the sphere of radius r has z such that the distance from the center of the Earth is r. Thus;
From this we can solve for z.
Knowing z we can find x and y, however the equations are starting to get fairly large. We can significantly reduce the number of calculations that need to performed to solve x, y and z by first solving sub components. For instance, let us define the following subcomponents.
We can then write simplified versions of x, y, and z:
We could then convert this to work directly with spherical coordinates; however the equations grow in size fairly quickly. It is better to convert the coordinates to Cartesian coordinates and then perform our calculations and then convert the points back into spherical coordinates.
Since the Earth is Spherical there will be two points on the surface of the Earth where these planes will intersect. This second point is the inverse coordinate.
Application
To calculate the intersection of two planes we have to define the planes with line segments. A line segment will be defined by two VELatLong objects. Since we want to find to intersection of two planes we need to pass four VELatLong objects into our function to represent the two line segments. We then have to convert these VELatLong objects into Cartesian coordinates. From there we can calculate two vectors that represent the two planes. We can then calculate the unit vectors of each plane and check that the planes are not equal. If the planes are equal we will return -1. If the planes are not equal we can calculate the vector director. The unit vector of the vector director is the first point of intersection. We can take the first intersection point and calculate its inverse coordinate to calculate the second point of intersection. We then have to convert our x, y, and z parameters into a VELatLong objects. We can then return an array containing these two VELatLongs. Note that this algorithm will calculate the point of intersections to within a distance of 1% of the largest polyline length used to define a plane.
function intersectionOfPlanes(latlong1,latlong2,latlong3,latlong4)
{
var point1 = convertSphericalToCartesian(latlong1);
var point2 = convertSphericalToCartesian(latlong2);
var point3 = convertSphericalToCartesian(latlong3);
var point4 = convertSphericalToCartesian(latlong4);
var vector1 = crossProduct(point1,point2);
var vector2 = crossProduct(point3,point4);
var unitVector1 = unitVector(vector1);
var unitVector2 = unitVector(vector2);
if((unitVector1.X != unitVector2.X )&& (unitVector1.Y != unitVector2.Y) &&
(unitVector1.Z != unitVector2.Z))
{
var vectorDirector = crossProduct(unitVector1,unitVector2);
var intersection1 = unitVector(vectorDirector);
var intersection2 = inverseCartesian(intersection1);
return new Array(convertCartesianToSpherical(intersection1),
convertCartesianToSpherical(intersection2));
}
return -1
}
Listing 1 Intersection Point of two Planes
Information on the convertSphericalToCartesian and convertCartesianToSpherical methods can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!280.entry
Information on the crossProduct method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!365.entry
Information on the unitVector method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!374.entry
Information on the inverseCartesian method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!342.entry