Rotazione attorno a un asse

[marcomasetti]

Per costruire oggetti con parti che seguano anche direzioni oblique è necessario fare uso di rotazioni.
I comandi predefiniti basati su ROT risultano agevoli da utilizzare soltanto per rotazioni attorno a un asse coordinato, Volendo ruotare un oggetto, o sue parti, attorno ad una asse qualsiasi, occorre definire una formula, da utilizzare eventualmente con il comando di trasformazione XFORM.
La formula veramente utile è quella in fondo, che non dimostro direttamente, ma credo che la cosa veramente importante sia comprendere il ragionamento che permette di definire nuovi sistemi di riferimento. Sapere introdurre ulteriori terne di riferimento è fondamentale per studiare non solo movimenti e rotazioni, ma anche per costruire superfici.                                               

In Meccanica razionale, per semplificare le formule e renderle generalizzate anche per spazi con dimensioni superiori a 3, si scrivono le coordinate utilizzando una unica lettera seguita da un indice, che rappresenta l’asse di riferimento:
x1, x2, x3   corrispondono pertanto a:   x, y, z
L’indice i applicato ad una variabile indica la sommatoria di questa variabile entro le dimensioni dello spazio dove si sta operando,
ad esempio nello spazio 3D:  fi    sta per:     f1 + f2 + f3
Tuttavia questo metodo molto sintetico si presta a fare confusione.
Vogliamo determinare la nuova posizione R in cui viene a trovarsi il punto P dopo la rotazione dell’angolo ro attorno alle retta r, con punto di applicazione A e versore K.
Preliminarmente occorre trovare la proiezione ortogonale Q del punto P su r, che si determina ponendo pari a 0 il prodotto scalare dei vettori QP e K, in modo che risultino perpendicolari:
(pi-qi)*ki=0    cioè:   (px-qx)*kx+(py-qy)*ky+(pz-qz)*kz =0 
Risulta:  AQ=q*K, dove q è la lunghezza del segmento AQ, essendo K un versore per definizione.
Inoltre:  Q=A+AQ   cioè  (qx,qy,qz)=(ax,ay,az)+q*(ky,ky,kz)  ovvero:   
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
Da cui si desume:  (px- ax-q*kx)*kx+(py- ay-q*ky)*ky+(pz- az-q*kz)*kz =0
q*(kx^2+ky^2+kz^2)=(px- ax)*kx+(py- ay)*ky+(pz- az)*kz
Risultando  (kx^2+ky^2+kz^2)=1, dato che K è versore, in conclusione si ha:
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
(valore che poteva essere ottenuto direttamente come proiezione di AP su k)

