FM - Here is an example on TIP 759 new ocean of possibilities.
Working with Quaternions would imply to write dozen of expr command enclosed in many list commands.
But, thanks to Arithmetic compilation, any switch body case can now be compiled as a multiline arithmetic script, and thanks to the list comma-parentheses operator, lists of calculation are made easy.
Hence, there is no need to multiply procs in a big namespace ensemble. One proc is enough. I named it : Quaternion.
Each subcommands are handled by this only proc, in an infix-manner, so that what happen can be immediately understood.
Explanation :
Each operator is then catch by a dedicated switch, in charge to return the result.
proc Quaternion {args} {
if {![llength $args]} {return {0.0 0.0 0.0 0.0}}
if {[llength $args] == 1} {
# parse from string and return a list
[( qs = qi = qj = qk = 0 )]
set qs_pattern {([-+]?(?:[0-9]+(?:\.[0-9]+)?|\.[0-9]+)(?:[eE][-+]?[0-9]+)?)(?![ijk])}
regexp $qs_pattern $args -> qs
set qi_pattern {([-+]?(?:[0-9]*(?:\.[0-9]+)?|\.[0-9]+)(?:[eE][-+]?[0-9]+)?)(i)}
regexp $qi_pattern $args -> qi
if {$qi eq "+"} (qi=1) elseif {$qi eq "-"} (qi=-1)
set qj_pattern {([-+]?(?:[0-9]*(?:\.[0-9]+)?|\.[0-9]+)(?:[eE][-+]?[0-9]+)?)(j)}
regexp $qj_pattern $args -> qj
if {$qj eq "+"} (qj=1) elseif {$qj eq "-"} (qj=-1)
set qk_pattern {([-+]?(?:[0-9]*(?:\.[0-9]+)?|\.[0-9]+)(?:[eE][-+]?[0-9]+)?)(k)}
regexp $qk_pattern $args -> qk
if {$qk eq "+"} (qk=1) elseif {$qk eq "-"} (qk=-1)
return [($qs, $qi, $qj, $qk)]
}
if {[llength $args] == 2} {
# unary operators
lassign $args op Q
lassign $Q qs qi qj qk
return [switch -- $op {
"-" {(-$qs, -$qi, -$qj, - $qk )}
conj - conjugate {($qs, -$qi, -$qj, - $qk )}
mod - module - norm {(sqrt($qs**2+$qi**2 +$qj**2 + $qk**2))}
inverse {(
{n^2}=$qs**2+$qi**2 +$qj**2 + $qk**2;
double($qs) / ${n^2},
double($qi) / ${n^2},
double($qj) / ${n^2},
double($qk) / ${n^2}
)}
unit {(
n=sqrt($qs**2 + $qi**2 + $qj**2 + $qk**2);
double($qs) / $n,
double($qi) / $n,
double($qj) / $n,
double($qk) / $n
)}
radian {(
n = sqrt($qs**2 + $qi**2 + $qj**2 + $qk**2);
$n == 0 ? [error "norm of quaternion must not be null"]
: 2 * acos( double($qs)/$n )
)}
degree {(
pi = atan(1)*4;
n=sqrt($qs**2 + $qi**2 + $qj**2 + $qk**2);
$n == 0 ? [error "norm of quaternion must not be null"]
: 2 * acos( double($a)/$n )*180/$pi
)}
toMatrix {(
($qs, -$qi, -$qj, -$qk),
($qi, $qs, -$qk, $qj),
($qj, $qk, $qs, -$qi),
($qk, -$qj, $qi, $qs)
)}
format {("$qs[($qi<0?"-":"+")][(abs($qi))]i[($qj<0?"-":"+")][(abs($qj)]j[($qk<0?"-":"+")][(abs($qk)]k")}
fromVect {(0, $qs, $qi, $qj)}
}]
}
if {[llength $args] == 3} {
# binary operator
lassign $args P op Q
lassign $P ps pi pj pk
lassign $Q qs qi qj qk
return [switch -- $op {
"+" {($ps+$qs, $pi+$qi, $pj+$qj, $pk+$qk)}
"-" {($ps-$qs, $pi-$qi, $pj-$qj, $pk-$qk)}
"*" {(
$ps*$qs - $pi*$qi - $pj*$qj - $pk*$qk,
$ps*$qi + $pi*$qs + $pj*$qk - $pk*$qj,
$ps*$qj + $pj*$qs + $pk*$qi - $pi*$qk,
$ps*$qk + $pk*$qs + $pi*$qj - $pj*$qi
)}
x {(
# multiplication par un scalaire
$pi eq {} ? (
s = $ps;
$qs*$s, $qi*$s, $qj*$s, $qk*$s
) : $qi eq {} ? (
s = $qs;
$ps*$s, $pi*$s, $pj*$s, $pk*$s
) : [error "unable to assign scalar value"]
)}
rotor {(
# create a rotor from a vector and an angle
angle = $qs;
{angle/2}=double($angle)/90.0 * atan(1) ;
# norme du vecteur :
n = sqrt($pi**2 + $pj**2 + $pk**2);
# Retourner le rotor-Quaternion
cos(${angle/2}),
double($pi)/$n * sin(${angle/2}),
double($pj)/$n * sin(${angle/2}),
double($pk)/$n * sin(${angle/2})
)}
}]
}
}
# Applications : ref https://mecaspa.cannes-aero-patrimoine.net/SCAO/QUATERN/complements/som_quat.htm
proc Rotate {vAxe angle V} {(
qAxe = [Quaternion fromVect $vAxe]
U = [Quaternion $qAxe rotor $angle];
Qv = [Quaternion fromVect $V];
W = [Quaternion $U * [Quaternion $Qv * [Quaternion conj $U]]];
[lrange $W 1 end]
)}
# 90° Rotation of {0 2 0} around {1 0 0}
Rotate {1 0 0} 90 {0 2 0}
# 0.0 4.440892098500626e-16 2.0; (y coord has been send to z axis)
proc ProjOnPlane {vPlan vector} {(
U = [Quaternion unit [Quaternion fromVect $vPlan]];
V = [Quaternion fromVect $vector];
"U*" = [Quaternion conj $U];
VU = [Quaternion $V * $U];
UV = [Quaternion $U * $V];
W = [Quaternion \
[Quaternion [Quaternion $VU - $UV] * ${U*}]\
x 0.5];
[lrange $W 1 end]
)}
# Project the vector {1 0 0} on the plane directed by vector {1 1 1}
ProjOnPlane {1 1 1} {1 0 0}
# 0.6666666666666669 -0.3333333333333334 -0.3333333333333334
ProjOnPlane {1 1 1} {-1 -1 -1}
# 0.0 0.0 0.0
proc ProjOnAxe {vAxe vector} {(
U = [Quaternion unit [Quaternion fromVect $vAxe]];
V = [Quaternion fromVect $vector];
"U*" = [Quaternion conj $U];
VU = [Quaternion $V * $U];
UV = [Quaternion $U * $V];
W = [Quaternion \
[Quaternion [Quaternion $VU + $UV] * ${U*}]\
x 0.5];
[lrange $W 1 end]
)}
ProjOnAxe {1 1 1} {-1 -1 -1}
# -1.0000000000000002 -1.0000000000000002 -1.0000000000000002
proc SymByPlane {vPlan vector} {(
U = [Quaternion unit [Quaternion fromVect $vPlan]];
V = [Quaternion fromVect $vector];
"U*" = [Quaternion conj $U];
W = [Quaternion - [Quaternion $U * [Quaternion $V * ${U*}]]];
[lrange $W 1 end]
)}
SymByPlane {1 1 1} {-1 -1 -1}
# 1.0000000000000002 1.0000000000000002 1.0000000000000002