J'ai bien mon idée sur le problème, mais la formule finale en fonction des latitudes et longitudes est extrêmement compliqué...
Posons O le centre de la Terre.
notons a,b,c les vecteurs normalisé des vecteur OA,OB et OC
1) On cherche un vecteur normal pour chacun des plans (OAB) et (OAC) :
n1 = a^b et n2 = a^c ("^" est le produit vectoriel, et ce sont des vecteurs)
2) On a :
||n1^n2|| = ||n1||x||n2|| sin(Â)
n1.n2 = ||n1||x||n2|| cos(Â)
Là on peut utiliser Arcos où Arcsin, mais je préfère utiliser atan ce qui simplifira les ||n1||x||n2|| (j'appelle atan(x,y) la fonction qui calcul l'angle en utilisant w = Arctan(x/y) mais en tenant compte du signe de x et y pour renvoyer w, w-Pi, ou w-Pi)
 = atan(||n1^n2|| , n1.n2)
3) Simplifions:
||n1^n2|| = ||(a^b)^ (a^c)|| = || (a.(a^b))c-(c.(a^b))b|| = || det(a,b,c)b|| = det(a,b,c) car ||b||=1
Je rappelle aussi que par définition (u^v).w=det(u,v,w)
n1.n2 = (a^b). (a^c) = det(a,b,a^c)
4) On obtient la belle formule  = atan( det(a,b,c) , det(a,b,a^c) )
5) On exprime les vecteurs a,b,c en coordonnées cartésienne :
On pose "ta" la latitude en radian de A, "la" la longitude en radian de A, idem pour B et C.
a = (cos(ta)sin(la) , cos(ta)cos(la), sin(ta) )
b = (cos(tb)sin(lb) , cos(tb)cos(lb), sin(tb) )
a = (cos(tc)sin(lc) , cos(tc)cos(lc), sin(tc) )
6) On utilise un logiciel de calcul formel pour avoir la formule moche
 = atan(cos(lc) cos(tc) (cos(tb) sin(lb) sin(ta) - cos(ta) sin(la) sin(tb)) +
cos(lb) cos(tb) (-cos(tc) sin(lc) sin(ta) + cos(ta) sin(la) sin(tc)) +
cos(la) cos(ta) (cos(tc) sin(lc) sin(tb) - cos(tb) sin(lb) sin(tc))) ,
(cos(ta)^2 cos(tb) cos(tc) (cos(lb) sin(la) - cos(la) sin(lb)) (cos(lc) sin(la) -
cos(la) sin(lc)) + (cos(lb) cos(tb) sin(ta) -
cos(la) cos(ta) sin(tb)) (cos(lc) cos(tc) sin(ta) -
cos(la) cos(ta) sin(tc)) + (cos(tb) sin(lb) sin(ta) -
cos(ta) sin(la) sin(tb)) (cos(tc) sin(lc) sin(ta) -
cos(ta) sin(la) sin(tc)))