Hallar la intersección de dos círculos dada la Lat/Lon y el radio?

Question

Estoy intentando calcular la intersección de dos círculos en la Tierra con una latitud, longitud y radio determinados. He empezado con esto ¿Calcular la intersección de dos círculos?

Estos círculos nunca tendrán más de 100 km de diámetro y normalmente estarán "cerca" el uno del otro, es decir, dos puntos de intersección para el par de círculos.

He añadido un comentario en el código de abajo en el lugar donde estoy teniendo el problema.

Originalmente hice esta pregunta en StackOverflow, pero siendo esto tan geoespacial, no obtuve ninguna respuesta.

var circle1 = {lat:37.673442, lon:-90.234036, r:107.5};
var circle2 = {lat:36.109997, lon:-90.953669, r:145};

function findInterscetionOfTwoCircles(circle1, circle2){
    var c1 = latLonToGeocentricCoords(circle1);
    var c2 = latLonToGeocentricCoords(circle2);
    var r1 = convertRadius(circle1.r);
    var r2 = convertRadius(circle2.r);
    var c1dotc2 = calculateDotProduct(c1, c2);

    // This is where the problem is occurring.
    // The value of 'a' should be 0.9735030 but I am getting 46.770178816522844.
    // I am certain that r1 and r2 are correct (match the values in the original post). Additionally, I believe the calculateDotProduct function to be
    // correct, therefore making c1dotc2 correct. I think I have ruled out
    // any Order of Operations error. Maybe it is a degree/radian conversion
    // issue or something like it? I tried different variations and I  and now out
    // of ideas.
    var a = (Math.cos(r1) - (Math.cos(r2) * c1dotc2)) / (1 - (c1dotc2 * c1dotc2)); 
    var b = (Math.cos(r2) - (Math.cos(r1) * c1dotc2)) / (1 - (c1dotc2 * c1dotc2)); 
    // -45.83805010474024, should be  0.0260194
    var n = calculateCrossProduct(c1, c2); 

    var x0 = calculateLinearCombination(a, b, c1, c2);

    var ndotn = calculateDotProduct(n, n);

    // var t = Math.sqrt((1 - calculateDotProduct(x0,x0))/calculateDotProduct(n,n));

    // Finish the algorithm
}

function latLonToGeocentricCoords(circle){
    var x = Math.cos(degToRad(circle.lon)) * Math.cos(degToRad(circle.lat));
    var y = Math.sin(degToRad(circle.lon)) * Math.cos(degToRad(circle.lat));
    var z = Math.sin(degToRad(circle.lat));

    return {x:x, y:y, z:z};
}

function geocentricCoordsToLatLon(x, y, z){
    var lon = Math.atan2(x,y)
    var lat = Math.atan2(Math.sqrt((x*x)+(y*y)), z)

    return { lat: lat, lon: lon};
}

function convertRadius(radius){
    // converts NM right now
    return degToRad(radius) / 60; 
}

function calculateDotProduct(input1, input2){
    if(arguments.length !== 2){
        throw new Error('You must supply two arguments to this method.');
    }

    if(Array.isArray(input1) !== Array.isArray(input2) || typeof input1 !== typeof input2){
        throw new Error('Inputs must be the same types.');
    }

    var dotProduct = 0;
    if(Array.isArray(input1)){
        // process as parallel arrays
        if(input1.length === input2.length){
            throw new Error('Input lengths must be the same.');
        }
        for(var i = 0, len = input1.length; i < len; i++){
            dotProduct += input1[i] * input2[i];
        }
    } else if(typeof input1 === 'object' && !Array.isArray(input1)){
        // process as key-value object
        for(var key in input1){
            if (input1.hasOwnProperty(key)) {
                if(!input2[key]){
                    throw new Error('Both inputs must have the same properties to be processed.');
                }
                dotProduct += input1[key] * input2[key];
            }

        }
    }

    return dotProduct; 
}

function calculateCrossProduct(circle1, circle2){
    var x = (circle1.y * circle2.z) - (circle1.z * circle2.y);
    var y = (circle1.z * circle2.x) - (circle1.x * circle2.z);
    var z = (circle1.x * circle2.y) - (circle1.y * circle2.x);

    return {x: x, y: y, z: z};
}

function calculateLinearCombination(a, b, c1, c2){
    var x = {x: a * c1.x, y: a * c1.y, z: a * c1.z};
    var y = {x: b * c2.x, y: b * c2.y, z: b * c2.z};

    return calculateDotProduct(x, y);
} 

function degToRad(degree){
    return degree * Math.PI / 180;
}

function radToDeg(radian){
    return radian / Math.PI * 180;
}

Answer 1

El círculo1 con r=107,5 millas náuticas y el círculo2 con r=145 millas náuticas se cruzan en los puntos que se muestran en la imagen:

Para este tipo de trabajo hay que adoptar el sistema de coordenadas proyectadas (cartesianas). He trabajado en la zona UTM 15 porque ambos centros del círculo se encuentran en esta zona. Si se intenta hacer esto manualmente/matemáticamente utilizando el sistema de coordenadas polares, será mucho más complicado debido a la forma imperfecta de la tierra. Se puede realizar un cálculo más preciso adoptando el mejor sistema de proyección/coordenadas para esta zona.