The computation the "hard" way can be simplified for the case r1 = r2 =: r. We still first have to convert the circle centers P1,P2 from (lat,lng) to Cartesian coordinates (x,y,z).
var DEG2RAD = Math.PI/180;
function LatLng2Cartesian(lat_deg,lng_deg)
{
var lat_rad = lat_deg*DEG2RAD;
var lng_rad = lng_deg*DEG2RAD;
var cos_lat = Math.cos(lat_rad);
return {x: Math.cos(lng_rad)*cos_lat,
y: Math.sin(lng_rad)*cos_lat,
z: Math.sin(lat_rad)};
}
var P1 = LatLng2Cartesian(lat1, lng1);
var P2 = LatLng2Cartesian(lat2, lng2);
But the intersection line of the planes holding the circles can be computed more easily. Let d
be the distance of the actual circle center (in the plane) to the corresponding point P1 or P2 on the surface. A simple derivation shows (with R the earth's radius):
var R = 6371; // earth radius in km
var r = 100; // the direct distance (in km) of the given points to the intersections points
// if the value rs for the distance along the surface is known, it has to be converted:
// var r = 2*R*Math.sin(rs/(2*R*Math.PI));
var d = r*r/(2*R);
Now let S1 and S2 be the intersections points and S their mid-point. With s = |OS|
and t = |SS1| = |SS2|
(where O = (0,0,0) is the earth's center) we get from simple derivations:
var a = Math.acos(P1.x*P2.x + P1.y*P2.y + P1.z*P2.z); // the angle P1OP2
var s = (R-d)/Math.cos(a/2);
var t = Math.sqrt(R*R - s*s);
Now since r1 = r2
the points S, S1, S2 are in the mid-plane between P1 and P2. For v_s = OS
we get:
function vecLen(v)
{ return Math.sqrt(v.x*v.x + v.y*v.y + v.z*v.z); }
function vecScale(scale,v)
{ return {x: scale*v.x, y: scale*v.y, z: scale*v.z}; }
var v = {x: P1.x+P2.x, y: P1.y+P2.y, z:P1.z+P2.z}; // P1+P2 is in the middle of OP1 and OP2
var S = vecScale(s/vecLen(v), v);
function crossProd(v1,v2)
{
return {x: v1.y*v2.z - v1.z*v2.y,
y: v1.z*v2.x - v1.x*v2.z,
z: v1.x*v2.y - v1.y*v2.x};
}
var n = crossProd(P1,P2); // normal vector to plane OP1P2 = vector along S1S2
var SS1 = vecScale(t/vecLen(n),n);
var S1 = {x: S.x+SS1.x, y: S.y+SS1.y, z: S.z+SS1.z}; // S + SS1
var S2 = {x: S.x-SS1.x, y: S.y-SS2.y, z: S.z-SS1.z}; // S - SS1
Finally we have to convert back to (lat,lng):
function Cartesian2LatLng(P)
{
var P_xy = {x: P.x, y:P.y, z:0}
return {lat: Math.atan2(P.y,P.x)/DEG2RAD, lng: Math.atan2(P.z,vecLen(P_xy))/DEG2RAD};
}
var S1_latlng = Cartesian2LatLng(S1);
var S2_latlng = Cartesian2LatLng(S2);