Per effettuare la rotazione occorre definire un nuovo sistema di riferimento con origine Q e versori K,I,J.
Il versore J è fornito da:   QP/mod(QP)
md= sqr( (px-qx)^2+(py-qy)^2(pz-qz)^2 )
Jx=(px-qx)/md
Jy=(py-qy)/md
Jz=(pz-qz)/md
Il terzo versore I si ricava dal prodotto vettoriale di J per K (regola della mano destra):
ix=jy*kz-ky*jz
iy=jz*kx-kz*jx
iz=jx*ky-kx*jy
A questo punto si effettua la rotazione di P verso R.  Rispetto il nuovo sistema R ha coordinate:
( md*sin(ro) , md*cos(ro), 0 )
OR=OQ+QR=OQ+I*md*sin(ro)+J*md*cos(ro)
rx=qx+ ix*md*sin(ro)+jx*md*cos(ro)
ry=qy+ iy*md*sin(ro)+jy*md*cos(ro)
rz=qz+ iz*md*sin(ro)+jz*md*cos(ro)
Per dedurre i termini di XFORM, le formule devono essere esplicitate in modo da avere una funzione lineare da spazio a spazio, ovvero da (px,py,pz) a (rx,ry,rz), i cui coefficienti servono per dedurre i termini di XFORM:
jx=(px-qx) /md
jy=(py-qy) /md
jz=(pz-qz) /md
ix=((py-qy) *kz-ky*(pz-qz))/md
iy=((pz-qz) *kx-kz*(px-qx)) /md
iz=((px-qx) *ky-kx*(py-qy)) /md
rx=qx+((py-qy)*kz-ky*(pz-qz))*sin(ro)+(px-qx)*cos(ro)
ry=qy+((pz-qz)*kx-kz*(px-qx))*sin(ro)+(py-qy)*cos(ro)
rz=qz+((px-qx)*ky-kx*(py-qy))*sin(ro)+(pz-qz)*cos(ro)
I termini qx,qy,qz contengono però ancora i parametri px,py,pz e vanno esplicitati per potere utilizzare la matrice di trasformazione.
Tutta la sequenza può essere comunque utilizzata per definire singoli vertici.
Sostituiamo  sin(ro)  con  –sin(ro) per avere una rotazione antioraria attorno a K:
!Rotazione di un angolo antiorario=ro attorno al versore K=(kx,ky,kz) del punto P=(px,py,pz) applicato in  A=(ax,ay,az).
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
rx=qx-((py-qy)*kz-ky*(pz-qz))*sin(ro)+(px-qx)*cos(ro)
ry=qy-((pz-qz)*kx-kz*(px-qx))*sin(ro)+(py-qy)*cos(ro)
rz=qz-((px-qx)*ky-kx*(py-qy))*sin(ro)+(pz-qz)*cos(ro)
Proviamo ad utilizzare un’altra strada:
la proiezione ortogonale di K su xy è definita dal vettore (kx,ky) di lunghezza sqr(kx^2+ ky^2), che ha come vettore ortogonale  (ky,-kx), corrispondente al versore 3D:   V= (ky,-kx,0)/sqr(kx^2+ ky^2)
Il versore ortogonale ad entrambi è contenuto in un piano verticale e può essere calcolato direttamente o con il prodotto vettoriale di K con V.  Lungo tale piano verticale per r, il versore K può essere scomposto lungo la orizzontale e la verticale in:  (sqr(kx^2+ ky^2),kz); lungo tali direzioni il versore ortogonale a questo si ricava pertanto da: (-kz, sqr(kx^2+ ky^2), quindi il corrispondente versore 3D è dato da:
W=(-kz*kx/sqr(kx^2+ ky^2),-kz*ky/sqr(kx^2+ ky^2), sqr(kx^2+ ky^2))
infatti kx,ky individuano la direzione in orizzontale di r, mentre deve risultare: sqr(wx^2+wy^2)= kz, inoltre wx^2+ wy^2+ wz^2= kz^2+kx^2+ ky^2=1
                                         
Dunque:
vx=ky/sqr(kx^2+ ky^2)
vy=-kx/sqr(kx^2+ ky^2)
vz=0
wx=-kz*kx/sqr(kx^2+ ky^2)
wy=-kz*ky/sqr(kx^2+ ky^2)
wz=sqr(kx^2+ ky^2)
Oppure:
V=K^(0,0,1)/ mod( K^(0,0,1) )  !prodotto vettoriale
vx=ky/ sqr(kx^2+ ky^2)
vy=-kx/ sqr(kx^2+ ky^2)
vz=0
W=V^K
wx=vy*kz =-kx*kz/ sqr(kx^2+ ky^2)
wy=-kz*vx=-kz*ky/ sqr(kx^2+ ky^2)
wz= vx*ky-kx*vy=( ky^2+ kx^2 )/ sqr(kx^2+ ky^2)= sqr(kx^2+ ky^2)
La proiezione ortogonale di QP su v è data da:
pv=(px-qx)*vx+(py-qy)*vy+(pz-qz)*vz=(px-qx)*vx+(py-qy)*vy
Analogamente la proiezione ortogonale di QP su w è data da:
pw=(px-qx)*wx+(py-qy)*wy+(pz-qz)*wz
Per ruotare P utilizziamo la formula ricavata dai numeri complessi:
(x+i*y)*(cos(ro)+i*sin(ro))=x*cos(ro)-y*sin(ro) + i*(y*cos(ro)+x*sin(ro))
xx= x*cos(ro)-y*sin(ro)
yy= y*cos(ro)+x*sin(ro))
che nel caso in esame diventa:
rv=  pv*cos(ro)-pw*sin(ro)
rw= pw*cos(ro)+pv*sin(ro))
R=Q+ V*rv+ W*rw
rx=qx+ vx*rv+wy*rw
ry=qy+ vy*rv+wy*rw
rz=qz+ vz*rv+wz*rw
Esplicitare questa formula è ancora più complesso, tuttavia anche in questo caso può essere utilizzata su singoli vertici. Pure qui nello script sottostante cambiamo segno a  sin(ro)  per avere una rotazione antioraria.
!Rotazione di un angolo antiorario=ro  rispetto K attorno al versore K=(kx,ky,kz) del punto P=(px,py,pz) applicato in  A=(ax,ay,az).
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
vx= ky/sqr(kx^2+ ky^2)
vy=-kx/sqr(kx^2+ ky^2)
wx=-kz*kx/sqr(kx^2+ ky^2)
wy=-kz*ky/sqr(kx^2+ ky^2)
wz=sqr(kx^2+ ky^2)
pv=(px-qx)*vx+(py-qy)*vy
pw=(px-qx)*wx+(py-qy)*wy+(pz-qz)*wz
rv= pv*cos(ro)+pw*sin(ro)
rw= pw*cos(ro)-pv*sin(ro)
rx=qx+ vx*rv+wx*rw
ry=qy+ vy*rv+wy*rw
rz=qz+ vz*rv+wz*rw

