Rotations en coordonnées sphériques

Aide à la résolution d'exercices ou de problèmes de niveau supérieur au baccalauréat.
[participation réservée aux membres inscrits]
Règles du forum
Merci de soigner la rédaction de vos messages et de consulter ce sujet avant de poster. Pensez également à utiliser la fonction recherche du forum.
FiReTiTi

Rotations en coordonnées sphériques

Message non lu par FiReTiTi »

Bonjour,

je souhaite effectuer les rotations des éléments d'un tableau 3D dans un repère sphérique.
Pour cela, j'ai créé cette fonction à partir des formules données ici :

Code : Tout sélectionner

public static void SphericalRotations(double[][][] array, double theta, double phi, double[][][] result)
	{
	int x, y, z, px, py, pz ;
	int sizex = result[0][0].length ; // Dimensions du tableau résultat.
	int sizey = result[0].length ;
	int sizez = result.length ;
	int ami = array[0][0].length >> 1 ; // Centre du tableau d'origine
	int amy = array[0].length >> 1 ;
	int amz = array.length >> 1 ;
	int resmx = sizex >> 1 ; // Centre du tableau résultat
	int resmy = sizey >> 1 ;
	int resmz = sizez >> 1 ;
	double cx, cy, cz, st, sp, r ;
	
	
	for (z=0 ; z < sizez ; z++)
		for (y=0 ; y < sizey ; y++)
			for (x=0 ; x < sizex ; x++)
				{
				cx = x - resmx ; // Coordonnées cartésienne centrées.
				cy = y - resmy ;
				cz = z - resmz ;
				
				r = Math.sqrt(cx*cx + cy*cy + cz*cz) ; // Transformation en coordonnées sphériques.
				sp = Math.acos(cz / r) ;
				st = Math.atan2(cy, cx) ;
				
				sp -= phi ; // Rotations en sphérique.
				st -= theta ;
				
				cx = r * Math.sin(sp) * Math.cos(st) ; // Retour en cartésien.
				cy = r * Math.sin(sp) * Math.sin(st) ;
				cz = r * Math.cos(sp) ;
				
				px = ami + (int)(cx+0.5) ; // Vérifications pour éviter les débordements.
				if ( px < 0 || array[0][0].length <= px ) continue ;
				py = amy + (int)(cy+0.5) ;
				if ( py < 0 || array[0].length <= py ) continue ;
				pz = amz + (int)(cz+0.5) ;
				if ( pz < 0 || array.length <= pz ) continue ;
				
				result[z][y][x] = array[pz][py][px] ;
				}
	}
En résumé :
- pour chaque élément de mon tableau
- je calcule les coordonnées sphériques (depuis les coordonnées cartésiennes)
- je fais la rotation en coordonnées sphérique (simple soustraction)
- je repasse en coordonnées cartésiennes
- je copie dans mon résultat la case du tableau de départ correspondante.


