Discussion :

[LadyWasky]

Bonjour à tous, je suis un poil dans la panade, noyé par des Sinus et des Cosinus à n'en plus finir...

J'ai réalisé un petit logiciel de cartographie qui me permet de tracer n'importe quelle zone du Globe terrestre à partir des cordonnées géographiques des cotes, des frontières, des fleuves, etc...

La projection que j'ai utilisé est on ne peu plus simple : on appelle celà une projection "cylindrique equidistante" (voir ici si vous voulez approfondir : http://mathworld.wolfram.com/Cylindr...rojection.html)

En fait, c'est plus précisemment une projection Equirectangulaire étant donné que ma Latitude de référence est l'équateur (Phi=0).

Donc, bref et en résumé, voici le système d'équation que j'utilise pour passer de mes (longitudes,latitudes) à des coordonnées euclidiennes (x,y) :

x=Longitude=Lambda
y=Latitude=Phi

Maintenant, imaginez que sur ma carte, je souhaites tracer un cercle (donc en forme de cercle sur ma sphere, mais pas sur ma carte) de rayon 50km autour d'un point (Paris par exemple, que l'on représentera par un point ayant pour Longitude et Latitude (Lamda1,Phi1), donc sur ma carte (x1,y1)).
Pour avoi la distance sur une sphère entre deux points (Lamda1,Phi1) et (Lamda2,Phi2) on utilise les distances de deux points sur un "Grand cercle) de la sphère (l'Equateur et le méridiens sont des grands cercles) :
D'après l'équation (5) de cette page (attention ils ont représenté la latitude par la lettre grecque Delta cette fois ci)
http://mathworld.wolfram.com/GreatCircle.html

...j'obtiens comme équation pour mon "cercle" de centre (x1,y1) et de rayon (r) sur ma carte :

les inconnues :
x et y

les constantes :
r=50 km
a=6378 km (rayon terrestre)

C=cos(r/a)
A=cos(Phi1)
B=Sin(Phi1)

et la relation :
A.cos(y).cos(x1-x)+B.sin(y)=C
x, y, x1 et y1 sont en degrés

Mon problème est que je voudrais tracer ma courbe correspondant à mon cercle, non pas avec y en fonction de x (y=f(x), ni avec x en fonction de y (x=f(y)) qui m'obligerais à déterminer les bornes inférieures et supérieures de mes fonctions, mais plutot avec une fonction paramétrique du genre :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
x=R.cos(alpha)+x1
y=R.sin(alpha)+y1
où R est une fonction de x1,y1, A, B et C, et alpha avec alpha compris entre 0 et 2.Pi

Pour avoir y=f(x) ou x=f(y), pas de problèmes, mais à tracer sur la carte c'est galère.
pour le système d'équations paramétrées, là, je sèche. (et ce serait pourtant tellement plus simple à coder)

En fait, tou mon problème consiste à isoler R dans l'équation
A.cos(y).cos(x1-x)+B.sin(y)=C, ou on a remplacé x et y par leur formes paramétrées respective.

Je n'y arrive pas.

Y aurait-il un mathématcien dans la salle qui pourrais m'aider siouplé ?

Merci par avance.

[FrancisSourd]

Ce qui m'inquiète, c'est que lorsque ton point initial est proche d'un pôle (plus précisément lorsque le pôle est à l'intérieur du cercle), la courbe a une forme "bizarre" et je pense qu'il faut plusieurs formules en fonction de alpha.

Sinon, j'ai essayé de calculer l'équation... et je suis tombé sur la même formule que toi. J'aurais mieux fait de lire la page web directement;-)

[LadyWasky]

Salut FrancisSourd,

Je te remercie de t'être penché sur ce problème, c'est vraiment sympa

Lorsque le Pôle est à l'intérieur du cercle ou lorsqu'on s'approche des Pôles, je te rassure, c'est normal si la courbe à une forme bizarre
La projection d'un cercle de la sphère ne donnera jamais un cercle. Le cercle projeté sur la carte devrait en toute logique s'élargir dans sa partie la plus près du pôle (une forme de poire on va dire). Quand le Pôle est dans le cercle, alors on aura carrément une courbe partant du bord gauche de la carte et allant vers le bord droit.

Quand le pôle est sur le cercle, il faudra sans doute que je fasse attention, car les 2 points de la sphère que sont les pôles sont transformés en droites lors de la projection (bords haut et bas du planisphère).

Sinon, pour l'instant, j'en suis là :

deux équations :
x=x1+ArcCos( (C-B.sin(y)/(A.cos(y) )
et
x=x1-ArcCos( (C-B.sin(y)/(A.cos(y) )

(Supposant que ArcCos, le cosinus inverse, renvoie des degrés compris entre 0 et 180)

avec :
y1-r/a<y<y1+r/a
et y différent de 90° et -90°
et |C-sin(y)|<=|A.cos(y)|

Pour l'instant je fais avec ça, mais il est tout de même étonnant que l'on ne puisse pas obtenir une courbe paramétrée en fonction d'un rayon et d'un angle...

[FrancisSourd]

Tu peux essayer

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
y = arctan((C-A.cos(x1-x))/B)

car arctan à l'avantage d'être défini sur R
Erreur! Voir ci-dessous!

[LadyWasky]

Tu peux essayer

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
y = arctan((C-A.cos(x1-x))/B)

car arctan à l'avantage d'être défini sur R

C'est vrai, j'ai essayé, mais tu t'es trompé dans ton calcul.

A la place, ce donne plutot ça :

y = arctan((C-A.cos(x1-x))/(B.sin(y))

qui n'est pas franchement exploitable...

[FrancisSourd]

[Jean-Marc.Bourguet]

Cela fait un peu trop longtemps que j'ai fait de la trigonométrie pour pousser ce genre de calculs loin. Une idée: partir des équations paramètrique d'un cercle autou de lambda=0 phi = 0. On a

x = a cos r/a
y = a sin r/a cos t
z = a sin r/a sin t

et appliquer les rotations d'angles phi et lambda. Si je ne me trompe pas on arrive a

x = a cos r/a cos p cos l - a sin r/a sin t sin p cos l - a sin r/a cos t sin l
y = a cos r/a cos p sin l - a sin r/a sin t sin p sin l + a sin r/a cos t cos l
z = a cos r/a sin p + a sin r/a sin t cos p

Comme on veut en fait:

x = a cos p' cos l'
y = a cos p' sin l'
z = a sin p'

Donc

p' = arcsin(cos r/a sin p + sin r/a sin t cos p)
l' = arctan2((cos r/a cos p sin l - sin r/a sin t sin p sin l + sin r/a cos t cos l), (cos r/a cos p cos l - sin r/a sin t sin p cos l - sin r/a cos t sin l))

Il vaut mieux disposer d'une primitive comme atan2 en C qui utilise les signes du rapport.

[LadyWasky]

Merci, merci, c'est ce que je cherchais.

Je regarde ça de près et vous tient informé !

[LadyWasky]

Jean-Marc, tu es un géni : ça marche !



Bref les mots me manquent... you are the best.


un remerciement spécial s'impose pour Francis qui s'est lui aussi penché sur la question.

Merci à vous deux !

Pour info, Il y a bien Arctan2 en Delphi