Per ricavare la formula seguente, ho utilizzato invece un’altra strada ancora: prima ho eseguito la rotazione lungo un asse per l’origine, in modo da annullare ax,ay,az, poi ho effettuato una traslazione, tenendo conto dei valori differenti dovuti alla rotazione. Anche in questo caso, però, sono stati necessari molti calcoli, che ho eseguito su uno scartafaccio.

a11=kx^2*(1-cos(ro))+cos(ro)
a21=kx*ky*(1-cos(ro))-kz*sin(ro)
a31=kx*kz*(1-cos(ro))+ky*sin(ro)
a41=ax-a11*ax-a21*ay-a31*az
a12=kx*ky*(1-cos(ro))+kz*sin(ro)
a22=ky^2*(1-cos(ro))+cos(ro)
a32=ky*kz*(1-cos(ro))-kx*sin(ro)
a42=ay-a12*ax-a22*ay-a32*az
a13=kx*kz*(1-cos(ro))-ky*sin(ro)
a23=ky*kz*(1-cos(ro))+kx*sin(ro)
a33=kz^2*(1-cos(ro))+cos(ro)
a43=az-a13*ax-a23*ay-a33*az

XFORM   a11,a21,a31,a41,
                a12,a22,a32,a42,
                a13,a23,a33,a43
La formula di trasformazione risulta:
rx= a11*px+a21*py+a31*pz+a41
ry= a12*px+a22*py+a32*pz+a42
rz= a13*px+a23*py+a33*pz+a43
Corrispondente al comando
XFORM  a11, a21, a31, a41,
               a12, a22, a32, a42,
               a13, a23, a33, a43

Come verifica, con lo script sottostante applicato a un cubo, si può vedere che i diversi metodi conducono ad uno stesso risultato:
!Definire i parametri:
ax=
ay=
az=       !punto per l’asse
vax=
vay=
vaz=     !vettore direzione asse
kx=vax/sqr(vax^2+vay^2+vaz^2)
ky=vay/sqr(vax^2+vay^2+vaz^2)
kz=vaz/sqr(vax^2+vay^2+vaz^2)   !versore direzione asse
ro=       !angolo di rotazione

