javascriptmathcoordinatesearthdistance

Find intersection coordinates of two circles on earth?


I'm trying to find a second intersection point of two circles. One of the points that I already know was used to calculate a distance and then used as the circle radius (exemple). The problem is that im not getting the know point, im getting two new coordinates, even thou they are similar. The problem is probably related to the earth curvature but I have searched for some solution and found nothing.

The circles radius are calculated with the earth curvature. And this is the code I have:

function GET_coordinates_of_circles(position1,r1, position2,r2) {
  var deg2rad = function (deg) { return deg * (Math.PI / 180); };
  x1=position1.lng;
  y1=position1.lat;
  x2=position2.lng;
  y2=position2.lat;
  var centerdx = deg2rad(x1 - x2); 
  var centerdy = deg2rad(y1 - y2); 
  var R = Math.sqrt(centerdx * centerdx + centerdy * centerdy);

  if (!(Math.abs(r1 - r2) <= R && R <= r1 + r2)) { // no intersection
    console.log("nope");
    return []; // empty list of results
  }

  // intersection(s) should exist
  var R2 = R*R;
  var R4 = R2*R2;
  var a = (r1*r1 - r2*r2) / (2 * R2);
  var r2r2 = (r1*r1 - r2*r2);
  var c = Math.sqrt(2 * (r1*r1 + r2*r2) / R2 - (r2r2 * r2r2) / R4 - 1);

  var fx = (x1+x2) / 2 + a * (x2 - x1);
  var gx = c * (y2 - y1) / 2;
  var ix1 = fx + gx;
  var ix2 = fx - gx;

  var fy = (y1+y2) / 2 + a * (y2 - y1);
  var gy = c * (x1 - x2) / 2;
  var iy1 = fy + gy;
  var iy2 = fy - gy;
  // note if gy == 0 and gx == 0 then the circles are tangent and there is only one solution
  // but that one solution will just be duplicated as the code is currently written
  return [[iy1, ix1], [iy2, ix2]];
}

The deg2rad variable it is suppose to adjust the other calculations with the earth curvature.

Thank you for any help.


Solution

  • Your calculations for R and so on are wrong because plane Pythagorean formula does not work for spherical trigonometry (for example - we can have triangle with all three right angles on the sphere!). Instead we should use special formulas. Some of them are taken from this page.

    At first find big circle arcs in radians for both radii using R = Earth radius = 6,371km

    a1 = r1 / R
    a2 = r2 / R
    

    And distance (again arc in radians) between circle center using haversine formula

    var R = 6371e3; // metres
    var φ1 = lat1.toRadians();
    var φ2 = lat2.toRadians();
    var Δφ = (lat2-lat1).toRadians();
    var Δλ = (lon2-lon1).toRadians();
    
    var a = Math.sin(Δφ/2) * Math.sin(Δφ/2) +
            Math.cos(φ1) * Math.cos(φ2) *
            Math.sin(Δλ/2) * Math.sin(Δλ/2);
    var ad = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
    

    And bearing from position1 to position 2:

     //where    φ1,λ1 is the start point, φ2,λ2 the end point 
     //(Δλ is the difference in longitude)
    var y = Math.sin(λ2-λ1) * Math.cos(φ2);
    var x = Math.cos(φ1)*Math.sin(φ2) -
            Math.sin(φ1)*Math.cos(φ2)*Math.cos(λ2-λ1);
    var brng = Math.atan2(y, x);
    

    Now look at the picture from my answer considering equal radii case. (Here circle radii might be distinct and we should use another approach to find needed arcs)

    enter image description here

    We have spherical right-angle triangles ACB and FCB (similar to plane case BD is perpendicular to AF in point C and BCA angle is right).
    Spherical Pythagorean theorem (from the book on sph. trig) says that

     cos(AB) = cos(BC) * cos(AC)
     cos(FB) = cos(BC) * cos(FC)
    

    or (using x for AC, y for BC and (ad-x) for FC)

     cos(a1) = cos(y) * cos(x)
     cos(a2) = cos(y) * cos(ad-x)
    

    divide equations to eliminate cos(y)

     cos(a1)*cos(ad-x) = cos(a2) * cos(x)
     cos(a1)*(cos(ad)*cos(x) + sin(ad)*sin(x)) = cos(a2) * cos(x)
     cos(ad)*cos(x) + sin(ad)*sin(x) = cos(a2) * cos(x) / cos(a1)
     sin(ad)*sin(x) = cos(a2) * cos(x) / cos(a1) - cos(ad)*cos(x)
     sin(ad)*sin(x) = cos(x) * (cos(a2) / cos(a1) - cos(ad))
     TAC = tg(x) = (cos(a2) / cos(a1) - cos(ad)) / sin(ad)
    

    Having hypotenuse and cathetus of ACB triangle we can find angle between AC and AB directions (Napier's rules for right spherical triangles) - note we already know TAC = tg(AC) and a1 = AB

    cos(CAB)= tg(AC) * ctg(AB)
    CAB = Math.acos(TAC * ctg(a1))
    

    Now we can calculate intersection points - they lie at arc distance a1 from position1 along bearings brng-CAB and brng+CAB

    B_bearing = brng - CAB
    D_bearing = brng + CAB
    

    Intersection points' coordinates:

    var latB = Math.asin( Math.sin(lat1)*Math.cos(a1) + 
                  Math.cos(lat1)*Math.sin(a1)*Math.cos(B_bearing) );
    var lonB = lon1.toRad() + Math.atan2(Math.sin(B_bearing)*Math.sin(a1)*Math.cos(lat1), 
                         Math.cos(a1)-Math.sin(lat1)*Math.sin(lat2));
    

    and the same for D_bearing

    latB, lonB are in radians