I do a lot of rotations, and usually use quaternions. Experts say
they do not suffer from "gimbal lock" which can occur with Euler
angles. I know this is not the form you asked about, but just to
throw this out there, here's what I use:
------------------------------------------------------------------
NB. Quaternion maths.
NB. Four-element quaternion is scalar followed by 3-vector
NB. Should work on (n 4) arrays of quaternions.
L=: +/"1 &. (*:"_) NB. Gives the lengths of vector arrays
dot=: [: +/"1 * NB. dot products of arrays of vectors
pm1=: 1&|."1 NB. 1st cyclic permutation
pm2=: 2&|."1 NB. 2nd cyclic permutation
pm12=: 4 : '(pm1 x)*pm2 y'
cross=: pm12 - pm12~ NB. cross products of vectors
NB. cross=: pm12~ - pm12 NB. for left-handed coordinate system
NB. Multiply quaternions: q1 QM q2 --> q3
QM=: 4 : 0
s1=. 0{"1 x
v1=. }."1 x
s2=. 0{"1 y
v2=. }."1 y
s=. (s1*s2) - v1 dot v2
v=. (s1*v2)+(s2*v1)+(v1 cross v2)
s ,&.|: v NB. catenate under transpose
)
NB. The conjuate of the quaternions: Qb q --> q_bar
Qb=: 3 : '(1 _1 _1 _1)*"1 y'
NB. Quaternion for rotation of points by theta about unit vector u
NB. (Rotates *points*; use -theta to rotate coordinate axes)
NB. In right handed coordinates, looking in the +z direction, a
NB. positive rotation (from x to y) is *clockwise* -- we must use
NB. -theta to rotate counter clockwise (or make the axis -z).
NB. Usage: u Qrot theta --> q then q Pr P --> P'
Qrot=: 4 : 0
u=. x
s=. cos -: y
f=. sin -: y
s ,&.|: (f*u) NB. catenate under transpose
)
NB. Rotate point using quaternion q
NB. Usage: q Pr (x y z) --> (x' y' z')
Pr=: 4 : 0
qb=. Qb q=. x
}."1 q QM (0 ,&.|: y) QM qb
)
norm=: ] % L NB. normalize to unit length
NB. Rotate points P (P[n,3]) about some vector u
NB. by angle theta: (0 1 0; 0.25p1) Qrtn P --> P'
NB. (theta in radians; u need not be a unit vector)
Qrtn=: 4 : 0
'u theta'=. x
((norm u) Qrot theta) & Pr"1 y
)
----------------------------------------------------------------
Example: Construct 8 unit vectors in random directions
]vects=. norm 5 3$_1+2*rand 15
0.831204 0.549015 _0.0876507
0.551173 0.818662 0.161251
0.77789 0.430158 _0.458096
0.843154 _0.534975 0.0537832
_0.803313 _0.415795 _0.426383
L vects
1 1 1 1 1
Lets rotate them by 45 degrees about an axis in direction (_1 1.5 2):
]rvects =. (_1 1.5 2;0.25p1) Qrtn vects
0.272258 0.790675 _0.548369
_0.0170958 0.971275 _0.237343
0.186041 0.528875 _0.828058
0.960505 _0.0145176 _0.277885
_0.490408 _0.868698 0.069747
L rvects
1 1 1 1 1
(they retain their unit lengths)
Rotate them back (-45 degrees):
(_1 1.5 2;_0.25p1) Qrtn rvects
0.831204 0.549015 _0.0876507
0.551173 0.818662 0.161251
0.77789 0.430158 _0.458096
0.843154 _0.534975 0.0537832
_0.803313 _0.415795 _0.426383
Cheers,
Patrick
On Sun, 10 Jul 2016, Raul Miller wrote:
I need a routine which makes rotation matrices for a given angle.
For example, let's say that my angle is 0.05p1, I would be wanting
these two rotation matrices:
3 3$(2 1 o./ (,-)0.05p1) 0 4 1 3}&,=i.3
0.987688 0.156434 0
_0.156434 0.987688 0
0 0 1
and
3 3$(2 1 o./ (,-)0.05p1) 4 8 5 7}&,=i.3
1 0 0
0 0.987688 0.156434
0 _0.156434 0.987688
Now, obviously, I could implement this directly:
rotmat=:dyad define
3 3$(2 1 o./ (,-)y) (0 4 1 3+4*x)}&,=i.3
)
But it seems to me that there ought to be a better way of doing this.
And by better I mean more concise and direct.
But I'm not seeing that.
The best alternative I have thought of so far is:
elcric=:dyad define"0
if. 1<|y do.
(_2+|y) o. x**y
else.
y
end.
)
rotmat=:dyad define
(1 1-x) |. y elcric 1 0 0,0 4 3,:0 _3 4
)
And while maybe that satisfies some concept of "prettier" it is not
more concise and it is not more direct. It is not even faster.
That said, what I really want to build would be this:
rotate=:monad define
+/ .*/ 0 1 0 rotmat"0 y+0 0 0
)
(I threw in the unnecessary +0 0 0 to emphasize that I expect three angles.)
In other words, I have three different euler angles and I want the
rotation matrix they represent. So maybe there is even something
already implemented for this purpose? (But when I search jsoftware.com
for mentions of euler angles, I don't see anything promising.)
Any ideas?
Thanks,
--
Raul
P.S. I imagine I could live with getting the rotation quaternion they
represent, instead of the rotation matrix. But I have not investigated
the details of that approach here.
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