a11=kx^2*(1-cos(ro))+cos(ro)
a21=kx*ky*(1-cos(ro))-kz*sin(ro)
a31=kx*kz*(1-cos(ro))+ky*sin(ro)
a41=ax-a11*ax-a21*ay-a31*az
a12=kx*ky*(1-cos(ro))+kz*sin(ro)
a22=ky^2*(1-cos(ro))+cos(ro)
a32=ky*kz*(1-cos(ro))-kx*sin(ro)
a42=ay-a12*ax-a22*ay-a32*az
a13=kx*kz*(1-cos(ro))-ky*sin(ro)
a23=ky*kz*(1-cos(ro))+kx*sin(ro)
a33=kz^2*(1-cos(ro))+cos(ro)
a43=az-a13*ax-a23*ay-a33*az
xform a11,a21,a31,a41,
     a12,a22,a32,a42,
     a13,a23,a33,a43
block 1,1,1
del 1

!!!    metodo 1


BASE
px=0
py=0
pz=0 !#1
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
rx=qx-((py-qy)*kz-ky*(pz-qz))*sin(ro)+(px-qx)*cos(ro)
ry=qy-((pz-qz)*kx-kz*(px-qx))*sin(ro)+(py-qy)*cos(ro)
rz=qz-((px-qx)*ky-kx*(py-qy))*sin(ro)+(pz-qz)*cos(ro)
VERT rx,ry,rz

px=1
py=0
pz=0 !#2
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
rx=qx-((py-qy)*kz-ky*(pz-qz))*sin(ro)+(px-qx)*cos(ro)
ry=qy-((pz-qz)*kx-kz*(px-qx))*sin(ro)+(py-qy)*cos(ro)
rz=qz-((px-qx)*ky-kx*(py-qy))*sin(ro)+(pz-qz)*cos(ro)
VERT rx,ry,rz

px=1
py=1
pz=0 !#3
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
rx=qx-((py-qy)*kz-ky*(pz-qz))*sin(ro)+(px-qx)*cos(ro)
ry=qy-((pz-qz)*kx-kz*(px-qx))*sin(ro)+(py-qy)*cos(ro)
rz=qz-((px-qx)*ky-kx*(py-qy))*sin(ro)+(pz-qz)*cos(ro)
VERT rx,ry,rz

px=0
py=1
pz=0 !#4
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
rx=qx-((py-qy)*kz-ky*(pz-qz))*sin(ro)+(px-qx)*cos(ro)
ry=qy-((pz-qz)*kx-kz*(px-qx))*sin(ro)+(py-qy)*cos(ro)
rz=qz-((px-qx)*ky-kx*(py-qy))*sin(ro)+(pz-qz)*cos(ro)
VERT rx,ry,rz

px=0
py=0
pz=1 !#5
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
rx=qx-((py-qy)*kz-ky*(pz-qz))*sin(ro)+(px-qx)*cos(ro)
ry=qy-((pz-qz)*kx-kz*(px-qx))*sin(ro)+(py-qy)*cos(ro)
rz=qz-((px-qx)*ky-kx*(py-qy))*sin(ro)+(pz-qz)*cos(ro)
VERT rx,ry,rz

px=1
py=0
pz=1 !#6
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
rx=qx-((py-qy)*kz-ky*(pz-qz))*sin(ro)+(px-qx)*cos(ro)
ry=qy-((pz-qz)*kx-kz*(px-qx))*sin(ro)+(py-qy)*cos(ro)
rz=qz-((px-qx)*ky-kx*(py-qy))*sin(ro)+(pz-qz)*cos(ro)
VERT rx,ry,rz

px=1
py=1
pz=1 !#7
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
rx=qx-((py-qy)*kz-ky*(pz-qz))*sin(ro)+(px-qx)*cos(ro)
ry=qy-((pz-qz)*kx-kz*(px-qx))*sin(ro)+(py-qy)*cos(ro)
rz=qz-((px-qx)*ky-kx*(py-qy))*sin(ro)+(pz-qz)*cos(ro)
VERT rx,ry,rz

px=0
py=1
pz=1 !#8
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
rx=qx-((py-qy)*kz-ky*(pz-qz))*sin(ro)+(px-qx)*cos(ro)
ry=qy-((pz-qz)*kx-kz*(px-qx))*sin(ro)+(py-qy)*cos(ro)
rz=qz-((px-qx)*ky-kx*(py-qy))*sin(ro)+(pz-qz)*cos(ro)
VERT rx,ry,rz

