lat1 et lon1 sont les coordonnées du premier point
lat2 et lon2 sont celle du deuxième point, j'assume que ces valeurs sont converties en radians
On veut la distance et le cap de 1 vers 2
Il nous faut dlon = lon1 - lon2, mais là il y a un piège subtil qui arrive si les distances sont très grandes (je sais je suis tombé dedans

Il faut limiter dlon entre -PI et + PI. Sur ma page web je le fais avec la fonction normalise() que voici:
/***************************************************************
* Function to convert a difference of radians to +/- pi
* The result is equal to x but from -pi to +pi
* Function added for Version 1.3
****************************************************************/
function normalise(x){
x = x % (2*pi) ; // difference modulo 2*pi
if (x < 0) x += 2*pi; // rendre l'angle positif
if (x > pi) x -= 2*pi; // Set dlon between -pi to +pi regime
return x;
}
donc dlon = normalise(lon1 - lon2)
J'ai aussi les signes des latitudes comme variables
las1 = 1 si la latitude du point 1 est au nord et -1 si elle est au sud
las2 = 1 si la latitude du point 2 est au nord et -1 si elle est au sud
Ensuite si lat1 n'est pas égal à lat2:
cap = atan( dlon /( las1*acosh(1/cos(lat1)) - las2*acosh(1/cos(lat2))))
et comme la tangente a deux possibilités vérifier s'il faut ajouter PI
:
if (lat2 < lat1) cap += pi;
la distance est:
dist = 60*180 /pi * ABS((lat2 -lat1)/cos(cap)) // distance en milles nautiques
Mais si lat2 == lat1 on fait autrement:
if (lon1 > lon2) {
cap = -pi/2 ;
} else {
cap = pi/2;
}
dist = 60*180 / pi * ABS( dlon * cos(lat1)) ; // en milles nautiques
Et finalement pour mettre le cap de 0 à 360
if ( cap < 0) cap += 2*pi
cap = cap * 180 /pi
Le piège dont je parlais est que pour joindre 2 points sur la projection Mercator avec une ligne droite (définition du loxodrome) on a plusieurs solutions (une infinité en fait) mais ce qu'on veut c'est la plus courte solution, or mon calculateur marchait très bien jusqu'au jour où quelqu'un a voulu traverser l'antiméridien du Japon jusqu'à San Francisco et mon programme l'a envoyé à travers la Russie. C'était mathématiquement correct mais pas ce qu'il voulait. D'où la fonction normalise() que j'ai ajouté en 2013.