Trigo sphérique - Navigation

[Sextant]

Bonjour,
Pouvez-vous m'aider à résoudre un problème de navigation ? :
Soient 2 cercles tracés à la surface d'une sphère (la terre), et dont on connait pour chacun d'eux :
- les coordonnées du centre (latitude / longitude)
- la longueur angulaire du rayon.
Si d'autre part la somme des 2 rayons est supérieure à l'écart entre les 2 centres, comment calculer les coordonnées des 2 poins d'intersection des cercles ?
Merci de votre aide.

[Dzp]

Bonjour,
La logique des formules en trigonométrie sphérique est très voisine que celle utilisée pour la trigonométrie plane.
Si les rayons sont trop petits, alors les cercles ne se coupent pas et il n'y a pas d'intersection.

[Robot]

On travaille dans un repère orthonormé direct d'origine le centre de la sphère (supposée de rayon 1) avec pour le point de latitude et longitude et pour le pôle nord. Les coordonnées du point de la sphère de longitude et latitude sont .
Soient et les coordonnées des centres des cercles dessinés sur la sphère, et leurs rayons angulaires respectifs. On a à résoudre le système

Les deux dernières équations permettent d'exprimer deux des inconnues en fonction de la troisième, et on obtient alors en portant ces expressions dans la première équation une équation du second degré, qui peut ou non avoir des solutions réelles.
Une fois déterminées les solutions , on peut revenir aux longitudes et latitudes par les formules données plus haut, avec , etc.

[Dzp]

Bonjour Robot,
Je ne doute pas que cette méthode que tu indiques ne fournisse le résultat, mais apparemment, ce n'est pas cela qu'attend Sextant.
On a un triangle sphérique dont on connait 3 côtés (la distance entre les 2 centres de cercle et les deux rayons). On a toute une panoplie de formules pour résoudre ce genre de chose.

[Robot]

Rien ne t'empêche de proposer une formule de la panoplie qui donne la solution du problème de Sextant. On pourra alors juger sur pièces, et voir en particulier si c'est vraiment différent de ce que je propose.

[Dzp]

Bonjour,
[HS]J'avais fait une réponse, mais j'ai probablement fait une fausse manipulation.[/HS]
J'ai fait un essai.
Soit A le centre d'un cercle situé au voisinage du village d'Astérix, et C le centre d'un cercle situé près de Calais.
Le trésor est situé près de Rouen. J'ai calculé les rayons des cercle pour être dans les hypothèses de base.
La précision du résultat est de l'ordre de 1/1000 de degré, ce qui n'est pas très bon.

J'ai commencé le calcul par la méthode cartésienne donnée par Robot, et j'ai renoncé, devant la complexité des calculs.
Il existe une méthode dite du point approché. Elle consiste à prendre un point approché, par n'importe quel moyen, calculer à partir de ce point approché les deux rayons, puis on fait un graphique à grande échelle et on applique les corrections.

[Ben314]

J'ai commencé le calcul par la méthode cartésienne donnée par Robot, et j'ai renoncé, devant la complexité des calculs.

??????
Ah, ça c'est sûr, résoudre un système linéaire puis une équation du second degrés, c'est "incroyablement complexe" (surtout à l'époque où tout est informatisé...)

[Dzp]

Bonjour Ben,
Je ne vois pas très bien ce que vient faire l'informatique dans l'histoire. L'informatique, dans le cas précis ne sert qu'à faire les opérations.
Le formule consiste à calculer x² et y². Ces deux termes sont des carrés de 4 termes
x²=(b'.cosp + b'cz - b.cosp' + bc'z)² / DC²
La forme de ces relations (des différences) plus la résolution de l'équation du second degré va énormément dégrader la précision.
Avec la méthode trigonométrique qui ne présente pas cet inconvénient, on a déjà une perte de précision importante.
Il est bien évident que dans le cas général, je n'aurais pas hésité à faire le calcul. ce n'est pas la complexité en elle même qui m'a arrêté, mais la perte de précision due aux types d'opération.
Autrefois, il y avait un exercice classique dont le titre était "rendre calculable par logarithme". J'ai été confronté à ce problème (et j'en ai encore honte), il y a quelques mois, il s'agissait de calculer un arc de cercle en connaissant sa "courbure".

