Helpers to get vertical and horizontal bearings given two GPS points with altitude data

bearing.js

function toRadian(num) {
    return num * (Math.PI / 180);
}

function toDegree(num) {
    return num * (180 / Math.PI);
}

// North is 0 degree, South is 180 degree
function getHorizontalBearing(fromLat, fromLon, toLat, toLon, currentBearing) {
    fromLat = toRadian(fromLat);
    fromLon = toRadian(fromLon);
    toLat = toRadian(toLat);
    toLon = toRadian(toLon);

    let dLon = toLon - fromLon;
    let x = Math.tan(toLat / 2 + Math.PI / 4);
    let y = Math.tan(fromLat / 2 + Math.PI / 4);
    let dPhi = Math.log(x / y);
    if (Math.abs(dLon) > Math.PI) {
        if (dLon > 0.0) {
            dLon = -(2 * Math.PI - dLon);
        } else {
            dLon = (2 * Math.PI + dLon);
        }
    }

    let targetBearing = (toDegree(Math.atan2(dLon, dPhi)) + 360) % 360;
    return targetBearing - currentBearing;
}

// Horizon is 0 degree, Up is 90 degree
function getVerticalBearing(fromLat, fromLon, fromAlt, toLat, toLon, toAlt, currentElevation) {
    fromLat = toRadian(fromLat);
    fromLon = toRadian(fromLon);
    toLat = toRadian(toLat);
    toLon = toRadian(toLon);

    let fromECEF = getECEF(fromLat, fromLon, fromAlt);
    let toECEF = getECEF(toLat, toLon, toAlt);
    let deltaECEF = getDeltaECEF(fromECEF, toECEF);

    let d = (fromECEF[0] * deltaECEF[0] + fromECEF[1] * deltaECEF[1] + fromECEF[2] * deltaECEF[2]);
    let a = ((fromECEF[0] * fromECEF[0]) + (fromECEF[1] * fromECEF[1]) + (fromECEF[2] * fromECEF[2]));
    let b = ((deltaECEF[0] * deltaECEF[0]) + (deltaECEF[2] * deltaECEF[2]) + (deltaECEF[2] * deltaECEF[2]));
    let elevation = toDegree(Math.acos(d / Math.sqrt(a * b)));
    elevation = 90 - elevation;

    return elevation - currentElevation;
}

function getDeltaECEF(from, to) {
    let X = to[0] - from[0];
    let Y = to[1] - from[1];
    let Z = to[2] - from[2];

    return [X, Y, Z];
}

function getECEF(lat, lon, alt) {
    let radius = 6378137;
    let flatteningDenom = 298.257223563;
    let flattening = 0.003352811;
    let polarRadius = 6356752.312106893;
    
    let asqr = radius * radius;
    let bsqr = polarRadius * polarRadius;
    let e = Math.sqrt((asqr-bsqr)/asqr);
    // let eprime = Math.sqrt((asqr-bsqr)/bsqr);

    let N = getN(radius, e, lat);
    let ratio = (bsqr / asqr);

    let X = (N + alt) * Math.cos(lat) * Math.cos(lon);
    let Y = (N + alt) * Math.cos(lat) * Math.sin(lon);
    let Z = (ratio * N + alt) * Math.sin(lat);

    return [X, Y, Z];
}

function getN(a, e, latitude) {
    let sinlatitude = Math.sin(latitude);
    let denom = Math.sqrt(1 - e * e * sinlatitude * sinlatitude);
    return a / denom;
}

let n = getHorizontalBearing(39.099912, -94.581213, 39.099912, -94.588032, 0.00); // should be 270 degree
console.info("Horizontal bearing:\t", n);

let m = getVerticalBearing(39.099912, -94.581213, 273.543, 39.099912, -94.588032, 873.543, 0.0); // should be around 40-50 degree
console.info("Vertical bearing:\t", m);