EDGE 1, 2, -1, -1, 0 !#1
EDGE 2, 3, -1, -1, 0 !#2
EDGE 3, 4, -1, -1, 0 !#3
EDGE 4, 1, -1, -1, 0 !#4
EDGE 5, 6, -1, -1, 0 !#5
EDGE 6, 7, -1, -1, 0 !#6
EDGE 7, 8, -1, -1, 0 !#7
EDGE 8, 5, -1, -1, 0 !#8
EDGE 1, 5, -1, -1, 0 !#9
EDGE 2, 6, -1, -1, 0 !#10
EDGE 3, 7, -1, -1, 0 !#11
EDGE 4, 8, -1, -1, 0 !#12
PGON 4, 0, -1, -1, -4, -3, -2 !#1
PGON 4, 0, -1, 5, 6, 7, 8 !#2
PGON 4, 0, -1, 1, 10, -5, -9 !#3
PGON 4, 0, -1, 2, 11, -6, -10 !#4
PGON 4, 0, -1, 3, 12, -7, -11 !#5
PGON 4, 0, -1, 4, 9, -8, -12 !#6
BODY -1

!!!    metodo 2

BASE
px=0
py=0
pz=0 !#1
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
vx= ky/sqr(kx^2+ ky^2)
vy=-kx/sqr(kx^2+ ky^2)
wx=-kz*kx/sqr(kx^2+ ky^2)
wy=-kz*ky/sqr(kx^2+ ky^2)
wz=sqr(kx^2+ ky^2)
pv=(px-qx)*vx+(py-qy)*vy
pw=(px-qx)*wx+(py-qy)*wy+(pz-qz)*wz
rv= pv*cos(ro)+pw*sin(ro)
rw= pw*cos(ro)-pv*sin(ro)
rx=qx+ vx*rv+wx*rw
ry=qy+ vy*rv+wy*rw
rz=qz+ vz*rv+wz*rw
VERT rx,ry,rz

px=1
py=0
pz=0 !#2
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
vx= ky/sqr(kx^2+ ky^2)
vy=-kx/sqr(kx^2+ ky^2)
wx=-kz*kx/sqr(kx^2+ ky^2)
wy=-kz*ky/sqr(kx^2+ ky^2)
wz=sqr(kx^2+ ky^2)
pv=(px-qx)*vx+(py-qy)*vy
pw=(px-qx)*wx+(py-qy)*wy+(pz-qz)*wz
rv= pv*cos(ro)+pw*sin(ro)
rw= pw*cos(ro)-pv*sin(ro)
rx=qx+ vx*rv+wx*rw
ry=qy+ vy*rv+wy*rw
rz=qz+ vz*rv+wz*rw
VERT rx,ry,rz

px=1
py=1
pz=0 !#3
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
vx= ky/sqr(kx^2+ ky^2)
vy=-kx/sqr(kx^2+ ky^2)
wx=-kz*kx/sqr(kx^2+ ky^2)
wy=-kz*ky/sqr(kx^2+ ky^2)
wz=sqr(kx^2+ ky^2)
pv=(px-qx)*vx+(py-qy)*vy
pw=(px-qx)*wx+(py-qy)*wy+(pz-qz)*wz
rv= pv*cos(ro)+pw*sin(ro)
rw= pw*cos(ro)-pv*sin(ro)
rx=qx+ vx*rv+wx*rw
ry=qy+ vy*rv+wy*rw
rz=qz+ vz*rv+wz*rw
VERT rx,ry,rz

