Bonjour à tous,

J'ai besoin d'extraire les 3 angles d'Euler à partir d'une matrice pseudoinverse de rotation d'un vecteur a vers b. (comme défini dans la discussion : http://www.developpez.net/forums/d12...ntre-vecteurs/)

J'ai implémenté pour ça l'extraction de quaternions à partir de la matrice, puis la conversion quaternion-angle d'euler proposées dans l'article suivant : http://en.wikipedia.org/wiki/Rotatio...29#Quaternions :

quat.f
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
 PROGRAM QUAT
       implicit none
       character*80 cmd
       real :: M11,M12,M13,M21,M22,M23,M31,M32,M33
       real :: nA(3),nB(3),pi,q(4),q2(4),phi,theta,psi
       read(5,*) nA(1),nA(2),nA(3),nB(1),nB(2),nB(3)
       open(unit=20,file='scilab.sci',form='formatted',status='unknown')
       write(20,*) "vecA = [",nA(1),";",nA(2),";",nA(3),"];"
       write(20,*) "vecB = [",nB(1),";",nB(2),";",nB(3),"];"
       write(20,*) "M = vecB*pinv(vecA);"
       write(20,*) 'write("OUT",M);'
       write(20,*) 'quit'
       close(20)
 
       cmd="scilab-adv-cli -f scilab.sci -nb"
       call system(cmd)
 
       open(unit=21,file='OUT',form='formatted',status='unknown')
       read(21,*) M11,M12,M13
       read(21,*) M21,M22,M23
       read(21,*) M31,M32,M33
       close(21)
 
       q(4)=0.5*sqrt(1+M11+M22+M33)
       q(1)=0.25*(1./q(4))*(M32-M23)
       q(2)=0.25*(1./q(4))*(M13-M31)
       q(3)=0.25*(1./q(4))*(M21-M12)
 
       pi=2.*acos(0.d0)
 
       psi=(180./pi)*atan2((q(1)*q(3) + q(2)*q(4)),
      &-(q(2)*q(3) - q(1)*q(4)))
       theta =(180./pi)*asin(-2.*(q(2)*q(4)-q(3)*q(1)))
       phi=(180./pi)*atan2((q(1)*q(3) - q(2)*q(4)),
      &(q(2)*q(3) + q(1)*q(4))) 
 
       write(6,*) "psi   (deg) --",psi
       write(6,*) "theta (deg) --",theta
       write(6,*) "phi   (deg) --",phi
       END PROGRAM QUAT
Cepdendant quand j'entre en input pour ce code une rotation d'un vecteur unitaire sur l'axe x vers l'axe y , c'est à dire : a = (1,0,0) b = (0,1,0) j'obtiens :

Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
 ./quat.f 
1. 0. 0. 0. 1. 0.
output :
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
 
 psi   (deg) --   180.00000    
 theta (deg) --  -0.0000000    
 phi   (deg) --   0.0000000
Alors que j'attends phi = 90 et theta = psi = 0.
Après cela j'ai un doute sur les formules proposées dans la page wiki, est ce que quelqu'un aurait une idée ?

Merci d'avance,