Spherical geometry question - directions and non-great circles

The earth being oblong/elliptical makes things difficult for flat maps.

On one app, using flat Leaflet maps, I usually show one thing (a city, port, etc…) and below that have a navigation row of buttons where the user can go to the next place and so forth. I use nautical miles and try to show a compass direction. I reckon the distance is about right yet the compass pointer is really only a general “go thataways”. I use what I believe is the haversine algorithm through PostGis as such:

                ST_DistanceSphere(
                    t.geom_point,
                    (SELECT p.geom_point FROM ports p WHERE p.id = target_id)
                ) / 1852 AS distance_nm,
                ROUND(
                    degrees(
                        ST_Azimuth(
                            (SELECT p.geom_point FROM ports p WHERE p.id = target_id),
                            t.geom_point
                        )
                    )
                ) AS degrees

It’s a horseshoes and handgrenades closeness.

Recently I started using spherical globes with THREE and globe.gl. My goal was to have airplanes of various kinds fly from airport to airport. I try to keep the airplane pointing in the direction it’s going and rather than keep doing diffs I calculate the curve at takeoff.

My flightcurve function:

    function makeFlightCurve(startLat, startLng, endLat, endLng) {
        const start = globe.getCoords(startLat, startLng);
        const end = globe.getCoords(endLat, endLng);
        const startVec = new THREE.Vector3(start.x, start.y, start.z);
        const endVec = new THREE.Vector3(end.x, end.y, end.z);
        const midVec = startVec.clone().add(endVec).normalize();
        const earthRadius = startVec.length();
        const distance = startVec.distanceTo(endVec);
        const arcHeight = Math.min(earthRadius * 0.2, distance * 0.3);
        midVec.multiplyScalar(earthRadius + arcHeight);
        return new THREE.QuadraticBezierCurve3(startVec, midVec, endVec);
    }

All the planes appear to be flying too high, and furthermore for long-haul flights I have to give them even more altitude in the middle else they dive into the Pacific and emerge near Hong Kong (etc..) so they often look like they’re in low-earth-orbit.

I’ve put both of the above on hiatus for a while to get other useful things done.

Next I wanted to find out the bearing of an object moving on a general great circle would be - which can be determined by calculating the dot product of the velocity vector and the northward vector along the Earth’s surface. Call that vector N. N starts at the point (x,y,z) and ends at (0, 0, z-R^2/z) (picture a right triangle with vertices at (0, 0, 0) and (x, y, z) - the third point is at (0, 0, R/sin(L), but we can simplify that).

So the bearing is the arccos of ((x,y,z) dot N)/(|N||(x,y,z)|).

Working all the algebra gets bearing = arc cos (sin(Lm)/sin(L) (1-sin^2(L)/sin^2(Lm))^(1/2)).

More later

In my “great circle” app, it would just take a few parameters (coordinates and height) and bam - great circle.

In the app I described above, I wanted to animate the same for an airplane. In my previous example, I rely heavily on established JavaScript functions - ending with the “simple” request for a QuadraticBezierCurve3.

My goal is (was?) to make the planes fly close enough to a (not so) spherical earth that they a) did not dip below the surface and b) look like they were in orbit. My tweaks and such for distance were okay enough, yet we’ll be flying at FL 100,000 between New York and Tokyo.

ETA: Not a pilot, yet I believe that would still be an extraordinary FL 100 which maybe a U2 could achieve

But the equation for bearing isn’t what I need - I need to know the maximum latitude of the great circle that goes through a given latitude Ld at a given bearing Bd.

So rearrange the equation to get Lm=arc cos (cos Ld sin Bd).

The last piece - what’s the time when the object goes through Ld at bearing Bd?

Using the previous equations, I can get to wt=arc cos (sin(Ld)/sin(Lm))

Getting there…

So now we want to figure out the derivatives of L(t) and l(t) at time t=arc cos (sin(Ld)/sin(Lm))/w

There’s a lot of algebra there, but the answers come out simple-ish

L’(arc cos (sin(Ld)/sin(Lm))/w) = w cos(Bd)
L’'(arc cos (sin(Ld)/sin(Lm))/w) =w^2 sin^2(Bd) tan (Ld)

On a rhumb line course,
L’(arc cos (sin(Ld)/sin(Lm))/w) = w cos(Bd) (same as on great circle), but
L’'(t) = 0 for all t.

So w^2 sin^2(Bd) tan (Ld) is the northward acceleration required to keep on the rhumb line (keep deviating from great circle (shortest-distance) path.

This makes sense - you can see that it’s right by imagining an airplane circling the north pole at high latitude - there’s the ordinary centripetal acceleration required to keep going in a circle, and the tan L factor accounts for the fact that we’re interested in the component of that centripetal acceleration that’s along the Earth’s surface.