px=0
py=1
pz=0 !#4
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
vx= ky/sqr(kx^2+ ky^2)
vy=-kx/sqr(kx^2+ ky^2)
wx=-kz*kx/sqr(kx^2+ ky^2)
wy=-kz*ky/sqr(kx^2+ ky^2)
wz=sqr(kx^2+ ky^2)
pv=(px-qx)*vx+(py-qy)*vy
pw=(px-qx)*wx+(py-qy)*wy+(pz-qz)*wz
rv= pv*cos(ro)+pw*sin(ro)
rw= pw*cos(ro)-pv*sin(ro)
rx=qx+ vx*rv+wx*rw
ry=qy+ vy*rv+wy*rw
rz=qz+ vz*rv+wz*rw
VERT rx,ry,rz

px=0
py=0
pz=1 !#5
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
vx= ky/sqr(kx^2+ ky^2)
vy=-kx/sqr(kx^2+ ky^2)
wx=-kz*kx/sqr(kx^2+ ky^2)
wy=-kz*ky/sqr(kx^2+ ky^2)
wz=sqr(kx^2+ ky^2)
pv=(px-qx)*vx+(py-qy)*vy
pw=(px-qx)*wx+(py-qy)*wy+(pz-qz)*wz
rv= pv*cos(ro)+pw*sin(ro)
rw= pw*cos(ro)-pv*sin(ro)
rx=qx+ vx*rv+wx*rw
ry=qy+ vy*rv+wy*rw
rz=qz+ vz*rv+wz*rw
VERT rx,ry,rz

px=1
py=0
pz=1 !#6
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
vx= ky/sqr(kx^2+ ky^2)
vy=-kx/sqr(kx^2+ ky^2)
wx=-kz*kx/sqr(kx^2+ ky^2)
wy=-kz*ky/sqr(kx^2+ ky^2)
wz=sqr(kx^2+ ky^2)
pv=(px-qx)*vx+(py-qy)*vy
pw=(px-qx)*wx+(py-qy)*wy+(pz-qz)*wz
rv= pv*cos(ro)+pw*sin(ro)
rw= pw*cos(ro)-pv*sin(ro)
rx=qx+ vx*rv+wx*rw
ry=qy+ vy*rv+wy*rw
rz=qz+ vz*rv+wz*rw
VERT rx,ry,rz

px=1
py=1
pz=1 !#7
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
vx= ky/sqr(kx^2+ ky^2)
vy=-kx/sqr(kx^2+ ky^2)
wx=-kz*kx/sqr(kx^2+ ky^2)
wy=-kz*ky/sqr(kx^2+ ky^2)
wz=sqr(kx^2+ ky^2)
pv=(px-qx)*vx+(py-qy)*vy
pw=(px-qx)*wx+(py-qy)*wy+(pz-qz)*wz
rv= pv*cos(ro)+pw*sin(ro)
rw= pw*cos(ro)-pv*sin(ro)
rx=qx+ vx*rv+wx*rw
ry=qy+ vy*rv+wy*rw
rz=qz+ vz*rv+wz*rw
VERT rx,ry,rz

px=0
py=1
pz=1 !#8
q= (px-ax)*kx+(py-ay)*ky+(pz-az)*kz
qx=ax+q*kx
qy=ay+q*ky
qz=az+q*kz
vx= ky/sqr(kx^2+ ky^2)
vy=-kx/sqr(kx^2+ ky^2)
wx=-kz*kx/sqr(kx^2+ ky^2)
wy=-kz*ky/sqr(kx^2+ ky^2)
wz=sqr(kx^2+ ky^2)
pv=(px-qx)*vx+(py-qy)*vy
pw=(px-qx)*wx+(py-qy)*wy+(pz-qz)*wz
rv= pv*cos(ro)+pw*sin(ro)
rw= pw*cos(ro)-pv*sin(ro)
rx=qx+ vx*rv+wx*rw
ry=qy+ vy*rv+wy*rw
rz=qz+ vz*rv+wz*rw
VERT rx,ry,rz