Les rotations en fonction de theta (autour de l'axe Z) fonctionnent parfaitement, mais j'ai un résultat totalement aberrant pour la rotation avec phi (inclinaison par rapport à l'axe Z).
Est ce que quelqu'un pourrait m'expliquer mon erreur ?
Framboise
Utilisateur chevronné
Utilisateur chevronné
Messages : 1161
Inscription : lundi 21 mai 2007, 13:57
Statut actuel : Autre
Localisation : Dordogne

Re: Rotations en coordonnées sphériques

Message non lu par Framboise »

Bonjour,

Un peu de détails sur le problème encontré seraient bienvenus.
Le résultat faux, des distorsions, un manque de précision,... ?
Quel langage ( C++ ? ) et compilateur ?

Quelques remarques:
Les équations sources utilisées sont exactes en théorie mais sources de problèmes pour une application "pro" numérique.
On a des pertes de précision dans certaines plages trigo ( -> distorsions ).
Des plantages avec les fonctions ACOS, ASIN pour des arguments du style 1.000 000 000 000 000 1 pouvant parfois résulter des limites numériques de précision des calculs.

Il vaut mieux tout concevoir en passant par des coordonnées rect. -> polaire 3D -> rotation -> polaire 3D -> rect. effectivement.

Les fonctions ACOS, ASIN ( et ATAN, ASEC, ACSC ) ne doivent pas apparaitre sauf s'il fallait les utiliser strictement pour implémenter une fonction ATAN2 ( présente dans notre cas, donc sans ce problème ici si elle est bien conçue ).
Je vise le: sp = Math.acos(cz / r) ; incongru. Il faut choisir ensuite le bon quadrant vu que acos est ambigu sur ce point:
sp = Math.acos(cz / r) ;
sp = 2*PI - Math.acos(cz / r) ;
selon les cas.
Je préfère utiliser atan2 qui n'a pas d’ambiguïté et remplace avantageusement acos, pire encore si l'on rencontre acos (1.000 000 000 000 000 1) sans soucis pour atan2 non plus...

Quelques références:

Jean Meeus
Astronomical Algorithms

Numerical Recipes : The Art of Scientific Computing
http://www.numerical-recipes.com/
LA 2ed est dispo sur Internet en freeware par chapitre.
J'ai le virus des sciences, ça se soigne ?
FiReTiTi

Re: Rotations en coordonnées sphériques

Message non lu par FiReTiTi »

Bonjour, merci pour la réponse.

C'est un code Java avec une JVM 1.7

Pour les résultats aberrants, lorsque j'essaie d'incliner un segment, je me retrouve avec une sorte de V.

Tu écris "Il vaut mieux tout concevoir en passant par des coordonnées rect. -> polaire 3D -> rotation -> polaire 3D -> rect. effectivement.", mais il me semble que c'est ce que j'essaie de faire.

Comment atan2 peut remplacer ados ?
Framboise
Utilisateur chevronné
Utilisateur chevronné
Messages : 1161
Inscription : lundi 21 mai 2007, 13:57
Statut actuel : Autre
Localisation : Dordogne

Re: Rotations en coordonnées sphériques

Message non lu par Framboise »

je me retrouve avec une sorte de V
Ce serait bien un problème de quadrant qui n'est pas le bon selon le signe d'un des paramètres comme je le soupçonne plus haut.
Je ne vois rien dans le code qui traite ce problème du bon quadrant, donc selon les cas cela fonctionne ou provoque un déraillage ( le V ).

Pour passer de acos ou asin à atan2 : [ extrait d'un help fortran ]
ATAN2(Y, X) computes the principal value of the argument function of the complex number X + i Y. This function can be used to transform from Cartesian into polar coordinates and allows to determine the angle in the correct quadrant.
Un coup de magie avec Pythagore et l'on a les 2 côtés du triangle rectangle permettant d'appliquer atan2 au lieu de acos.
Cela évite les risques d'un acos ( 1.0....01 ) et l’ambiguïté de quadrant.
mais il me semble que c'est ce que j'essaie de faire.
Oui, je ne fais que confirmer que la voie est bonne [ d'où le 'effectivement' ].

Note: je ne pratique pas Java, mais le C qui en est très proche.
Dernière modification par Framboise le lundi 29 décembre 2014, 14:46, modifié 1 fois.
J'ai le virus des sciences, ça se soigne ?
FiReTiTi

Re: Rotations en coordonnées sphériques

Message non lu par FiReTiTi »

Merci pour ta réponse.
Je viens donc de remplacer

Code : Tout sélectionner

sp = Math.acos(cz / r)
par :

Code : Tout sélectionner

double d = cz / r ;
sp = Math.atan2(Math.sqrt(1-d*d), d) ;
J'ai malheureusement exactement le même "mauvais" résultat :-(
Je continue à regarder pour ces histoires de quadrants.
Framboise
Utilisateur chevronné
Utilisateur chevronné
Messages : 1161
Inscription : lundi 21 mai 2007, 13:57
Statut actuel : Autre
Localisation : Dordogne

Re: Rotations en coordonnées sphériques

Message non lu par Framboise »

Il semble que l'on perde en route l'info du bon quadrant.

Je vais regarder cela.

Un aperçu de mon code en C:

Code : Tout sélectionner

void RP3 ( double *x, double *y, double *z )
{ /* Conv. RECT(x,y,z)->POL(rho,a,d) 3-dim */
   RP( x, y );
   RP( x, z );
   return;
}

void RP( double *x, double *y )
{ /* RECT(X,Y) => POL(RHO,THETA) */
   double t;
   t = ANG( *x, *y );
   *x = RHO( *x, *y );
   *y = t;
   return;
}

double RHO( double x, double y )
{ /* Distance (Pythagore) */
   return sqrt( x * x + y * y );
}
avec ANG ( x, y ) correspondant à atan2 ( y, x )

C'est un fragment de programme à usage astro de calculs de planètes.
Les fonctions RP sont la transposition de cette fonction des calculatrices en RPN ( HP41,... ).
J'ai le virus des sciences, ça se soigne ?
FiReTiTi

Re: Rotations en coordonnées sphériques

Message non lu par FiReTiTi »

Il semblerait bien que ce soit un problème de quadrant.
Si j'ajoute

Code : Tout sélectionner

if ( cx < 0.0 ) sp = 2.0*Math.PI - sp ;
, alors cela règle le problème dans certains cas, mais pas tous.
Framboise
Utilisateur chevronné
Utilisateur chevronné
Messages : 1161
Inscription : lundi 21 mai 2007, 13:57
Statut actuel : Autre
Localisation : Dordogne

Re: Rotations en coordonnées sphériques

Message non lu par Framboise »

cx ... ou cz < 0.0 ?

Mais le problème ne devrait pas apparaitre avec atan2.
C'est :
sp = Math.acos(cz / r) ;
initial qui ne va pas.
sp = atan2( sqrt(cx*cx +cy*cy) , cz )
me semble la solution.

J'ai imaginé la figure:
axes Ox, Oy, Oz
Le point M, sa projection M' sur le plan xOy
cx, cy, cz les proj. sur les axes xyz.
sp = angle zOM "distance zénitale"
st = angle xOM'

donc atan2 doit prendre le triangle rect. cz O M' et non pas cz O M.

Sans la figure cela reste abscons.
J'ai le virus des sciences, ça se soigne ?
FiReTiTi

Re: Rotations en coordonnées sphériques

Message non lu par FiReTiTi »

Merci pour ton aide, mais malheureusement cela ne fonctionne pas :-(

Tu écris : "donc atan2 doit prendre le triangle rect. cz O M' et non pas cz O M."
Mais si cz est la projection de M sur l'axe Z et M' la projection de M sur le plan xOy, alors cet angle est toujours 90°.

En regardant ces formules ( http://en.wikipedia.org/wiki/Inverse_tr ... _functions ), j'avais plutôt conclu qu'il fallait remplacer
sp = Math.acos(cz / r)
par
double d = cz / r
sp = Math.atan2(Math.sqrt(1.0-d*d), d)
Mais ça ne fonctionne pas mieux :-(
FiReTiTi

Re: Rotations en coordonnées sphériques

Message non lu par FiReTiTi »

En fait, les trois formules sont identiques (je viens de tester) :

- 1 - double d = cz / r ;
phi = Math.atan2(Math.sqrt(1.0-d*d), d) ;

- 2 - phi = Math.atan2(Math.sqrt(cx*cx+cy*cy), cz) ;

- 3 - phi = Math.acos(cz / r) ;

Sauf que sur Wikipedia, la formule 3 est couplée avec "if (cy<0) phi = 2*PI-phi". Ce problème ne semble donc pas réglé en utilisant atan2.
De toute façon en ajoutant le test ou non, le résultat est exactement le même.
kojak
Modérateur général
Modérateur général
Messages : 10410
Inscription : samedi 18 novembre 2006, 19:50

Re: Rotations en coordonnées sphériques

Message non lu par kojak »

bonjour,

FiReTiTi a écrit : Sauf que sur Wikipedia, la formule 3 est couplée avec "if (cy<0) phi = 2*PI-phi".
Non, ce n'est pas l'angle $\varphi$ mais l'angle $\theta$.

Tu as bien $\varphi=\arccos\dfrac{z}{r}$ ce qui te donnera toujours un angle entre $0$ et $\pi$, alors qu'avec l'arctangente, tu auras n angle entre $]-\pi/2,\pi/2[$ donc ça ne colle pas.

Ton $\theta$ vaut bien $\arccos\dfrac{x}{\sqrt{x^2+y^2}$ si $y\geq 0$ car l'arccosinus te sort un angle entre $0$ et $\pi$ et l'autre $2\pi-\arccos\dfrac{x}{\sqrt{x^2+y^2}$ si $y\leq 0$. comme ça tu auras bien un angle entre $\pi$ et $2\pi$.

Pour finir avec ces coordonnées sphériques, il faut faire attention car il y a plusieurs définitions : celles des matheux et celle des physiciens. Là, j'ai pris ta référence du début.
Pas d'aide par MP.
Framboise
Utilisateur chevronné
Utilisateur chevronné
Messages : 1161
Inscription : lundi 21 mai 2007, 13:57
Statut actuel : Autre
Localisation : Dordogne

Re: Rotations en coordonnées sphériques

Message non lu par Framboise »

Mais si cz est la projection de M sur l'axe Z et M' la projection de M sur le plan xOy, alors cet angle est toujours 90°.
Le TR cz O M est bien rectangle, mais c'est l'angle cz O M ( non rect. dans le cas général ) qui nous intéresse et qui doit être fourni par atan2. cz O M' est lui un angle rect. effectivement.

Le choix des angle est un peu inhabituel, on utilise plutôt un système longitude - latitude, alors que là, sp est la distance au pôle N.

Ça change un peu des formules habituelles.

Pour info atan( y, x ) est un alias de atan2( y, x ) en C depuis quelques années.
2 - phi = Math.atan2(Math.sqrt(cx*cx+cy*cy), cz) ;

- 3 - phi = Math.acos(cz / r) ;
ne peuvent pas être identique car r = sqrt(cx*cx+cy*cy + cz*cz)

PS: Je ne trouve pas le Ç majuscule sur mon clavier QWERTZ-CH, j'ai bien le ç avec MAJ 4, mais rien en format majuscule. Une astuce autre que alt + pavé num. ?
En QWERTY-US int., c'est facile avec ' puis C, mais ici cela ne fonctionne pas.

Image
J'ai le virus des sciences, ça se soigne ?
kojak
Modérateur général
Modérateur général
Messages : 10410
Inscription : samedi 18 novembre 2006, 19:50

Re: Rotations en coordonnées sphériques

Message non lu par kojak »

Framboise a écrit :Le choix des angle est un peu inhabituel, on utilise plutôt un système longitude - latitude, alors que là
Tout dépend de quel point de vue on se place : matheux, physiciens ou géographes.

PS : sur ton image jointe, ton repère n'est pas direct, et ton cx n'est pas correct, car la projection de $M'$ sur les abscisses n'est pas parallèle à l’axe des ordonnées.
Pas d'aide par MP.
Framboise
Utilisateur chevronné
Utilisateur chevronné
Messages : 1161
Inscription : lundi 21 mai 2007, 13:57
Statut actuel : Autre
Localisation : Dordogne

Re: Rotations en coordonnées sphériques

Message non lu par Framboise »

Hello kojak,

Je me place d'un point de vue astronomique ascension droite et déclinaison mais j'ai préféré la terminologie longitude/latitude mieux connue et similaire.
C'est vrai, ma // est un peu de travers...
L'inversion "miroir": heureusement cela ne change pas le calcul.
Je n'ai pas d'excuse :roll: , je ne bois quasiment jamais d'alcool ( -> migraines :evil: ).

Bonnes fêtes à toi et à tout le monde
J'ai le virus des sciences, ça se soigne ?
FiReTiTi

Re: Rotations en coordonnées sphériques

Message non lu par FiReTiTi »

Le problème avec la rotation en Phi est en fait que ce n'est pas une rotation complète : un point est éloigné ou rapproché de l'origine Oz.
Donc actuellement, lorsque je fais une rotation suivant Phi, tout les points se rapprochent ou s'éloignent de l'axe, au lieu de tourner autour du point O.
Pour que ce soit une rotation, il faudrait que Phi soit codé sur [0, 2Pi[, mais je ne pense pas que ce soit possible.
FiReTiTi

Re: Rotations en coordonnées sphériques

Message non lu par FiReTiTi »

En fait, si c'est possible. Il faut juste que je fasse une rotation afin que le point à tourner soit dans un plan ( xOz par exemple ), une simple rotation autour de Y, puis retour dans le plan d'origine.