Rotations en coordonnées sphériques

Aide à la résolution d'exercices ou de problèmes de niveau Supérieur.

Modérateur: gdm_aidesco

Règles du forum
Merci d'éviter le style SMS dans vos messages et de penser à utiliser la fonction Recherche avant de poster un message. Pour joindre des fichiers à vos messages, consulter ce sujet.
> Penser à utiliser le mode LaTeX (voir ici) afin de rendre vos formules plus lisibles.
> Ne poster qu'un exercice (ou problème) par sujet et indiquer son niveau précis dans le titre du message.

Rotations en coordonnées sphériques

Messagepar FiReTiTi » Mardi 23 Décembre 2014, 16:45

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 [url=http://fr.wikipedia.org/wiki/Coordonnées_sphériques]données ici[/url] :
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 ?
FiReTiTi
FiReTiTi
Hecto-utilisateur
 
Messages: 63
Inscription: Mercredi 26 Juillet 2006, 12:49
Localisation: Portland, OR, USA
Statut actuel: Actif et salarié

Publicité

Re: Rotations en coordonnées sphériques

Messagepar Framboise » Dimanche 28 Décembre 2014, 00:19

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 ?
Framboise
Téra-utilisateur
 
Messages: 1151
Inscription: Lundi 21 Mai 2007, 12:57
Localisation: Dordogne
Statut actuel: Post-bac | Doctorat

Re: Rotations en coordonnées sphériques

Messagepar FiReTiTi » Lundi 29 Décembre 2014, 00:43

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 ?
FiReTiTi
FiReTiTi
Hecto-utilisateur
 
Messages: 63
Inscription: Mercredi 26 Juillet 2006, 12:49
Localisation: Portland, OR, USA
Statut actuel: Actif et salarié

Re: Rotations en coordonnées sphériques

Messagepar Framboise » Lundi 29 Décembre 2014, 11:49

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 édition par Framboise le Lundi 29 Décembre 2014, 14:46, édité 1 fois.
J'ai le virus des sciences, ça se soigne ?
Framboise
Téra-utilisateur
 
Messages: 1151
Inscription: Lundi 21 Mai 2007, 12:57
Localisation: Dordogne
Statut actuel: Post-bac | Doctorat

Re: Rotations en coordonnées sphériques

Messagepar FiReTiTi » Lundi 29 Décembre 2014, 12:27

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.
FiReTiTi
FiReTiTi
Hecto-utilisateur
 
Messages: 63
Inscription: Mercredi 26 Juillet 2006, 12:49
Localisation: Portland, OR, USA
Statut actuel: Actif et salarié

Re: Rotations en coordonnées sphériques

Messagepar Framboise » Lundi 29 Décembre 2014, 15:21

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 ?
Framboise
Téra-utilisateur
 
Messages: 1151
Inscription: Lundi 21 Mai 2007, 12:57
Localisation: Dordogne
Statut actuel: Post-bac | Doctorat

Re: Rotations en coordonnées sphériques

Messagepar FiReTiTi » Mardi 30 Décembre 2014, 01:43

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.
FiReTiTi
FiReTiTi
Hecto-utilisateur
 
Messages: 63
Inscription: Mercredi 26 Juillet 2006, 12:49
Localisation: Portland, OR, USA
Statut actuel: Actif et salarié

Re: Rotations en coordonnées sphériques

Messagepar Framboise » Mardi 30 Décembre 2014, 14:02

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 ?
Framboise
Téra-utilisateur
 
Messages: 1151
Inscription: Lundi 21 Mai 2007, 12:57
Localisation: Dordogne
Statut actuel: Post-bac | Doctorat

Re: Rotations en coordonnées sphériques

Messagepar FiReTiTi » Mercredi 31 Décembre 2014, 01:36

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_trigonometric_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
FiReTiTi
Hecto-utilisateur
 
Messages: 63
Inscription: Mercredi 26 Juillet 2006, 12:49
Localisation: Portland, OR, USA
Statut actuel: Actif et salarié

Re: Rotations en coordonnées sphériques

Messagepar FiReTiTi » Mercredi 31 Décembre 2014, 02:12

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.
FiReTiTi
FiReTiTi
Hecto-utilisateur
 
Messages: 63
Inscription: Mercredi 26 Juillet 2006, 12:49
Localisation: Portland, OR, USA
Statut actuel: Actif et salarié

Re: Rotations en coordonnées sphériques

Messagepar kojak » Mercredi 31 Décembre 2014, 10:08

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
kojak
Modérateur
 
Messages: 10304
Inscription: Samedi 18 Novembre 2006, 19:50
Statut actuel: Actif et salarié | Enseignant

Re: Rotations en coordonnées sphériques

Messagepar Framboise » Mercredi 31 Décembre 2014, 11:54

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 ?
Framboise
Téra-utilisateur
 
Messages: 1151
Inscription: Lundi 21 Mai 2007, 12:57
Localisation: Dordogne
Statut actuel: Post-bac | Doctorat

Re: Rotations en coordonnées sphériques

Messagepar kojak » Mercredi 31 Décembre 2014, 16:38

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
kojak
Modérateur
 
Messages: 10304
Inscription: Samedi 18 Novembre 2006, 19:50
Statut actuel: Actif et salarié | Enseignant

Re: Rotations en coordonnées sphériques

Messagepar Framboise » Mercredi 31 Décembre 2014, 16:59

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 ?
Framboise
Téra-utilisateur
 
Messages: 1151
Inscription: Lundi 21 Mai 2007, 12:57
Localisation: Dordogne
Statut actuel: Post-bac | Doctorat

Re: Rotations en coordonnées sphériques

Messagepar FiReTiTi » Vendredi 02 Janvier 2015, 17:03

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
FiReTiTi
Hecto-utilisateur
 
Messages: 63
Inscription: Mercredi 26 Juillet 2006, 12:49
Localisation: Portland, OR, USA
Statut actuel: Actif et salarié

Re: Rotations en coordonnées sphériques

Messagepar FiReTiTi » Vendredi 02 Janvier 2015, 18:12

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.
FiReTiTi
FiReTiTi
Hecto-utilisateur
 
Messages: 63
Inscription: Mercredi 26 Juillet 2006, 12:49
Localisation: Portland, OR, USA
Statut actuel: Actif et salarié


Retourner vers Exercices et problèmes : Supérieur

 


  • Articles en relation
    Réponses
    Vus
    Dernier message

Qui est en ligne

Utilisateurs parcourant ce forum: Aucun utilisateur enregistré et 2 invités

cron