D'ailleurs, j'ai dans ma machine un module qui calcule les solutions d'un système linéaire, sauf une équation qui est du second degré. J'ai un autre module qui calcule une solution d'un système de N équations du second degré. Alors, tu vois, c'est pas les calculs qui me font peur. Evidemment je ne fais appel à aucun logiciel extérieur, type Wolframe ou formel.
Ceci dit, je n'ai pas très bien compris l'intérêt de ton message.

[Ben314]

Ceci dit, je n'ai pas très bien compris l'intérêt de ton message.

Uniquement à signaler que, contrairement à ce que tu affirme, les calculs en question sont on ne peut plus simple.
Evidement, il faut avoir un truc pour faire les calculs (calculette ou ordi), mais de toute façon, comme plus personne n'utilise de table pour lire dedans les valeurs des fonctions trigo, tu ne coupera pas à l'utilisation d'une calculette.

Et sinon, je peut te retourner la question : où voit tu ne serait ce qu'un tout petit peu de "complexité" dans la résolution d'un bête système linéaire et d'une équation du second degré ?

[Dzp]

Bon, je vais les faire, et je donnerai les résultats avec explications nécessaires. Rassure-toi, sauf cas particulier, je ne fais plus mes calculs à la table de log, mais avec mon ordinateur.
La complexité réside dans le fait qu'on soustrait des nombres de même ordre. Donc on perd toute la précision. Je pensais avoir été clair dans mon message précédent.

[Dzp]