EDGE 1, 2, -1, -1, 0 !#1
EDGE 2, 3, -1, -1, 0 !#2
EDGE 3, 4, -1, -1, 0 !#3
EDGE 4, 1, -1, -1, 0 !#4
EDGE 5, 6, -1, -1, 0 !#5
EDGE 6, 7, -1, -1, 0 !#6
EDGE 7, 8, -1, -1, 0 !#7
EDGE 8, 5, -1, -1, 0 !#8
EDGE 1, 5, -1, -1, 0 !#9
EDGE 2, 6, -1, -1, 0 !#10
EDGE 3, 7, -1, -1, 0 !#11
EDGE 4, 8, -1, -1, 0 !#12
PGON 4, 0, -1, -1, -4, -3, -2 !#1
PGON 4, 0, -1, 5, 6, 7, 8 !#2
PGON 4, 0, -1, 1, 10, -5, -9 !#3
PGON 4, 0, -1, 2, 11, -6, -10 !#4
PGON 4, 0, -1, 3, 12, -7, -11 !#5
PGON 4, 0, -1, 4, 9, -8, -12 !#6
BODY -1

[ZioBob]

Premetto che, molto prima di arrivare a metà della trattazione le mie competenze non mi permettono di comprendere le argomentazioni, che comunque leggo sempre fino in fondo.
Mi chiedevo se hai preso in considerazione un altro comando del GDL, meno complesso di XFORM ma comunque poco usato. Si tratta della quarta variante di ROT, quello "con asse generico". Oltre ai classici e comodi ROTx a, ROTy a e ROTz a c'è anche ROT x, y, z, a.
In questo caso x,y,z non definiscono un vero e proprio versore, ma un asse orientato relativo alla posizione corrente dell'origine (che quindi va prima traslata opportunamente).

[marcomasetti]

Si, ma non mi sembra molto flessibile. Poi c'è sempre il problema che ad ogni rotazione viene trascinato anche il sistema di riferimento e ci si scervella nel caso di rototraslazioni. Nonostante si utilizzi uno script un po' lunghetto, il comando che usa direttamente XFORM è molto più semplice, una volta che si sia chiarito come inserire i 6 valori:
3 per un punto di passaggio dell'asse di rotazione,
3 per la definizione della direzione, tramite un vettore oppure il versore,
ricavabile da questo dividendo le coordinate per il modulo.
E' poi possibile esprimere questo versore di direzione tramite le coordinate sferiche, utilizzando l'azimut o longitudine e la latitudine:
vx=cos(latitudine)*cos(azimut)
vy=cos(latitudine)*sin(azimut)
vz=sin(latitudine)

Lo script indiretto, presentato nella prima parte, diventa utile per inserire valori già ruotati con il comando PUT

[marcomasetti]

In effetti si può utilizzare pure il comando ROT che permette di ruotare dell’angolo ro un oggetto attorno ad un asse per O definito dal vettore VA=(vax, vay, vaz) :
ROT  vax, vay, vaz, ro
Se però l’asse passasse per un punto (ax,ay,az), come nel nostro caso, le cose si complicano. Occorrerebbe, infatti, traslare l’oggetto in modo che venga a trovarsi rispetto l’asse per O nella stessa posizione in cui si trovava inizialmente rispetto l’asse per A.  Per cui occorre prima determinare l’intersezione x,y sul piano orizzontale xy della retta r per A con direzione VA:
x=ax+vax*t
y=ay+vay*t
0=az+vaz*t     t=-az/vaz
x= ax-az*vax/vaz
y= ay-az*vay/vaz
Per cui l’oggetto va traslato di  add –x,-y,0 per averlo nella mutua posizione rispetto la retta poi il tutto va traslato nel senso opposto per riportare il sistema nella posizione iniziale.  La procedura sotto riportata conduce pertanto al medesimo risultato di rotazione di ro attorno ad r per A con direzione VA:

add ax -az*vax/vaz, ay -az*vay/vaz,0
ROT  vax, vay, vaz, ro
add az*vax/vaz-ax, az*vay/vaz-ay,0
block 1,1,1
del 3

Mi sembra che questo metodo in effetti funzioni, ma come succede con l'uso del calcolatore più che la logica ho usato l'intuizione ed invece di dimostrare la cosa la ho verificata per questo caso particolare: un metodo che farebbe indignare Euclide e i classici.