Bonjour,
J'ai fait le calcul en mode cartésien. Voila les données de base :
/*Intersection de deux cercles
Asterix Lat=49°37'36"551 Lon=-1°37'19"801
Calais Lat=50°57'40"767 Lon=1°44'49"926
Rouen Lat=49°28'25"545 Lon=1°02'33"772
*/
Avec la méthode cartésienne, j'arrive à l'équation du second degré suivante :
z²*((zB*yA - zA*yB)²/DC² + (xB*zA - xA*zB)²/DC² +1) +
2*z *((zB*yA - zA*yB)*(rA*yB - rB*yA)/DC² + (xB*zA - xA*zB)*(xA*rB - xB*rA)/DC²) +
((rA*yB - rB*yA)² + (xA*rB - xB*rA)²)/DC²) -1 = 0
delta = b'² - a*c
z1=(-b' + sqrt(delta))/a
z2=(-b' - sqrt(delta))/a
On remarque que les coefficients sont constitués de produits, puis de différences, et on élève tout ça au carré.
Par ailleurs, on sait que l'informatique ne sait que résoudre très rapidement des opérations que l'on peut faire à la main. En l'occurrence, l'informatique ne change par vraiment le problème de fond qui est celui de la précision des calculs qui est directement liée au nombre de chiffres significatifs.
Sauf erreur de ma part, cette méthode donne une précision de l'ordre du 1/10 de degré.

[Ben314]

Sauf erreur de ma part, cette méthode donne une précision de l'ordre du 1/10 de degré.

Je peut te rassurer : il y a très clairement ENORME erreur de ta part : avec la calculatrice/ordi la plus pourrie qui soit, tu as une précision infiniment plus grande que le dixième de degré.

Je ferais bien le calcul EXACT de l'erreur commise par une machine pour te démontrer que tu es (de beaucoup) dans le faux, sauf que je vois pas le rapport entre tes données initiales, à savoir 3 points et la question posée qui demande de connaitre deux points et un (ou deux) rayon(s).

[Dzp]

Pour ton premier paragraphe, j'ai fait le calcul en double précision, c'est à dire les nombres sont sur 64 bits.
Pour le résultat de ton calcul, tu n'oublieras pas de préciser le nombre de chiffres significatifs utilisé par ta machine.
Pour le second paragraphe, le seul moyen de vérifier les calculs est de fixer le résultat, puis calculer les deux distances angulaires, celles-ci servant de base au calcul de l'intersection de cercles. On est donc dans l'hypothèse initiale, on connait les deux points (Astérix et Calais) et on a calculé des deux rayons. Il me parait évident que le problème n'est défini que si on connait DEUX rayons et non un seul.

[Ben314]

Tapé à la va vite et sans aucun test (delta<0 etc...) ça donne ça :

Code: Tout sélectionner

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct coord coord;
struct coord
{ double x;
  double y;
  double z;
};

double Rayon(coord A, coord B)
{ return acos(A.x*B.x+A.y*B.y+A.z*B.z);
}

coord Coordonnees(double latD, double latM, double latS, double lonD, double lonM, double lonS)
{ double lat,lon; coord cM;
  lat=(latD+latM/60+latS/3600)*2*M_PI/360; lon=(lonD+lonM/60+lonS/3600)*2*M_PI/360;
  cM.x=cos(lat)*cos(lon); cM.y=cos(lat)*sin(lon); cM.z=sin(lat);
  return cM;
} 

void Solve(coord A, double rA, coord B, double rB, coord *M1, coord *M2)
{ double det,coefX,coefY,cstX,cstY,a,b,c,delta,z1,z2;
  det=A.x*B.y-B.x*A.y;
  coefX=(-A.z*B.y+B.z*A.y)/det; cstX=(cos(rA)*B.y-cos(rB)*A.y)/det;
  coefY=(-A.x*B.z+B.x*A.z)/det; cstY=(A.x*cos(rB)-B.x*cos(rA))/det;
  a=coefX*coefX+coefY*coefY+1; b=coefX*cstX+coefY*cstY; c=cstX*cstX+cstY*cstY-1;
  delta=b*b-a*c; z1=(-b-sqrt(delta))/a; z2=(-b+sqrt(delta))/a;
  M1->x=coefX*z1+cstX; M1->y=coefY*z1+cstY; M1->z=z1;
  M2->x=coefX*z2+cstX; M2->y=coefY*z2+cstY; M2->z=z2;
}

void Affiche(coord M)
{ double lat,lon; int latD,latM,lonD,lonM;
  lat=asin(M.z)*360/(2*M_PI); lon=atan(M.y/M.x)*360/(2*M_PI);
  latD=(int)lat; lat=60*(lat-latD); lonD=(int)lon; lon=60*(lon-lonD);
  latM=(int)lat; lat=60*(lat-latM); lonM=(int)lon; lon=60*(lon-lonM); 
  printf("(%d:%d:%f ; %d:%d:%f)\n",latD,latM,lat,lonD,lonM,lon);
}

int main()
{ coord Asterix,Calais,Rouen,Point1,Point2; double rAsterix,rCalais;
  Asterix = Coordonnees( 49, 37, 36.551 , -1,-37,-19.801);
  Calais  = Coordonnees( 50, 57, 40.767 ,  1, 44, 49.926);
  Rouen   = Coordonnees( 49, 28, 25.545 ,  1, 02, 33.772);
  rAsterix = Rayon(Asterix,Rouen);
  rCalais  = Rayon(Calais ,Rouen);
  Solve(Asterix,rAsterix,Calais,rCalais,&Point1,&Point2);
  Affiche(Rouen);
  Affiche(Point1);
  Affiche(Point2);
  return 0;
}

Qui donne ça :
(49:28:25.545000 ; 1:2:33.772000)
(49:28:25.545000 ; 1:2:33.772000)
(51:15:22.607505 ; 0:-41:-1.056099)
C'est a dire le bon résultat avec une erreur inférieure au millionième de seconde d'arc

Et même si tu pousse le bouchon jusqu'à faire les mêmes calculs avec des réels simple précision (32 bits), ça donne ça :
(49:28:25.545044 ; 1:2:33.772118)
(49:28:26.053162 ; 1:2:33.439095)
(51:15:22.151184 ; 0:-41:-0.460052)
Donc une erreur de l'ordre de la demi seconde d'arc (soit environ 15 mètres)

[Robot]

Super, Ben. Juste un petit détail : j'aurais plutôt utilisé

Code: Tout sélectionner

lon=atan2(M.y, M.x)*360/(2*M_PI);

[Ben314]

De toute façon, c'est tout fait "à la pas propre" :
- Le système, pour pas perdre trop de précision surtout éviter les cas où l'intersection des deux plans est (quasi) parallèle au plan z=0 (i.e. le déterminant est quasi nul), il faudrait calculer les 3 déterminants possibles et prendre le plus grand.
- Je vérifie rien du tout concernant Delta donc ça plante s'il est négatif et ça renvoie soit un truc faux, soit un truc absolument pas précis du tout lorsque Delta est proche de 0 (le deux cercles sont quasi tangents et la moindre variation de rayon change grandement le résultat)
- Et effectivement, pour l'affichage, ça ne marche que si la longitude est entre -90° et 90° vu que j'ai pas pris l'arc moitié pour le calcul (je connaissait pas la fonction atan2...)
- Plus sans doute pas mal d'autres trucs à améliorer...

[Dzp]

Bonjour Ben,
Oui, bravo, mais il y a juste un petit problème concernant l'énoncé.
Tu calcules les rayons à partir des coordonnées cartésiennes de Rouen alors que justement, c'est cela que l'on doit calculer. J'ai indiqué les coordonnées géographiques pour être sûr des valeurs de base, Asterix et Calais connu, Rouen à calculer.
Dans l'énoncé, les rayons sont connus par leur angle géographique : " la longueur angulaire du rayon."
Sauf erreur de calcul (mais disons que ce sont les hypothèses de l'énoncé du problème), les longueurs angulaires des rayons sont :
RayAsrerix= 2.6692700° ; RayCalais = 1.6459133°
Je les ai donnés en degrés décimaux.

[Ben314]

Tu calcules les rayons à partir des coordonnées cartésiennes de Rouen alors que justement, c'est cela que l'on doit calculer.

Heuuuuuu.
Pour toi, ça signifie quoi le verbe CONNAITRE ?

Soient 2 cercles tracés à la surface d'une sphère (la terre), et dont on connait pour chacun d'eux :
- les coordonnées du centre (latitude / longitude)
- la longueur angulaire du rayon.

[Dzp]

Si on était dans le plan, le problème serait de calculer un point défini par les distances à deux points donnés. Il s'agit de l'intersection de deux cercles, ça se calcule par la relation des cos (a² = b² + c² - 2bc cosA)
Là on est sur une sphère (la terre assimilée à un sphère), on connait les deux points par leur coordonnées géographiques et les deux distances sont déterminées par les longueurs des deux arcs, ce que Sextant appelle la longueur angulaire. Le terme employé n'est très correct, mais, étant donné le contexte, il n'y a pas d'ambiguïté.
Sans autre précision, lorsqu'on parle de distance mesurée sur la terre, il s'agit d'un arc de grand cercle. qu'on appelle géodésique, et naturellement pas de la distance qui correspondrait à une corde. Si on veut calculer en géométrie plane, alors on fait intervenir la notion de projection, ce qui n'est pas le cas ici.

[Robot]

Dzp, tu n'as pas l'air d'avoir compris que, si on se donne deux points de coordonnées cartésiennes et sur la sphère unité (autrement dit et ), alors leur distance angulaire (= plus petite longueur d'arc de grand cercle qui les joint) est si et seulement si .