RRRR dual synthesis

1. Synthesis

Counting: nR not limited, nmax = 9 complete positions.

1. Definition of the RRRR chain in dual quaternions

quatmult[r_, 
        t_] := {r[[4]]*t[[1]] + t[[4]]*r[[1]] + 
          r[[2]]*t[[3]] - t[[2]]*r[[3]], 
        r[[4]]*t[[2]] + t[[4]]*r[[2]] + 
          r[[3]]*t[[1]] - t[[3]]*r[[1]], 
        r[[4]]*t[[3]] + t[[4]]*r[[3]] + 
          r[[1]]*t[[2]] - 
          t[[1]]*r[[2]], -r[[1]]*t[[1]] - 
          t[[2]]*r[[2]] - r[[3]]*t[[3]] + 
          t[[4]]*r[[4]]}

G = {Sin[θ/2] gx, Sin[θ/2] gy, Sin[θ/2] gz, 
          Cos[θ/2]} + ϵ {Sin[θ/2] gx0, 
            Sin[θ/2] gy0, Sin[θ/2] gz0, 0}

W = {Sin[ϕ/2] wx, Sin[ϕ/2] wy, Sin[ϕ/2] wz, 
          Cos[ϕ/2]} + ϵ {Sin[ϕ/2] wx0, Sin[ϕ/2] wy0, 
            Sin[ϕ/2] wz0, 0}

F = {Sin[ψ/2] fx, Sin[ψ/2] fy, Sin[ψ/2] fz, 
          Cos[ψ/2]} + ϵ {Sin[ψ/2] fx0, Sin[ψ/2] fy0, 
            Sin[ψ/2] fz0, 0}

U = {Sin[ξ/2] ux, Sin[ξ/2] uy, Sin[ξ/2] uz, 
          Cos[ξ/2]} + ϵ {Sin[ξ/2] ux0, Sin[ξ/2] uy0, 
            Sin[ξ/2] uz0, 0}

Qrr1 = 
      Simplify[Expand[quatmult[G, W]] /. ϵ^2 -> 0]

Qrr2 = 
      Simplify[Expand[quatmult[F, U]] /. ϵ^2 -> 0]

So that we don't have to repeat the following step, which takes forever, we just write the result below:

This is the result:

Q4r3 = {Cos[(θ/2)] (Cos[(ϕ/2)] (ux Cos[(ψ/2)] Sin[(ξ/2)] + (fx Cos[(ξ/2)] + (-fz uy + 
                                  fy uz) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] (Cos[(ψ/2)] (wx Cos[(ξ/2)] + (uz wy - uy wz) Sin[(ξ/2)]) + ((fz wy - 
                                  fy wz) Cos[(ξ/2)] + (-fx ux wx - fy uy wx - fz uz wx - fy ux wy + fx uy wy - fz ux wz + 
                                  fx uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
          Sin[(θ/2)] (Cos[(ϕ/2)] (Cos[(ψ/2)] (gx Cos[(ξ/2)] + (-gz uy + gy uz) Sin[(ξ/2)]) + ((fz gy - 
                                  fy gz) Cos[(ξ/2)] + (-fx gx ux - fy gy ux - fz gz ux - fy gx uy + fx gy uy - fz gx uz + 
                                  fx gz uz) Sin[(ξ/2)]) Sin[(ψ/2)]) + Sin[(ϕ/2)] (Cos[(ψ/2)] ((-gz wy + 
                                  gy wz) Cos[(ξ/2)] + (-gx ux wx + gy uy wx + gz uz wx - gy ux wy - gx uy wy - gz ux wz - 
                                  gx uz wz) Sin[(ξ/2)]) + ((-fx gx wx + fy gy wx + fz gz wx - fy gx wy - fx gy wy - 
                                  fz gx wz - 
                                  fx gz wz) Cos[(ξ/2)] + (fz gy ux wx - fy gz ux wx + fz gx uy wx + fx gz uy wx - fy gx uz wx - 
                                  fx gy uz wx - fz gx ux wy + 
                                  fx gz ux wy + fz gy uy wy + 
                                  fy gz uy wy + fx gx uz wy - 
                                  fy gy uz wy + fz gz uz wy + 
                                  fy gx ux wz - fx gy ux wz - 
                                  fx gx uy wz - fy gy uy wz + 
                                  fz gz uy wz - fz gy uz wz - 
                                  fy gz uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])) + ϵ (Cos[(θ/2)] (Cos[(ϕ/2)] (ux0 Cos[(ψ/2)] Sin[(ξ/2)] + (fx0 Cos[(ξ/2)] + (-fz0 uy - 
                                        fz uy0 + fy0 uz + 
                                        fy uz0) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                      Sin[(ϕ/2)] (Cos[(ψ/2)] (wx0 Cos[(ξ/2)] + (uz0 wy + uz wy0 - uy0 wz - 
                                        uy wz0) Sin[(ξ/2)]) + ((fz0 wy + fz wy0 - fy0 wz - 
                                        fy wz0) Cos[(ξ/2)] + (-fx0 ux wx - fx ux0 wx - fy0 uy wx - fy uy0 wx - fz0 uz wx - 
                                        fz uz0 wx - fx ux wx0 - 
                                        fy uy wx0 - fz uz wx0 - 
                                        fy0 ux wy - fy ux0 wy + 
                                        fx0 uy wy + fx uy0 wy - 
                                        fy ux wy0 + fx uy wy0 - 
                                        fz0 ux wz - fz ux0 wz + 
                                        fx0 uz wz + fx uz0 wz - 
                                        fz ux wz0 + 
                                        fx uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
                Sin[(θ/2)] (Cos[(ϕ/2)] (Cos[(ψ/2)] (gx0 Cos[(ξ/2)] + (-gz0 uy - gz uy0 + gy0 uz + 
                                        gy uz0) Sin[(ξ/2)]) + ((fz0 gy + fz gy0 - fy0 gz - 
                                        fy gz0) Cos[(ξ/2)] + (-fx0 gx ux - fx gx0 ux - fy0 gy ux - fy gy0 ux - fz0 gz ux - 
                                        fz gz0 ux - fx gx ux0 - 
                                        fy gy ux0 - fz gz ux0 - 
                                        fy0 gx uy - fy gx0 uy + 
                                        fx0 gy uy + fx gy0 uy - 
                                        fy gx uy0 + fx gy uy0 - 
                                        fz0 gx uz - fz gx0 uz + 
                                        fx0 gz uz + fx gz0 uz - 
                                        fz gx uz0 + 
                                        fx gz uz0) Sin[(ξ/2)]) Sin[(ψ/2)]) + Sin[(ϕ/2)] (Cos[(ψ/2)] ((-gz0 wy - 
                                        gz wy0 + gy0 wz + 
                                        gy wz0) Cos[(ξ/2)] + (-gx0 ux wx - gx ux0 wx + gy0 uy wx + gy uy0 wx + gz0 uz wx + 
                                        gz uz0 wx - gx ux wx0 + 
                                        gy uy wx0 + gz uz wx0 - 
                                        gy0 ux wy - gy ux0 wy - 
                                        gx0 uy wy - gx uy0 wy - 
                                        gy ux wy0 - gx uy wy0 - 
                                        gz0 ux wz - gz ux0 wz - 
                                        gx0 uz wz - gx uz0 wz - 
                                        gz ux wz0 - 
                                        gx uz wz0) Sin[(ξ/2)]) + ((-fx0 gx wx - fx gx0 wx + fy0 gy wx + fy gy0 wx + 
                                        fz0 gz wx + fz gz0 wx - 
                                        fx gx wx0 + fy gy wx0 + 
                                        fz gz wx0 - fy0 gx wy - 
                                        fy gx0 wy - fx0 gy wy - 
                                        fx gy0 wy - fy gx wy0 - 
                                        fx gy wy0 - fz0 gx wz - 
                                        fz gx0 wz - fx0 gz wz - 
                                        fx gz0 wz - fz gx wz0 - 
                                        fx gz wz0) Cos[(ξ/2)] + (fz0 gy ux wx + fz gy0 ux wx - fy0 gz ux wx - fy gz0 ux wx + 
                                        fz gy ux0 wx - fy gz ux0 wx + 
                                        fz0 gx uy wx + fz gx0 uy wx + 
                                        fx0 gz uy wx + fx gz0 uy wx + 
                                        fz gx uy0 wx + fx gz uy0 wx - 
                                        fy0 gx uz wx - fy gx0 uz wx - 
                                        fx0 gy uz wx - fx gy0 uz wx - 
                                        fy gx uz0 wx - fx gy uz0 wx + 
                                        fz gy ux wx0 - fy gz ux wx0 + 
                                        fz gx uy wx0 + fx gz uy wx0 - 
                                        fy gx uz wx0 - fx gy uz wx0 - 
                                        fz0 gx ux wy - fz gx0 ux wy + 
                                        fx0 gz ux wy + fx gz0 ux wy - 
                                        fz gx ux0 wy + fx gz ux0 wy + 
                                        fz0 gy uy wy + fz gy0 uy wy + 
                                        fy0 gz uy wy + fy gz0 uy wy + 
                                        fz gy uy0 wy + fy gz uy0 wy + 
                                        fx0 gx uz wy + fx gx0 uz wy - 
                                        fy0 gy uz wy - fy gy0 uz wy + 
                                        fz0 gz uz wy + fz gz0 uz wy + 
                                        fx gx uz0 wy - fy gy uz0 wy + 
                                        fz gz uz0 wy - fz gx ux wy0 + 
                                        fx gz ux wy0 + fz gy uy wy0 + 
                                        fy gz uy wy0 + fx gx uz wy0 - 
                                        fy gy uz wy0 + fz gz uz wy0 + 
                                        fy0 gx ux wz + fy gx0 ux wz - 
                                        fx0 gy ux wz - fx gy0 ux wz + 
                                        fy gx ux0 wz - fx gy ux0 wz - 
                                        fx0 gx uy wz - fx gx0 uy wz - 
                                        fy0 gy uy wz - fy gy0 uy wz + 
                                        fz0 gz uy wz + fz gz0 uy wz - 
                                        fx gx uy0 wz - fy gy uy0 wz + 
                                        fz gz uy0 wz - fz0 gy uz wz - 
                                        fz gy0 uz wz - fy0 gz uz wz - 
                                        fy gz0 uz wz - fz gy uz0 wz - 
                                        fy gz uz0 wz + fy gx ux wz0 - 
                                        fx gy ux wz0 - fx gx uy wz0 - 
                                        fy gy uy wz0 + fz gz uy wz0 - 
                                        fz gy uz wz0 - 
                                        fy gz uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)]))), 
        Cos[(θ/2)] (Cos[(ϕ/2)] (uy Cos[(ψ/2)] Sin[(ξ/2)] + (fy Cos[(ξ/2)] + (fz ux - 
                                  fx uz) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] (Cos[(ψ/2)] (wy Cos[(ξ/2)] + (-uz wx + ux wz) Sin[(ξ/2)]) + ((-fz wx + 
                                  fx wz) Cos[(ξ/2)] + (fy ux wx - 
                                  fx uy wx - fx ux wy - fy uy wy - 
                                  fz uz wy - fz uy wz + 
                                  fy uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
          Sin[(θ/2)] (Cos[(ϕ/2)] (Cos[(ψ/2)] (gy Cos[(ξ/2)] + (gz ux - gx uz) Sin[(ξ/2)]) + ((-fz gx + 
                                  fx gz) Cos[(ξ/2)] + (fy gx ux - 
                                  fx gy ux - fx gx uy - fy gy uy - 
                                  fz gz uy - fz gy uz + 
                                  fy gz uz) Sin[(ξ/2)]) Sin[(ψ/2)]) + Sin[(ϕ/2)] (Cos[(ψ/2)] ((gz wx - 
                                  gx wz) Cos[(ξ/2)] + (-gy ux wx - gx uy wx + gx ux wy - gy uy wy + gz uz wy - gz uy wz - 
                                  gy uz wz) Sin[(ξ/2)]) + ((-fy gx wx - fx gy wx + fx gx wy - fy gy wy + fz gz wy - 
                                  fz gy wz - 
                                  fy gz wz) Cos[(ξ/2)] + (-fz gx ux wx - fx gz ux wx + fz gy uy wx - fy gz uy wx + 
                                  fx gx uz wx - fy gy uz wx - 
                                  fz gz uz wx - fz gy ux wy - 
                                  fy gz ux wy - fz gx uy wy + 
                                  fx gz uy wy + fy gx uz wy + 
                                  fx gy uz wy + fx gx ux wz + 
                                  fy gy ux wz - fz gz ux wz + 
                                  fy gx uy wz - fx gy uy wz + 
                                  fz gx uz wz + 
                                  fx gz uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])) + ϵ (Cos[(θ/2)] (Cos[(ϕ/2)] (uy0 Cos[(ψ/2)] Sin[(ξ/2)] + (fy0 Cos[(ξ/2)] + (fz0 ux + 
                                        fz ux0 - fx0 uz - 
                                        fx uz0) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                      Sin[(ϕ/2)] (Cos[(ψ/2)] (wy0 Cos[(ξ/2)] + (-uz0 wx - uz wx0 + ux0 wz + 
                                        ux wz0) Sin[(ξ/2)]) + ((-fz0 wx - fz wx0 + fx0 wz + 
                                        fx wz0) Cos[(ξ/2)] + (fy0 ux wx + fy ux0 wx - fx0 uy wx - fx uy0 wx + fy ux wx0 - 
                                        fx uy wx0 - fx0 ux wy - 
                                        fx ux0 wy - fy0 uy wy - 
                                        fy uy0 wy - fz0 uz wy - 
                                        fz uz0 wy - fx ux wy0 - 
                                        fy uy wy0 - fz uz wy0 - 
                                        fz0 uy wz - fz uy0 wz + 
                                        fy0 uz wz + fy uz0 wz - 
                                        fz uy wz0 + 
                                        fy uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
                Sin[(θ/2)] (Cos[(ϕ/2)] (Cos[(ψ/2)] (gy0 Cos[(ξ/2)] + (gz0 ux + gz ux0 - gx0 uz - 
                                        gx uz0) Sin[(ξ/2)]) + ((-fz0 gx - fz gx0 + fx0 gz + 
                                        fx gz0) Cos[(ξ/2)] + (fy0 gx ux + fy gx0 ux - fx0 gy ux - fx gy0 ux + fy gx ux0 - 
                                        fx gy ux0 - fx0 gx uy - 
                                        fx gx0 uy - fy0 gy uy - 
                                        fy gy0 uy - fz0 gz uy - 
                                        fz gz0 uy - fx gx uy0 - 
                                        fy gy uy0 - fz gz uy0 - 
                                        fz0 gy uz - fz gy0 uz + 
                                        fy0 gz uz + fy gz0 uz - 
                                        fz gy uz0 + 
                                        fy gz uz0) Sin[(ξ/2)]) Sin[(ψ/2)]) + Sin[(ϕ/2)] (Cos[(ψ/2)] ((gz0 wx + 
                                        gz wx0 - gx0 wz - 
                                        gx wz0) Cos[(ξ/2)] + (-gy0 ux wx - gy ux0 wx - gx0 uy wx - gx uy0 wx - gy ux wx0 - 
                                        gx uy wx0 + gx0 ux wy + 
                                        gx ux0 wy - gy0 uy wy - 
                                        gy uy0 wy + gz0 uz wy + 
                                        gz uz0 wy + gx ux wy0 - 
                                        gy uy wy0 + gz uz wy0 - 
                                        gz0 uy wz - gz uy0 wz - 
                                        gy0 uz wz - gy uz0 wz - 
                                        gz uy wz0 - 
                                        gy uz wz0) Sin[(ξ/2)]) + ((-fy0 gx wx - fy gx0 wx - fx0 gy wx - fx gy0 wx - 
                                        fy gx wx0 - fx gy wx0 + 
                                        fx0 gx wy + fx gx0 wy - 
                                        fy0 gy wy - fy gy0 wy + 
                                        fz0 gz wy + fz gz0 wy + 
                                        fx gx wy0 - fy gy wy0 + 
                                        fz gz wy0 - fz0 gy wz - 
                                        fz gy0 wz - fy0 gz wz - 
                                        fy gz0 wz - fz gy wz0 - 
                                        fy gz wz0) Cos[(ξ/2)] + (-fz0 gx ux wx - fz gx0 ux wx - fx0 gz ux wx - 
                                        fx gz0 ux wx - fz gx ux0 wx - 
                                        fx gz ux0 wx + fz0 gy uy wx + 
                                        fz gy0 uy wx - fy0 gz uy wx - 
                                        fy gz0 uy wx + fz gy uy0 wx - 
                                        fy gz uy0 wx + fx0 gx uz wx + 
                                        fx gx0 uz wx - fy0 gy uz wx - 
                                        fy gy0 uz wx - fz0 gz uz wx - 
                                        fz gz0 uz wx + fx gx uz0 wx - 
                                        fy gy uz0 wx - fz gz uz0 wx - 
                                        fz gx ux wx0 - fx gz ux wx0 + 
                                        fz gy uy wx0 - fy gz uy wx0 + 
                                        fx gx uz wx0 - fy gy uz wx0 - 
                                        fz gz uz wx0 - fz0 gy ux wy - 
                                        fz gy0 ux wy - fy0 gz ux wy - 
                                        fy gz0 ux wy - fz gy ux0 wy - 
                                        fy gz ux0 wy - fz0 gx uy wy - 
                                        fz gx0 uy wy + fx0 gz uy wy + 
                                        fx gz0 uy wy - fz gx uy0 wy + 
                                        fx gz uy0 wy + fy0 gx uz wy + 
                                        fy gx0 uz wy + fx0 gy uz wy + 
                                        fx gy0 uz wy + fy gx uz0 wy + 
                                        fx gy uz0 wy - fz gy ux wy0 - 
                                        fy gz ux wy0 - fz gx uy wy0 + 
                                        fx gz uy wy0 + fy gx uz wy0 + 
                                        fx gy uz wy0 + fx0 gx ux wz + 
                                        fx gx0 ux wz + fy0 gy ux wz + 
                                        fy gy0 ux wz - fz0 gz ux wz - 
                                        fz gz0 ux wz + fx gx ux0 wz + 
                                        fy gy ux0 wz - fz gz ux0 wz + 
                                        fy0 gx uy wz + fy gx0 uy wz - 
                                        fx0 gy uy wz - fx gy0 uy wz + 
                                        fy gx uy0 wz - fx gy uy0 wz + 
                                        fz0 gx uz wz + fz gx0 uz wz + 
                                        fx0 gz uz wz + fx gz0 uz wz + 
                                        fz gx uz0 wz + fx gz uz0 wz + 
                                        fx gx ux wz0 + fy gy ux wz0 - 
                                        fz gz ux wz0 + fy gx uy wz0 - 
                                        fx gy uy wz0 + fz gx uz wz0 + 
                                        fx gz uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)]))), 
        Cos[(θ/2)] (Cos[(ϕ/2)] (uz Cos[(ψ/2)] Sin[(ξ/2)] + (fz Cos[(ξ/2)] + (-fy ux + 
                                  fx uy) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] (Cos[(ψ/2)] (wz Cos[(ξ/2)] + (uy wx - ux wy) Sin[(ξ/2)]) + ((fy wx - 
                                  fx wy) Cos[(ξ/2)] + (fz ux wx - 
                                  fx uz wx + fz uy wy - fy uz wy - 
                                  fx ux wz - fy uy wz - 
                                  fz uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
          Sin[(θ/2)] (Cos[(ϕ/2)] (Cos[(ψ/2)] (gz Cos[(ξ/2)] + (-gy ux + gx uy) Sin[(ξ/2)]) + ((fy gx - 
                                  fx gy) Cos[(ξ/2)] + (fz gx ux - 
                                  fx gz ux + fz gy uy - fy gz uy - 
                                  fx gx uz - fy gy uz - 
                                  fz gz uz) Sin[(ξ/2)]) Sin[(ψ/2)]) + Sin[(ϕ/2)] (Cos[(ψ/2)] ((-gy wx + 
                                  gx wy) Cos[(ξ/2)] + (-gz ux wx - gx uz wx - gz uy wy - gy uz wy + gx ux wz + gy uy wz - 
                                  gz uz wz) Sin[(ξ/2)]) + ((-fz gx wx - fx gz wx - fz gy wy - fy gz wy + fx gx wz + 
                                  fy gy wz - 
                                  fz gz wz) Cos[(ξ/2)] + (fy gx ux wx + fx gy ux wx - fx gx uy wx + fy gy uy wx + fz gz uy wx + 
                                  fz gy uz wx - fy gz uz wx - 
                                  fx gx ux wy + fy gy ux wy - 
                                  fz gz ux wy - fy gx uy wy - 
                                  fx gy uy wy - fz gx uz wy + 
                                  fx gz uz wy + fz gy ux wz + 
                                  fy gz ux wz - fz gx uy wz - 
                                  fx gz uy wz + fy gx uz wz - 
                                  fx gy uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])) + ϵ (Cos[(θ/2)] (Cos[(ϕ/2)] (uz0 Cos[(ψ/2)] Sin[(ξ/2)] + (fz0 Cos[(ξ/2)] + (-fy0 ux - 
                                        fy ux0 + fx0 uy + 
                                        fx uy0) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                      Sin[(ϕ/2)] (Cos[(ψ/2)] (wz0 Cos[(ξ/2)] + (uy0 wx + uy wx0 - ux0 wy - 
                                        ux wy0) Sin[(ξ/2)]) + ((fy0 wx + fy wx0 - fx0 wy - 
                                        fx wy0) Cos[(ξ/2)] + (fz0 ux wx + fz ux0 wx - fx0 uz wx - fx uz0 wx + fz ux wx0 - 
                                        fx uz wx0 + fz0 uy wy + 
                                        fz uy0 wy - fy0 uz wy - 
                                        fy uz0 wy + fz uy wy0 - 
                                        fy uz wy0 - fx0 ux wz - 
                                        fx ux0 wz - fy0 uy wz - 
                                        fy uy0 wz - fz0 uz wz - 
                                        fz uz0 wz - fx ux wz0 - 
                                        fy uy wz0 - 
                                        fz uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
                Sin[(θ/2)] (Cos[(ϕ/2)] (Cos[(ψ/2)] (gz0 Cos[(ξ/2)] + (-gy0 ux - gy ux0 + gx0 uy + 
                                        gx uy0) Sin[(ξ/2)]) + ((fy0 gx + fy gx0 - fx0 gy - 
                                        fx gy0) Cos[(ξ/2)] + (fz0 gx ux + fz gx0 ux - fx0 gz ux - fx gz0 ux + fz gx ux0 - 
                                        fx gz ux0 + fz0 gy uy + 
                                        fz gy0 uy - fy0 gz uy - 
                                        fy gz0 uy + fz gy uy0 - 
                                        fy gz uy0 - fx0 gx uz - 
                                        fx gx0 uz - fy0 gy uz - 
                                        fy gy0 uz - fz0 gz uz - 
                                        fz gz0 uz - fx gx uz0 - 
                                        fy gy uz0 - 
                                        fz gz uz0) Sin[(ξ/2)]) Sin[(ψ/2)]) + Sin[(ϕ/2)] (Cos[(ψ/2)] ((-gy0 wx - 
                                        gy wx0 + gx0 wy + 
                                        gx wy0) Cos[(ξ/2)] + (-gz0 ux wx - gz ux0 wx - gx0 uz wx - gx uz0 wx - gz ux wx0 - 
                                        gx uz wx0 - gz0 uy wy - 
                                        gz uy0 wy - gy0 uz wy - 
                                        gy uz0 wy - gz uy wy0 - 
                                        gy uz wy0 + gx0 ux wz + 
                                        gx ux0 wz + gy0 uy wz + 
                                        gy uy0 wz - gz0 uz wz - 
                                        gz uz0 wz + gx ux wz0 + 
                                        gy uy wz0 - 
                                        
                                        gz uz wz0) Sin[(ξ/2)]) + ((-fz0 gx wx - fz gx0 wx - fx0 gz wx - fx gz0 wx - 
                                        fz gx wx0 - fx gz wx0 - 
                                        fz0 gy wy - fz gy0 wy - 
                                        fy0 gz wy - fy gz0 wy - 
                                        fz gy wy0 - fy gz wy0 + 
                                        fx0 gx wz + fx gx0 wz + 
                                        fy0 gy wz + fy gy0 wz - 
                                        fz0 gz wz - fz gz0 wz + 
                                        fx gx wz0 + fy gy wz0 - 
                                        fz gz wz0) Cos[(ξ/2)] + (fy0 gx ux wx + fy gx0 ux wx + fx0 gy ux wx + fx gy0 ux wx + 
                                        fy gx ux0 wx + fx gy ux0 wx - 
                                        fx0 gx uy wx - fx gx0 uy wx + 
                                        fy0 gy uy wx + fy gy0 uy wx + 
                                        fz0 gz uy wx + fz gz0 uy wx - 
                                        fx gx uy0 wx + fy gy uy0 wx + 
                                        fz gz uy0 wx + fz0 gy uz wx + 
                                        fz gy0 uz wx - fy0 gz uz wx - 
                                        fy gz0 uz wx + fz gy uz0 wx - 
                                        fy gz uz0 wx + fy gx ux wx0 + 
                                        fx gy ux wx0 - fx gx uy wx0 + 
                                        fy gy uy wx0 + fz gz uy wx0 + 
                                        fz gy uz wx0 - fy gz uz wx0 - 
                                        fx0 gx ux wy - fx gx0 ux wy + 
                                        fy0 gy ux wy + fy gy0 ux wy - 
                                        fz0 gz ux wy - fz gz0 ux wy - 
                                        fx gx ux0 wy + fy gy ux0 wy - 
                                        fz gz ux0 wy - fy0 gx uy wy - 
                                        fy gx0 uy wy - fx0 gy uy wy - 
                                        fx gy0 uy wy - fy gx uy0 wy - 
                                        fx gy uy0 wy - fz0 gx uz wy - 
                                        fz gx0 uz wy + fx0 gz uz wy + 
                                        fx gz0 uz wy - fz gx uz0 wy + 
                                        fx gz uz0 wy - fx gx ux wy0 + 
                                        fy gy ux wy0 - fz gz ux wy0 - 
                                        fy gx uy wy0 - fx gy uy wy0 - 
                                        fz gx uz wy0 + fx gz uz wy0 + 
                                        fz0 gy ux wz + fz gy0 ux wz + 
                                        fy0 gz ux wz + fy gz0 ux wz + 
                                        fz gy ux0 wz + fy gz ux0 wz - 
                                        fz0 gx uy wz - fz gx0 uy wz - 
                                        fx0 gz uy wz - fx gz0 uy wz - 
                                        fz gx uy0 wz - fx gz uy0 wz + 
                                        fy0 gx uz wz + fy gx0 uz wz - 
                                        fx0 gy uz wz - fx gy0 uz wz + 
                                        fy gx uz0 wz - fx gy uz0 wz + 
                                        fz gy ux wz0 + fy gz ux wz0 - 
                                        fz gx uy wz0 - fx gz uy wz0 + 
                                        fy gx uz wz0 - 
                                        fx gy uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)]))), 
        Cos[(θ/2)] (Cos[(ϕ/2)] (Cos[(ξ/2)] Cos[(ψ/2)] + (-fx ux - fy uy - fz uz) Sin[(ξ/2)] Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] ((-ux wx - uy wy - 
                            uz wz) Cos[(ψ/2)] Sin[(ξ/2)] + ((-fx wx - fy wy - fz wz) Cos[(ξ/2)] + (fz uy wx - 
                                  fy uz wx - fz ux wy + fx uz wy + 
                                  fy ux wz - 
                                  fx uy wz) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
          Sin[(θ/2)] (Cos[(ϕ/2)] ((-gx ux - gy uy - 
                            gz uz) Cos[(ψ/2)] Sin[(ξ/2)] + ((-fx gx - fy gy - fz gz) Cos[(ξ/2)] + (-fz gy ux + 
                                  fy gz ux + fz gx uy - fx gz uy - 
                                  fy gx uz + 
                                  fx gy uz) Sin[(ξ/2)]) Sin[(ψ/2)]) + Sin[(ϕ/2)] (Cos[(ψ/2)] ((-gx wx - gy wy - 
                                  gz wz) Cos[(ξ/2)] + (-gz uy wx + gy uz wx + gz ux wy - gx uz wy - gy ux wz + 
                                  gx uy wz) Sin[(ξ/2)]) + ((fz gy wx - fy gz wx - fz gx wy + fx gz wy + fy gx wz - 
                                  fx gy wz) Cos[(ξ/2)] + (fx gx ux wx - fy gy ux wx - fz gz ux wx + fy gx uy wx + fx gy uy wx + 
                                  fz gx uz wx + fx gz uz wx + 
                                  fy gx ux wy + fx gy ux wy - 
                                  fx gx uy wy + fy gy uy wy - 
                                  fz gz uy wy + fz gy uz wy + 
                                  fy gz uz wy + fz gx ux wz + 
                                  fx gz ux wz + fz gy uy wz + 
                                  fy gz uy wz - fx gx uz wz - 
                                  fy gy uz wz + 
                                  fz gz uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])) + ϵ (Cos[(θ/2)] ((-fx0 ux - 
                            fx ux0 - fy0 uy - fy uy0 - fz0 uz - 
                            fz uz0) Cos[(ϕ/2)] Sin[(ξ/2)] Sin[(ψ/2)] + Sin[(ϕ/2)] ((-ux0 wx - ux wx0 - uy0 wy - 
                                  uy wy0 - uz0 wz - 
                                  uz wz0) Cos[(ψ/2)] Sin[(ξ/2)] + ((-fx0 wx - fx wx0 - fy0 wy - fy wy0 - fz0 wz - 
                                        fz wz0) Cos[(ξ/2)] + (fz0 uy wx + fz uy0 wx - fy0 uz wx - fy uz0 wx + fz uy wx0 - 
                                        fy uz wx0 - fz0 ux wy - 
                                        fz ux0 wy + fx0 uz wy + 
                                        fx uz0 wy - fz ux wy0 + 
                                        fx uz wy0 + fy0 ux wz + 
                                        fy ux0 wz - fx0 uy wz - 
                                        fx uy0 wz + fy ux wz0 - 
                                        fx uy wz0) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
                Sin[(θ/2)] (Cos[(ϕ/2)] ((-gx0 ux - 
                                  gx ux0 - gy0 uy - gy uy0 - gz0 uz - 
                                  gz uz0) Cos[(ψ/2)] Sin[(ξ/2)] + ((-fx0 gx - fx gx0 - fy0 gy - fy gy0 - fz0 gz - 
                                        fz gz0) Cos[(ξ/2)] + (-fz0 gy ux - fz gy0 ux + fy0 gz ux + fy gz0 ux - fz gy ux0 + 
                                        fy gz ux0 + fz0 gx uy + 
                                        fz gx0 uy - fx0 gz uy - 
                                        fx gz0 uy + fz gx uy0 - 
                                        fx gz uy0 - fy0 gx uz - 
                                        fy gx0 uz + fx0 gy uz + 
                                        fx gy0 uz - fy gx uz0 + 
                                        fx gy uz0) Sin[(ξ/2)]) Sin[(ψ/2)]) + Sin[(ϕ/2)] (Cos[(ψ/2)] ((-gx0 wx - 
                                        gx wx0 - gy0 wy - gy wy0 - 
                                        gz0 wz - 
                                        gz wz0) Cos[(ξ/2)] + (-gz0 uy wx - gz uy0 wx + gy0 uz wx + gy uz0 wx - gz uy wx0 + 
                                        gy uz wx0 + gz0 ux wy + 
                                        gz ux0 wy - gx0 uz wy - 
                                        gx uz0 wy + gz ux wy0 - 
                                        gx uz wy0 - gy0 ux wz - 
                                        gy ux0 wz + gx0 uy wz + 
                                        gx uy0 wz - gy ux wz0 + 
                                        gx uy wz0) Sin[(ξ/2)]) + ((fz0 gy wx + fz gy0 wx - fy0 gz wx - fy gz0 wx + fz gy wx0 - 
                                        fy gz wx0 - fz0 gx wy - 
                                        fz gx0 wy + fx0 gz wy + 
                                        fx gz0 wy - fz gx wy0 + 
                                        fx gz wy0 + fy0 gx wz + 
                                        fy gx0 wz - fx0 gy wz - 
                                        fx gy0 wz + fy gx wz0 - 
                                        fx gy wz0) Cos[(ξ/2)] + (fx0 gx ux wx + fx gx0 ux wx - fy0 gy ux wx - fy gy0 ux wx - 
                                        fz0 gz ux wx - fz gz0 ux wx + 
                                        fx gx ux0 wx - fy gy ux0 wx - 
                                        fz gz ux0 wx + fy0 gx uy wx + 
                                        fy gx0 uy wx + fx0 gy uy wx + 
                                        fx gy0 uy wx + fy gx uy0 wx + 
                                        fx gy uy0 wx + fz0 gx uz wx + 
                                        fz gx0 uz wx + fx0 gz uz wx + 
                                        fx gz0 uz wx + fz gx uz0 wx + 
                                        fx gz uz0 wx + fx gx ux wx0 - 
                                        fy gy ux wx0 - fz gz ux wx0 + 
                                        fy gx uy wx0 + fx gy uy wx0 + 
                                        fz gx uz wx0 + fx gz uz wx0 + 
                                        fy0 gx ux wy + fy gx0 ux wy + 
                                        fx0 gy ux wy + fx gy0 ux wy + 
                                        fy gx ux0 wy + fx gy ux0 wy - 
                                        fx0 gx uy wy - fx gx0 uy wy + 
                                        fy0 gy uy wy + fy gy0 uy wy - 
                                        fz0 gz uy wy - fz gz0 uy wy - 
                                        fx gx uy0 wy + fy gy uy0 wy - 
                                        fz gz uy0 wy + fz0 gy uz wy + 
                                        fz gy0 uz wy + fy0 gz uz wy + 
                                        fy gz0 uz wy + fz gy uz0 wy + 
                                        fy gz uz0 wy + fy gx ux wy0 + 
                                        fx gy ux wy0 - fx gx uy wy0 + 
                                        fy gy uy wy0 - fz gz uy wy0 + 
                                        fz gy uz wy0 + fy gz uz wy0 + 
                                        fz0 gx ux wz + fz gx0 ux wz + 
                                        fx0 gz ux wz + fx gz0 ux wz + 
                                        fz gx ux0 wz + fx gz ux0 wz + 
                                        fz0 gy uy wz + fz gy0 uy wz + 
                                        fy0 gz uy wz + fy gz0 uy wz + 
                                        fz gy uy0 wz + fy gz uy0 wz - 
                                        fx0 gx uz wz - fx gx0 uz wz - 
                                        fy0 gy uz wz - fy gy0 uz wz + 
                                        fz0 gz uz wz + fz gz0 uz wz - 
                                        fx gx uz0 wz - fy gy uz0 wz + 
                                        fz gz uz0 wz + fz gx ux wz0 + 
                                        fx gz ux wz0 + fz gy uy wz0 + 
                                        fy gz uy wz0 - fx gx uz wz0 - 
                                        fy gy uz wz0 + 
                                        fz gz uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)])))}

2. Separating the 8 equations

Q4rsep = {Coefficient[Q4r3[[1]], ϵ, 0], 
        Coefficient[Q4r3[[2]], ϵ, 0], 
        Coefficient[Q4r3[[3]], ϵ, 0], 
        Coefficient[Q4r3[[4]], ϵ, 0], 
        Coefficient[Q4r3[[1]], ϵ, 1], 
        Coefficient[Q4r3[[2]], ϵ, 1], 
        Coefficient[Q4r3[[3]], ϵ, 1], 
        Coefficient[Q4r3[[4]], ϵ, 1]}

Q4rsep2 = 
      Collect[Expand[
          Q4rsep], {Cos[(θ/2)]^(_ : 1) Sin[(ϕ/2)]^(_ : 1) , 
          Sin[(θ/2)]^(_ : 1) Cos[(ϕ/2)]^(_ : 1) , 
          Sin[(θ/2)]^(_ : 1) Sin[(ϕ/2)]^(_ : 1) , 
          Cos[(θ/2)]^(_ : 1) Cos[(ϕ/2)]^(_ : 1) , 
          Sin[(ψ/2)],  Cos[(ψ/2)], Sin[(ξ/2)], Cos[(ξ/2)]}, 
        Simplify]

3a. The parameterized equations

deseq = (Q4rsep2 /. {θ -> v1, ϕ -> 
                  v2, ψ -> v3, ξ -> v4})[[{1, 2, 3, 5, 
              6, 7}]] - {px, py, pz, px0, py0, pz0}

extraeqs = {-1 + (gx^2) + (gy^2) + (gz^2), -1 + (wx^2) + (wy^2) + 
          (wz^2), -1 + (fx^2) + (fy^2) + (fz^2), -1 + (ux^2) + (uy^2) + 
          (uz^2), gx gx0 + gy gy0 + gz gz0, wx wx0 + wy wy0 + wz wz0, 
        fx fx0 + fy fy0 + fz fz0, ux ux0 + uy uy0 + uz uz0}

Create the data file

puf = OpenWrite["RRRR.dat", PageWidth -> 200]

This writes out an expression, then closes the stream.

"********The design equations*******" >>> "RRRR.dat"

"********equ1*******" >>> "RRRR.dat"
deseq[[1]] >>> "RRRR.dat"
"********equ2*******" >>> "RRRR.dat"
deseq[[2]] >>> "RRRR.dat"
"********equ3*******" >>> "RRRR.dat"
deseq[[3]] >>> "RRRR.dat"
"********equ4*******" >>> "RRRR.dat"
deseq[[4]] >>> "RRRR.dat"
"********equ5*******" >>> "RRRR.dat"
deseq[[5]] >>> "RRRR.dat"
"********equ6*******" >>> "RRRR.dat"
deseq[[6]] >>> "RRRR.dat"

"********The extra equations*******" >>> "RRRR.dat"

extraeqs[[1]] >>> "RRRR.dat"
extraeqs[[2]] >>> "RRRR.dat"
extraeqs[[3]] >>> "RRRR.dat"
extraeqs[[4]] >>> "RRRR.dat"
extraeqs[[5]] >>> "RRRR.dat"
extraeqs[[6]] >>> "RRRR.dat"
extraeqs[[7]] >>> "RRRR.dat"
extraeqs[[8]] >>> "RRRR.dat"

Close[puf]

3. Reduction: Solving for the rotations

first4R = {Q4rsep2[[1]], Q4rsep2[[2]], Q4rsep2[[3]], 
        Q4rsep2[[4]]}

col1 = 
      Coefficient[first4R, Cos[(θ/2)] Sin[(ϕ/2)]]

col2 = 
      Coefficient[first4R, 
        Sin[(θ/2)] Cos[(ϕ/2)]]
col3 = 
      Coefficient[first4R, 
        Sin[(θ/2)] Sin[(ϕ/2)]]
col4 = 
      Coefficient[first4R, Cos[(θ/2)] Cos[(ϕ/2)]]

firstmat = Transpose[{col1, col2, col3, col4}]

MatrixForm[firstmat]

invmat = Inverse[firstmat];

den = Simplify[Denominator[invmat]]

den2 = 
      den /. {(fx^2) -> 1 - (fy^2) - (fz^2), 
          ux^2 -> 1 - uy^2 - uz^2}

Notice this is -(g×w).(g×w)=-((gz wx-gx wz)^2+(gz wy-gy wz)^2+(gy wx-gx wy)^2)  , it is different from zero if G and W are not parallel!!

Simplify[
      den2 + ((gz wx - gx wz)^2 + (gz wy - 
                  gy wz)^2 + (gy wx - gx wy)^2)]

We consider a degenerate case when two of the rotation axes are parallel. Hence we can eliminate the denominator g×w.g×w

Simplify[
        Numerator[invmat[[1, 1]]]] /. {fx^2 -> 1 - fy^2 - fz^2, 
        ux^2 -> 1 - uy^2 - uz^2}

To avoid this calculation....

invmatN2 = {{-Cos[(ξ/2)] (((gy^2) wx - gx gy wy + 
                        gz (gz wx - 
                              gx wz)) Cos[(ψ/2)] + (fy (gx gz wx + gy gz wy - (gx^2) wz - (gy^2) wz) + 
                        fz (-gx gy wx + (gx^2) wy + (gz^2) wy - 
                              gy gz wz)) Sin[(ψ/2)]) - 
            Sin[(ξ/2)] ((gx (gz uy - gy uz) wx - (gy uy + 
                              gz uz) (-gz wy + gy wz) + 
                        (gx^2) (uz wy - 
                              uy wz)) Cos[(ψ/2)] + (fy (-(gy^2) uy wx - (gx^2) ux wy - (gz^2) (uy wx + ux wy) + 
                              gy gz ux wz + 
                              gx (gy ux wx + gy uy wy + 
                                    gz uy wz)) + 
                        fx ((gz^2) (-ux wx + uy wy) + 
                              gx gz (-uz wx + ux wz) - 
                              gy (gx uy wx - gx ux wy + gz uz wy + 
                                    gz uy wz) + 
                              (gy^2) (-ux wx + uz wz) + 
                              (gx^2) (uy wy + uz wz)) + 
                        fz (-(gz^2) uz wx + gy gz ux wy - 
                              (gx^2) ux wz - (gy^2) (uz wx + ux wz) + 
                              gx (gz ux wx + gy uz wy + 
                                    gz uz wz))) Sin[(ψ/2)]), -Cos[(ξ/2)] ((-gx gy wx + (gx^2) wy + 
                        gz (gz wy - 
                              gy wz)) Cos[(ψ/2)] + (fx (-gx gz wx - gy gz wy + (gx^2) wz + (gy^2) wz) + 
                        fz (-(gy^2) wx - (gz^2) wx + gx gy wy + 
                              gx gz wz)) Sin[(ψ/2)]) - 
            Sin[(ξ/2)] ((-(gy^2) uz wx - (gz^2) uz wx - 
                        gy gz ux wy + (gx^2) ux wz + (gy^2) ux wz + 
                        gx (-gz ux wx + gy uz wy + 
                              gz uz wz)) Cos[(ψ/2)] + (fz (gx (gz uy + gy uz) wx - (gy uy - gz uz) (-gz wy + 
                                    gy wz) - 
                              (gx^2) (uz wy + uy wz)) + 
                        fx (-(gy^2) uy wx - (gx^2) ux wy - 
                              (gz^2) (uy wx + ux wy) + gy gz ux wz + 
                              gx (gy ux wx + gy uy wy + 
                                    gz uy wz)) + 
                        fy ((gz^2) (ux wx - uy wy) - 
                              gx gz (uz wx + ux wz) + 
                              gy (gx uy wx - gx ux wy - gz uz wy + 
                                    gz uy wz) + 
                              (gy^2) (ux wx + uz wz) + 
                              (gx^2) (-uy wy + 
                                    uz wz))) Sin[(ψ/2)]), 
          Cos[(ξ/2)] ((gx gz wx - (gx^2) wz + 
                        gy (gz wy - 
                              gy wz)) Cos[(ψ/2)] + (fy (-(gy^2) wx - (gz^2) wx + gx gy wy + gx gz wz) + 
                        fx (-gx gy wx + (gx^2) wy + (gz^2) wy - 
                              gy gz wz)) Sin[(ψ/2)]) + 
            Sin[(ξ/2)] ((-(gy^2) uy wx - (gz^2) uy wx + 
                        (gx^2) ux wy + (gz^2) ux wy - gy gz ux wz + 
                        gx (-gy ux wx + gy uy wy + 
                              gz uy wz)) Cos[(ψ/2)] + (fy (-gx (gz uy + gy uz) wx + (gy uy - 
                                    gz uz) (-gz wy + gy wz) + 
                              (gx^2) (uz wy + uy wz)) + 
                        fz (-(gz^2) ux wx - gx gz uz wx - 
                              (gx^2) uy wy - (gz^2) uy wy + 
                              gx gz ux wz + (gx^2) uz wz + 
                              gy (gx uy wx + gx ux wy - gz uz wy + 
                                    gz uy wz) + 
                              (gy^2) (-ux wx + uz wz)) + 
                        fx ((gz^2) uz wx - gy gz ux wy + 
                              (gx^2) ux wz + (gy^2) (uz wx + ux wz) - 
                              gx (gz ux wx + gy uz wy + 
                                    gz uz wz))) Sin[(ψ/2)]), ((gz^2) (ux wx + uy wy) - gx gz (uz wx + ux wz) - 
                  gy (gx uy wx + gx ux wy + gz uz wy + 
                        gz uy wz) + (gy^2) (ux wx + uz wz) + 
                  (gx^2) (uy wy + 
                        uz wz)) Cos[(ψ/2)] Sin[(ξ/2)] + ((fz (-gx gz wx - gy gz wy + (gx^2) wz + (gy^2) wz) + 
                        fx ((gy^2) wx + (gz^2) wx - gx gy wy - 
                              gx gz wz) + 
                        fy (-gx gy wx + (gx^2) wy + 
                              gz (gz wy - 
                                    gy wz))) Cos[(ξ/2)] + (fx (gx (-gz uy wx + gy uz wx) + (gy uy + 
                                    gz uz) (-gz wy + gy wz) + 
                              (gx^2) (-uz wy + uy wz)) + 
                        fz (-(gy^2) uy wx - (gz^2) uy wx + 
                              (gx^2) ux wy + (gz^2) ux wy - 
                              gy gz ux wz + 
                              gx (-gy ux wx + gy uy wy + 
                                    gz uy wz)) + 
                        fy ((gz^2) uz wx + gy gz ux wy - 
                              (gx^2) ux wz + (gy^2) (uz wx - ux wz) + 
                              gx (gz ux wx - gy uz wy - 
                                    gz uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]}, {Cos[(ξ/2)] ((gy wx wy + gz wx wz - 
                        gx ((wy^2) + 
                              (wz^2))) Cos[(ψ/2)] + (fy (gz ((wx^2) + (wy^2)) - (gx wx + gy wy) wz) + 
                        fz (wy (gx wx + gz wz) - 
                              gy ((wx^2) + 
                                    (wz^2)))) Sin[(ψ/2)]) + 
            Sin[(ξ/2)] ((-gy uz (wx^2) + gx uz wx wy - 
                        gx uy wx wz - gy uy wy wz - gy uz (wz^2) + 
                        gz (uy ((wx^2) + (wy^2)) + 
                              uz wy wz)) Cos[(ψ/2)] + (fx (-gz (uz ((wx^2) + (wy^2)) + (ux wx - uy wy) wz) - 
                              gy (wy (ux wx - uz wz) + 
                                    uy ((wx^2) + (wz^2))) + 
                              gx (uy wx wy + uz wx wz + 
                                    ux ((wy^2) + (wz^2)))) + 
                        fy (-gz (uy wx + ux wy) wz + 
                              gy (-uy wx wy + 
                                    ux ((wx^2) + (wz^2))) + 
                              gx (-ux wx wy + 
                                    uy ((wy^2) + (wz^2)))) + 
                        fz (-gy wy (uz wx + ux wz) + 
                              gz (ux ((wx^2) + (wy^2)) - 
                                    uz wx wz) + 
                              gx (-ux wx wz + 
                                    uz ((wy^2) + 
                                        (wz^2))))) Sin[(ψ/2)]), 
          Cos[(ξ/2)] ((wy (gx wx + gz wz) - 
                        gy ((wx^2) + 
                              (wz^2))) Cos[(ψ/2)] + (fx (-gz ((wx^2) + (wy^2)) + (gx wx + gy wy) wz) + 
                        fz (-gy wx wy - gz wx wz + 
                              gx ((wy^2) + 
                                    (wz^2)))) Sin[(ψ/2)]) + 
            Sin[(ξ/2)] ((-gy uz wx wy + gx uz (wy^2) + 
                        gx ux wx wz + gy ux wy wz + gx uz (wz^2) - 
                        gz (ux ((wx^2) + (wy^2)) + 
                              uz wx wz)) Cos[(ψ/2)] + (fz (-gx wx (uz wy + uy wz) + 
                              gz (uy ((wx^2) + (wy^2)) - 
                                    uz wy wz) + 
                              gy (-uy wy wz + 
                                    uz ((wx^2) + (wz^2)))) + 
                        fy (-gz (uz ((wx^2) + (wy^2)) - 
                                    ux wx wz + uy wy wz) + 
                              gy (wy (ux wx + uz wz) + 
                                    uy ((wx^2) + (wz^2))) - 
                              gx (uy wx wy - uz wx wz + 
                                    ux ((wy^2) + (wz^2)))) + 
                        fx (-gz (uy wx + ux wy) wz + 
                              gy (-uy wx wy + 
                                    ux ((wx^2) + (wz^2))) + 
                              gx (-ux wx wy + 
                                    uy ((wy^2) + 
                                        (wz^2))))) Sin[(ψ/2)]), -Cos[(ξ/2)] ((gz ((wx^2) + (wy^2)) - (gx wx + 
                              gy wy) wz) Cos[(ψ/2)] + (fx (wy (gx wx + gz wz) - gy ((wx^2) + (wz^2))) + 
                        fy (-gy wx wy - gz wx wz + 
                              gx ((wy^2) + 
                                    (wz^2)))) Sin[(ψ/2)]) - 
            Sin[(ξ/2)] ((gz (-uy wx + ux wy) wz - 
                        gy (uy wx wy + ux ((wx^2) + (wz^2))) + 
                        gx (ux wx wy + 
                              uy ((wy^2) + 
                                    (wz^2)))) Cos[(ψ/2)] + (fy (gx wx (uz wy + uy wz) - 
                              gz (uy ((wx^2) + (wy^2)) - 
                                    uz wy wz) - 
                              gy (-uy wy wz + 
                                    uz ((wx^2) + (wz^2)))) + 
                        fz (-gz (uz ((wx^2) + 
                                        (wy^2)) + (ux wx + 
                                        uy wy) wz) + 
                              gy (-ux wx wy + uz wy wz + 
                                    uy ((wx^2) + (wz^2))) + 
                              gx (-uy wx wy + uz wx wz + 
                                    ux ((wy^2) + (wz^2)))) + 
                        fx (gy wy (uz wx + ux wz) - 
                              gz (ux ((wx^2) + (wy^2)) - 
                                    uz wx wz) - 
                              gx (-ux wx wz + 
                                    uz ((wy^2) + 
                                        (wz^2))))) Sin[(ψ/2)]), (gz (uz ((wx^2) + (wy^2)) - (ux wx + uy wy) wz) + 
                  gy (-wy (ux wx + uz wz) + 
                        uy ((wx^2) + (wz^2))) + 
                  gx (-uy wx wy - uz wx wz + 
                        ux ((wy^2) + 
                              (wz^2)))) Cos[(ψ/2)] Sin[(ξ/2)] + ((fz (gz ((wx^2) + (wy^2)) - (gx wx + gy wy) wz) + 
                        fy (-wy (gx wx + gz wz) + 
                              gy ((wx^2) + (wz^2))) + 
                        fx (-gy wx wy - gz wx wz + 
                              gx ((wy^2) + 
                                    (wz^2)))) Cos[(ξ/2)] + (fy (-gy uz wx wy + gx uz (wy^2) + gx ux wx wz + gy ux wy wz + 
                              gx uz (wz^2) - 
                              gz (ux ((wx^2) + (wy^2)) + 
                                    uz wx wz)) + 
                        fx (-gy uz (wx^2) + gx uz wx wy - 
                              gx uy wx wz - gy uy wy wz - 
                              gy uz (wz^2) + 
                              gz (uy ((wx^2) + (wy^2)) + 
                                    uz wy wz)) + 
                        fz (-gx ux wx wy - gx uy (wy^2) + 
                              gz uy wx wz - gz ux wy wz - 
                              gx uy (wz^2) + 
                              gy (uy wx wy + 
                                    ux ((wx^2) + 
                                        (wz^2))))) Sin[(ξ/2)]) Sin[(ψ/2)]}, {Cos[(ξ/2)] ((gz wy - 
                        gy wz) Cos[(ψ/2)] + (-fy gy wx - 
                        fz gz wx + fy gx wy + 
                        fz gx wz) Sin[(ψ/2)]) + 
            Sin[(ξ/2)] ((-gy uy wx - gz uz wx + gx uy wy + 
                        gx uz wz) Cos[(ψ/2)] + (fx (gy uz wx - gx uz wy - gz (uy wx + ux wy) + gy ux wz + gx uy wz) + 
                        fy (gz ux wx - gz uy wy - gx ux wz + 
                              gy uy wz) + 
                        fz (-gy ux wx + gx ux wy - gz uz wy + 
                              gy uz wz)) Sin[(ψ/2)]), 
          Cos[(ξ/2)] ((-gz wx + 
                        gx wz) Cos[(ψ/2)] + (fx gy wx - 
                        fx gx wy - fz gz wy + 
                        fz gy wz) Sin[(ψ/2)]) + 
            Sin[(ξ/2)] ((-(gx ux + gz uz) wy + 
                        gy (ux wx + 
                              uz wz)) Cos[(ψ/2)] + (fy (gz uy wx + gy uz wx + gz ux wy - gx uz wy - gy ux wz - gx uy wz) + 
                        fx (gz ux wx - gz uy wy - gx ux wz + 
                              gy uy wz) + 
                        fz (-gy uy wx + gz uz wx + gx uy wy - 
                              gx uz wz)) Sin[(ψ/2)]), 
          Cos[(ξ/2)] ((gy wx - 
                        gx wy) Cos[(ψ/2)] + (fx gz wx + 
                        fy gz wy - fx gx wz - 
                        fy gy wz) Sin[(ψ/2)]) + 
            Sin[(ξ/2)] ((gz (ux wx + uy wy) - (gx ux + 
                              gy uy) wz) Cos[(ψ/2)] + (fz (-gz uy wx - gy uz wx + gz ux wy + gx uz wy - gy ux wz + 
                              gx uy wz) + 
                        fy (-gy uy wx + gz uz wx + gx uy wy - 
                              gx uz wz) + 
                        fx (-gy ux wx + gx ux wy - gz uz wy + 
                              gy uz wz)) Sin[(ψ/2)]), (gz uy wx - gy uz wx - gz ux wy + gx uz wy + gy ux wz - 
                  gx uy wz) Cos[(ψ/2)] Sin[(ξ/2)] + ((-fz gy wx + fy gz wx + fz gx wy - fx gz wy - fy gx wz + 
                        fx gy wz) Cos[(ξ/2)] + (fz (gz ux wx + 
                              gz uy wy - gx ux wz - gy uy wz) + 
                        fx (-gy uy wx - gz uz wx + gx uy wy + 
                              gx uz wz) + 
                        fy (-(gx ux + gz uz) wy + 
                              gy (ux wx + 
                                    uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]}, {-(-(ux wx + uy wy + 
                            uz wz) Cos[(ψ/2)] Sin[(ξ/2)] + (-(fx wx + fy wy + fz wz) Cos[(ξ/2)] + (fz uy wx - 
                                fy uz wx - fz ux wy + fx uz wy + 
                                fy ux wz - 
                                fx uy wz) Sin[(ξ/2)]) Sin[(ψ/2)]) ((Cos[(ψ/2)] ((-gy wx + 
                                    gx wy) Cos[(ξ/2)] - (gx uz wx + gy uz wy - gx ux wz - gy uy wz + 
                                    gz (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + (-(fx gz wx + fy gz wy - fx gx wz - fy gy wz + 
                                      fz (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fz (gz uy wx + gy uz wx - gz ux wy - gx uz wy + gy ux wz - 
                                        gx uy wz) + 
                                    fx (gy ux wx - gx uy wx - 
                                        gx ux wy - gy uy wy + 
                                        gz uz wy - gz uy wz - 
                                        gy uz wz) + 
                                    fy (gy uy wx - gz uz wx + 
                                        gy ux wy + gz ux wz + 
                                        gx (ux wx - uy wy + 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (gy Cos[(ψ/2)] + (-fz gx + 
                                    fx gz) Sin[(ψ/2)]) + 
                        Sin[(ξ/2)] ((gz ux - 
                                    gx uz) Cos[(ψ/2)] + (fy gx ux - fx gy ux - fx gx uy - fy gy uy - fz gz uy - fz gy uz + 
                                    fy gz uz) Sin[(ψ/2)])) - (Cos[(ψ/2)] ((gz wx - gx wz) Cos[(ξ/2)] - (gx uy wx - 
                                    gx ux wy - gz uz wy + gz uy wz + 
                                    gy (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + (-(fx gy wx - fx gx wy - fz gz wy + fz gy wz + 
                                      fy (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fy (-gy uz wx + gx uz wy - gz (uy wx + ux wy) + gy ux wz + 
                                        gx uy wz) - 
                                    fz (-gy uy wx + gz uz wx + 
                                        gy ux wy + gz ux wz + 
                                        gx (ux wx + uy wy - 
                                        uz wz)) + 
                                    fx (gx uz wx + gy uz wy + 
                                        gx ux wz - gy uy wz + 
                                        gz (-ux wx + uy wy + 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (gz Cos[(ψ/2)] + (fy gx - 
                                    fx gy) Sin[(ψ/2)]) - 
                        Sin[(ξ/2)] ((gy ux - 
                                    gx uy) Cos[(ψ/2)] + (fx gz ux + fy gz uy + fx gx uz + fy gy uz - 
                                    fz (gx ux + gy uy - 
                                        gz uz)) Sin[(ψ/2)]))) + (-(gx ux + gy uy + 
                          gz uz) Cos[(ψ/2)] Sin[(ξ/2)] + (-(fx gx + fy gy + fz gz) Cos[(ξ/2)] + (-fz gy ux + 
                              fy gz ux + fz gx uy - fx gz uy - 
                              fy gx uz + 
                              fx gy uz) Sin[(ξ/2)]) Sin[(ψ/2)]) (-(Cos[(ψ/2)] ((gz wx - 
                                      gx wz) Cos[(ξ/2)] - (gx uy wx - gx ux wy - gz uz wy + gz uy wz + 
                                      gy (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + (-(fx gy wx - fx gx wy - fz gz wy + fz gy wz + 
                                        fy (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fy (-gy uz wx + gx uz wy - gz (uy wx + ux wy) + gy ux wz + 
                                        gx uy wz) - 
                                      fz (-gy uy wx + gz uz wx + 
                                        gy ux wy + gz ux wz + 
                                        gx (ux wx + uy wy - 
                                        uz wz)) + 
                                      fx (gx uz wx + gy uz wy + 
                                        gx ux wz - gy uy wz + 
                                        gz (-ux wx + uy wy + 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (wz Cos[(ψ/2)] + (fy wx - 
                                    fx wy) Sin[(ψ/2)]) + 
                        Sin[(ξ/2)] ((uy wx - 
                                    ux wy) Cos[(ψ/2)] + (fz ux wx - fx uz wx + fz uy wy - fy uz wy - fx ux wz - fy uy wz - 
                                    fz uz wz) Sin[(ψ/2)])) + (Cos[(ψ/2)] ((-gy wx + 
                                    gx wy) Cos[(ξ/2)] - (gx uz wx + gy uz wy - gx ux wz - gy uy wz + 
                                    gz (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + (-(fx gz wx + fy gz wy - fx gx wz - fy gy wz + 
                                      fz (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fz (gz uy wx + gy uz wx - gz ux wy - gx uz wy + gy ux wz - 
                                        gx uy wz) + 
                                    fx (gy ux wx - gx uy wx - 
                                        gx ux wy - gy uy wy + 
                                        gz uz wy - gz uy wz - 
                                        gy uz wz) + 
                                    fy (gy uy wx - gz uz wx + 
                                        gy ux wy + gz ux wz + 
                                        gx (ux wx - uy wy + 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (wy Cos[(ψ/2)] + (-fz wx + 
                                    fx wz) Sin[(ψ/2)]) - 
                        Sin[(ξ/2)] ((uz wx - 
                                    ux wz) Cos[(ψ/2)] + (fx uy wx + fx ux wy + fz uz wy + fz uy wz - 
                                    fy (ux wx - uy wy + 
                                        uz wz)) Sin[(ψ/2)]))) - (Cos[(ψ/2)] (-(gx wx + gy wy + 
                                gz wz) Cos[(ξ/2)] + (-gz uy wx + gy uz wx + gz ux wy - gx uz wy - gy ux wz + 
                              gx uy wz) Sin[(ξ/2)]) + ((fz gy wx - fy gz wx - fz gx wy + fx gz wy + fy gx wz - 
                              fx gy wz) Cos[(ξ/2)] + (fy (-gy ux wx + gx uy wx + gx ux wy + gy uy wy + gz uz wy + 
                                    gz uy wz - gy uz wz) + 
                              fz (-gz ux wx + gx uz wx - 
                                    gz uy wy + gy uz wy + gx ux wz + 
                                    gy uy wz + gz uz wz) + 
                              fx (gy uy wx + gz uz wx + gy ux wy + 
                                    gz ux wz + 
                                    gx (ux wx - uy wy - 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) (-(Cos[(ξ/2)] (gy Cos[(ψ/2)] + (-fz gx + fx gz) Sin[(ψ/2)]) + 
                          Sin[(ξ/2)] ((gz ux - 
                                      gx uz) Cos[(ψ/2)] + (fy gx ux - fx gy ux - fx gx uy - fy gy uy - fz gz uy - fz gy uz + 
                                      fy gz uz) Sin[(ψ/2)])) (Cos[(ξ/2)] (wz Cos[(ψ/2)] + (fy wx - 
                                    fx wy) Sin[(ψ/2)]) + 
                        Sin[(ξ/2)] ((uy wx - 
                                    ux wy) Cos[(ψ/2)] + (fz ux wx - fx uz wx + fz uy wy - fy uz wy - fx ux wz - fy uy wz - 
                                    fz uz wz) Sin[(ψ/2)])) + (Cos[(ξ/2)] (gz Cos[(ψ/2)] + (fy gx - 
                                    fx gy) Sin[(ψ/2)]) - 
                        Sin[(ξ/2)] ((gy ux - 
                                    gx uy) Cos[(ψ/2)] + (fx gz ux + fy gz uy + fx gx uz + fy gy uz - 
                                    fz (gx ux + gy uy - 
                                        gz uz)) Sin[(ψ/2)])) (Cos[(ξ/2)] (wy Cos[(ψ/2)] + (-fz wx + 
                                    fx wz) Sin[(ψ/2)]) - 
                        Sin[(ξ/2)] ((uz wx - 
                                    ux wz) Cos[(ψ/2)] + (fx uy wx + fx ux wy + fz uz wy + fz uy wz - 
                                    fy (ux wx - uy wy + 
                                        uz wz)) Sin[(ψ/2)]))), (-(ux wx + uy wy + 
                          uz wz) Cos[(ψ/2)] Sin[(ξ/2)] + (-(fx wx + fy wy + fz wz) Cos[(ξ/2)] + (fz uy wx - 
                              fy uz wx - fz ux wy + fx uz wy + 
                              fy ux wz - 
                              fx uy wz) Sin[(ξ/2)]) Sin[(ψ/2)]) ((Cos[(ψ/2)] (gx Cos[(ξ/2)] + (-gz uy + gy uz) Sin[(ξ/2)]) + ((fz gy - 
                                    fy gz) Cos[(ξ/2)] - (fy gy ux + fz gz ux + fy gx uy + fz gx uz + 
                                    fx (gx ux - gy uy - 
                                        gz uz)) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ψ/2)] ((-gy wx + 
                                    gx wy) Cos[(ξ/2)] - (gx uz wx + gy uz wy - gx ux wz - gy uy wz + 
                                    gz (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + (-(fx gz wx + fy gz wy - fx gx wz - fy gy wz + 
                                      fz (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fz (gz uy wx + gy uz wx - gz ux wy - gx uz wy + gy ux wz - 
                                        gx uy wz) + 
                                    fx (gy ux wx - gx uy wx - 
                                        gx ux wy - gy uy wy + 
                                        gz uz wy - gz uy wz - 
                                        gy uz wz) + 
                                    fy (gy uy wx - gz uz wx + 
                                        gy ux wy + gz ux wz + 
                                        gx (ux wx - uy wy + 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) - (Cos[(ψ/2)] ((-gz wy + 
                                    gy wz) Cos[(ξ/2)] - (-gy uy wx - gz uz wx + gy ux wy + gz ux wz + 
                                    gx (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + ((fy gy wx + fz gz wx - fy gx wy - fz gx wz - 
                                    fx (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fx (gz uy wx - gy uz wx + gz ux wy + gx uz wy - gy ux wz - 
                                        gx uy wz) - 
                                    fy (gz ux wx + gx uz wx - 
                                        gz uy wy + gy uz wy - 
                                        gx ux wz + gy uy wz + 
                                        gz uz wz) + 
                                    fz (gx uy wx - gx ux wy + 
                                        gz uz wy + gz uy wz + 
                                        gy (ux wx + uy wy - 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (gz Cos[(ψ/2)] + (fy gx - 
                                    fx gy) Sin[(ψ/2)]) - 
                        Sin[(ξ/2)] ((gy ux - 
                                    gx uy) Cos[(ψ/2)] + (fx gz ux + fy gz uy + fx gx uz + fy gy uz - 
                                    fz (gx ux + gy uy - 
                                        gz uz)) Sin[(ψ/2)]))) + (Cos[(ψ/2)] (-(gx wx + gy wy + 
                                gz wz) Cos[(ξ/2)] + (-gz uy wx + gy uz wx + gz ux wy - gx uz wy - gy ux wz + 
                              gx uy wz) Sin[(ξ/2)]) + ((fz gy wx - fy gz wx - fz gx wy + fx gz wy + fy gx wz - 
                              fx gy wz) Cos[(ξ/2)] + (fy (-gy ux wx + gx uy wx + gx ux wy + gy uy wy + gz uz wy + 
                                    gz uy wz - gy uz wz) + 
                              fz (-gz ux wx + gx uz wx - 
                                    gz uy wy + gy uz wy + gx ux wz + 
                                    gy uy wz + gz uz wz) + 
                              fx (gy uy wx + gz uz wx + gy ux wy + 
                                    gz ux wz + 
                                    gx (ux wx - uy wy - 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) ((Cos[(ψ/2)] (wx Cos[(ξ/2)] + (uz wy - 
                                    uy wz) Sin[(ξ/2)]) + ((fz wy - fy wz) Cos[(ξ/2)] - (fy uy wx + fz uz wx + fy ux wy + 
                                    fz ux wz + 
                                    fx (ux wx - uy wy - 
                                        uz wz)) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (gz Cos[(ψ/2)] + (fy gx - 
                                    fx gy) Sin[(ψ/2)]) - 
                        Sin[(ξ/2)] ((gy ux - 
                                    gx uy) Cos[(ψ/2)] + (fx gz ux + fy gz uy + fx gx uz + fy gy uz - 
                                    fz (gx ux + gy uy - 
                                        gz uz)) Sin[(ψ/2)])) - (Cos[(ψ/2)] (gx Cos[(ξ/2)] + (-gz uy + 
                                    gy uz) Sin[(ξ/2)]) + ((fz gy - fy gz) Cos[(ξ/2)] - (fy gy ux + fz gz ux + fy gx uy + 
                                    fz gx uz + 
                                    fx (gx ux - gy uy - 
                                        gz uz)) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (wz Cos[(ψ/2)] + (fy wx - 
                                    fx wy) Sin[(ψ/2)]) + 
                        Sin[(ξ/2)] ((uy wx - 
                                    ux wy) Cos[(ψ/2)] + (fz ux wx - fx uz wx + fz uy wy - fy uz wy - fx ux wz - fy uy wz - 
                                    fz uz wz) Sin[(ψ/2)]))) - (-(gx ux + gy uy + 
                          gz uz) Cos[(ψ/2)] Sin[(ξ/2)] + (-(fx gx + fy gy + fz gz) Cos[(ξ/2)] + (-fz gy ux + 
                              fy gz ux + fz gx uy - fx gz uy - 
                              fy gx uz + 
                              fx gy uz) Sin[(ξ/2)]) Sin[(ψ/2)]) ((Cos[(ψ/2)] (wx Cos[(ξ/2)] + (uz wy - 
                                    uy wz) Sin[(ξ/2)]) + ((fz wy - fy wz) Cos[(ξ/2)] - (fy uy wx + fz uz wx + fy ux wy + 
                                    fz ux wz + 
                                    fx (ux wx - uy wy - 
                                        uz wz)) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ψ/2)] ((-gy wx + 
                                    gx wy) Cos[(ξ/2)] - (gx uz wx + gy uz wy - gx ux wz - gy uy wz + 
                                    gz (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + (-(fx gz wx + fy gz wy - fx gx wz - fy gy wz + 
                                      fz (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fz (gz uy wx + gy uz wx - gz ux wy - gx uz wy + gy ux wz - 
                                        gx uy wz) + 
                                    fx (gy ux wx - gx uy wx - 
                                        gx ux wy - gy uy wy + 
                                        gz uz wy - gz uy wz - 
                                        gy uz wz) + 
                                    fy (gy uy wx - gz uz wx + 
                                        gy ux wy + gz ux wz + 
                                        gx (ux wx - uy wy + 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) - (Cos[(ψ/2)] ((-gz wy + 
                                    gy wz) Cos[(ξ/2)] - (-gy uy wx - gz uz wx + gy ux wy + gz ux wz + 
                                    gx (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + ((fy gy wx + fz gz wx - fy gx wy - fz gx wz - 
                                    fx (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fx (gz uy wx - gy uz wx + gz ux wy + gx uz wy - gy ux wz - 
                                        gx uy wz) - 
                                    fy (gz ux wx + gx uz wx - 
                                        gz uy wy + gy uz wy - 
                                        gx ux wz + gy uy wz + 
                                        gz uz wz) + 
                                    fz (gx uy wx - gx ux wy + 
                                        gz uz wy + gz uy wz + 
                                        gy (ux wx + uy wy - 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (wz Cos[(ψ/2)] + (fy wx - 
                                    fx wy) Sin[(ψ/2)]) + 
                        Sin[(ξ/2)] ((uy wx - 
                                    ux wy) Cos[(ψ/2)] + (fz ux wx - fx uz wx + fz uy wy - fy uz wy - fx ux wz - fy uy wz - 
                                    fz uz wz) Sin[(ψ/2)]))), -(-(ux wx + uy wy + 
                            uz wz) Cos[(ψ/2)] Sin[(ξ/2)] + (-(fx wx + fy wy + fz wz) Cos[(ξ/2)] + (fz uy wx - 
                                fy uz wx - fz ux wy + fx uz wy + 
                                fy ux wz - 
                                fx uy wz) Sin[(ξ/2)]) Sin[(ψ/2)]) ((Cos[(ψ/2)] (gx Cos[(ξ/2)] + (-gz uy + gy uz) Sin[(ξ/2)]) + ((fz gy - 
                                    fy gz) Cos[(ξ/2)] - (fy gy ux + fz gz ux + fy gx uy + fz gx uz + 
                                    fx (gx ux - gy uy - 
                                        gz uz)) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ψ/2)] ((gz wx - 
                                    gx wz) Cos[(ξ/2)] - (gx uy wx - gx ux wy - gz uz wy + gz uy wz + 
                                    gy (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + (-(fx gy wx - fx gx wy - fz gz wy + fz gy wz + 
                                      fy (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fy (-gy uz wx + gx uz wy - gz (uy wx + ux wy) + gy ux wz + 
                                        gx uy wz) - 
                                    fz (-gy uy wx + gz uz wx + 
                                        gy ux wy + gz ux wz + 
                                        gx (ux wx + uy wy - 
                                        uz wz)) + 
                                    fx (gx uz wx + gy uz wy + 
                                        gx ux wz - gy uy wz + 
                                        gz (-ux wx + uy wy + 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) - (Cos[(ψ/2)] ((-gz wy + 
                                    gy wz) Cos[(ξ/2)] - (-gy uy wx - gz uz wx + gy ux wy + gz ux wz + 
                                    gx (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + ((fy gy wx + fz gz wx - fy gx wy - fz gx wz - 
                                    fx (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fx (gz uy wx - gy uz wx + gz ux wy + gx uz wy - gy ux wz - 
                                        gx uy wz) - 
                                    fy (gz ux wx + gx uz wx - 
                                        gz uy wy + gy uz wy - 
                                        gx ux wz + gy uy wz + 
                                        gz uz wz) + 
                                    fz (gx uy wx - gx ux wy + 
                                        gz uz wy + gz uy wz + 
                                        gy (ux wx + uy wy - 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (gy Cos[(ψ/2)] + (-fz gx + 
                                    fx gz) Sin[(ψ/2)]) + 
                        Sin[(ξ/2)] ((gz ux - 
                                    gx uz) Cos[(ψ/2)] + (fy gx ux - fx gy ux - fx gx uy - fy gy uy - fz gz uy - fz gy uz + 
                                    fy gz uz) Sin[(ψ/2)]))) - (Cos[(ψ/2)] (-(gx wx + gy wy + 
                                gz wz) Cos[(ξ/2)] + (-gz uy wx + gy uz wx + gz ux wy - gx uz wy - gy ux wz + 
                              gx uy wz) Sin[(ξ/2)]) + ((fz gy wx - fy gz wx - fz gx wy + fx gz wy + fy gx wz - 
                              fx gy wz) Cos[(ξ/2)] + (fy (-gy ux wx + gx uy wx + gx ux wy + gy uy wy + gz uz wy + 
                                    gz uy wz - gy uz wz) + 
                              fz (-gz ux wx + gx uz wx - 
                                    gz uy wy + gy uz wy + gx ux wz + 
                                    gy uy wz + gz uz wz) + 
                              fx (gy uy wx + gz uz wx + gy ux wy + 
                                    gz ux wz + 
                                    gx (ux wx - uy wy - 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) ((Cos[(ψ/2)] (wx Cos[(ξ/2)] + (uz wy - 
                                    uy wz) Sin[(ξ/2)]) + ((fz wy - fy wz) Cos[(ξ/2)] - (fy uy wx + fz uz wx + fy ux wy + 
                                    fz ux wz + 
                                    fx (ux wx - uy wy - 
                                        uz wz)) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (gy Cos[(ψ/2)] + (-fz gx + 
                                    fx gz) Sin[(ψ/2)]) + 
                        Sin[(ξ/2)] ((gz ux - 
                                    gx uz) Cos[(ψ/2)] + (fy gx ux - fx gy ux - fx gx uy - fy gy uy - fz gz uy - fz gy uz + 
                                    fy gz uz) Sin[(ψ/2)])) - (Cos[(ψ/2)] (gx Cos[(ξ/2)] + (-gz uy + 
                                    gy uz) Sin[(ξ/2)]) + ((fz gy - fy gz) Cos[(ξ/2)] - (fy gy ux + fz gz ux + fy gx uy + 
                                    fz gx uz + 
                                    fx (gx ux - gy uy - 
                                        gz uz)) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (wy Cos[(ψ/2)] + (-fz wx + 
                                    fx wz) Sin[(ψ/2)]) - 
                        Sin[(ξ/2)] ((uz wx - 
                                    ux wz) Cos[(ψ/2)] + (fx uy wx + fx ux wy + fz uz wy + fz uy wz - 
                                    fy (ux wx - uy wy + 
                                        uz wz)) Sin[(ψ/2)]))) + (-(gx ux + gy uy + 
                          gz uz) Cos[(ψ/2)] Sin[(ξ/2)] + (-(fx gx + fy gy + fz gz) Cos[(ξ/2)] + (-fz gy ux + 
                              fy gz ux + fz gx uy - fx gz uy - 
                              fy gx uz + 
                              fx gy uz) Sin[(ξ/2)]) Sin[(ψ/2)]) ((Cos[(ψ/2)] (wx Cos[(ξ/2)] + (uz wy - 
                                    uy wz) Sin[(ξ/2)]) + ((fz wy - fy wz) Cos[(ξ/2)] - (fy uy wx + fz uz wx + fy ux wy + 
                                    fz ux wz + 
                                    fx (ux wx - uy wy - 
                                        uz wz)) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ψ/2)] ((gz wx - 
                                    gx wz) Cos[(ξ/2)] - (gx uy wx - gx ux wy - gz uz wy + gz uy wz + 
                                    gy (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + (-(fx gy wx - fx gx wy - fz gz wy + fz gy wz + 
                                      fy (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fy (-gy uz wx + gx uz wy - gz (uy wx + ux wy) + gy ux wz + 
                                        gx uy wz) - 
                                    fz (-gy uy wx + gz uz wx + 
                                        gy ux wy + gz ux wz + 
                                        gx (ux wx + uy wy - 
                                        uz wz)) + 
                                    fx (gx uz wx + gy uz wy + 
                                        gx ux wz - gy uy wz + 
                                        gz (-ux wx + uy wy + 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) - (Cos[(ψ/2)] ((-gz wy + 
                                    gy wz) Cos[(ξ/2)] - (-gy uy wx - gz uz wx + gy ux wy + gz ux wz + 
                                    gx (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + ((fy gy wx + fz gz wx - fy gx wy - fz gx wz - 
                                    fx (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fx (gz uy wx - gy uz wx + gz ux wy + gx uz wy - gy ux wz - 
                                        gx uy wz) - 
                                    fy (gz ux wx + gx uz wx - 
                                        gz uy wy + gy uz wy - 
                                        gx ux wz + gy uy wz + 
                                        gz uz wz) + 
                                    fz (gx uy wx - gx ux wy + 
                                        gz uz wy + gz uy wz + 
                                        gy (ux wx + uy wy - 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (wy Cos[(ψ/2)] + (-fz wx + 
                                    fx wz) Sin[(ψ/2)]) - 
                        Sin[(ξ/2)] ((uz wx - 
                                    ux wz) Cos[(ψ/2)] + (fx uy wx + fx ux wy + fz uz wy + fz uy wz - 
                                    fy (ux wx - uy wy + 
                                        uz wz)) Sin[(ψ/2)]))), (Cos[(ξ/2)] (wz Cos[(ψ/2)] + (fy wx - fx wy) Sin[(ψ/2)]) + 
                  Sin[(ξ/2)] ((uy wx - 
                              ux wy) Cos[(ψ/2)] + (fz ux wx - 
                              fx uz wx + fz uy wy - fy uz wy - 
                              fx ux wz - fy uy wz - 
                              fz uz wz) Sin[(ψ/2)])) ((Cos[(ψ/2)] (gx Cos[(ξ/2)] + (-gz uy + 
                                    gy uz) Sin[(ξ/2)]) + ((fz gy - fy gz) Cos[(ξ/2)] - (fy gy ux + fz gz ux + fy gx uy + 
                                    fz gx uz + 
                                    fx (gx ux - gy uy - 
                                        gz uz)) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ψ/2)] ((gz wx - 
                                    gx wz) Cos[(ξ/2)] - (gx uy wx - gx ux wy - gz uz wy + gz uy wz + 
                                    gy (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + (-(fx gy wx - fx gx wy - fz gz wy + fz gy wz + 
                                      fy (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fy (-gy uz wx + gx uz wy - gz (uy wx + ux wy) + gy ux wz + 
                                        gx uy wz) - 
                                    fz (-gy uy wx + gz uz wx + 
                                        gy ux wy + gz ux wz + 
                                        gx (ux wx + uy wy - 
                                        uz wz)) + 
                                    fx (gx uz wx + gy uz wy + 
                                        gx ux wz - gy uy wz + 
                                        gz (-ux wx + uy wy + 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) - (Cos[(ψ/2)] ((-gz wy + 
                                    gy wz) Cos[(ξ/2)] - (-gy uy wx - gz uz wx + gy ux wy + gz ux wz + 
                                    gx (ux wx + uy wy + 
                                        
                                        uz wz)) Sin[(ξ/2)]) + ((fy gy wx + fz gz wx - fy gx wy - fz gx wz - 
                                    fx (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fx (gz uy wx - gy uz wx + gz ux wy + gx uz wy - gy ux wz - 
                                        gx uy wz) - 
                                    fy (gz ux wx + gx uz wx - 
                                        gz uy wy + gy uz wy - 
                                        gx ux wz + gy uy wz + 
                                        gz uz wz) + 
                                    fz (gx uy wx - gx ux wy + 
                                        gz uz wy + gz uy wz + 
                                        gy (ux wx + uy wy - 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (gy Cos[(ψ/2)] + (-fz gx + 
                                    fx gz) Sin[(ψ/2)]) + 
                        Sin[(ξ/2)] ((gz ux - 
                                    gx uz) Cos[(ψ/2)] + (fy gx ux - fx gy ux - fx gx uy - fy gy uy - fz gz uy - fz gy uz + 
                                    fy gz uz) Sin[(ψ/2)]))) + (Cos[(ψ/2)] ((-gy wx + 
                              gx wy) Cos[(ξ/2)] - (gx uz wx + 
                              gy uz wy - gx ux wz - gy uy wz + 
                              gz (ux wx + uy wy + 
                                    uz wz)) Sin[(ξ/2)]) + (-(fx gz wx + fy gz wy - fx gx wz - fy gy wz + 
                                fz (gx wx + gy wy + 
                                      gz wz)) Cos[(ξ/2)] + (fz (gz uy wx + gy uz wx - gz ux wy - gx uz wy + gy ux wz - 
                                    gx uy wz) + 
                              fx (gy ux wx - gx uy wx - gx ux wy - 
                                    gy uy wy + gz uz wy - gz uy wz - 
                                    gy uz wz) + 
                              fy (gy uy wx - gz uz wx + gy ux wy + 
                                    gz ux wz + 
                                    gx (ux wx - uy wy + 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) ((Cos[(ψ/2)] (wx Cos[(ξ/2)] + (uz wy - 
                                    uy wz) Sin[(ξ/2)]) + ((fz wy - fy wz) Cos[(ξ/2)] - (fy uy wx + fz uz wx + fy ux wy + 
                                    fz ux wz + 
                                    fx (ux wx - uy wy - 
                                        uz wz)) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (gy Cos[(ψ/2)] + (-fz gx + 
                                    fx gz) Sin[(ψ/2)]) + 
                        Sin[(ξ/2)] ((gz ux - 
                                    gx uz) Cos[(ψ/2)] + (fy gx ux - fx gy ux - fx gx uy - fy gy uy - fz gz uy - fz gy uz + 
                                    
                                    fy gz uz) Sin[(ψ/2)])) - (Cos[(ψ/2)] (gx Cos[(ξ/2)] + (-gz uy + 
                                    gy uz) Sin[(ξ/2)]) + ((fz gy - fy gz) Cos[(ξ/2)] - (fy gy ux + fz gz ux + fy gx uy + 
                                    fz gx uz + 
                                    fx (gx ux - gy uy - 
                                        gz uz)) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (wy Cos[(ψ/2)] + (-fz wx + 
                                    fx wz) Sin[(ψ/2)]) - 
                        Sin[(ξ/2)] ((uz wx - 
                                    ux wz) Cos[(ψ/2)] + (fx uy wx + fx ux wy + fz uz wy + fz uy wz - 
                                    fy (ux wx - uy wy + 
                                        uz wz)) Sin[(ψ/2)]))) - (Cos[(ξ/2)] (gz Cos[(ψ/2)] + (fy gx - 
                              fx gy) Sin[(ψ/2)]) - 
                  Sin[(ξ/2)] ((gy ux - 
                              gx uy) Cos[(ψ/2)] + (fx gz ux + 
                              fy gz uy + fx gx uz + fy gy uz - 
                              fz (gx ux + gy uy - 
                                    gz uz)) Sin[(ψ/2)])) ((Cos[(ψ/2)] (wx Cos[(ξ/2)] + (uz wy - 
                                    uy wz) Sin[(ξ/2)]) + ((fz wy - fy wz) Cos[(ξ/2)] - (fy uy wx + fz uz wx + fy ux wy + 
                                    fz ux wz + 
                                    fx (ux wx - uy wy - 
                                        uz wz)) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ψ/2)] ((gz wx - 
                                    gx wz) Cos[(ξ/2)] - (gx uy wx - gx ux wy - gz uz wy + gz uy wz + 
                                    gy (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + (-(fx gy wx - fx gx wy - fz gz wy + fz gy wz + 
                                      fy (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fy (-gy uz wx + gx uz wy - gz (uy wx + ux wy) + gy ux wz + 
                                        gx uy wz) - 
                                    fz (-gy uy wx + gz uz wx + 
                                        gy ux wy + gz ux wz + 
                                        gx (ux wx + uy wy - 
                                        uz wz)) + 
                                    fx (gx uz wx + gy uz wy + 
                                        gx ux wz - gy uy wz + 
                                        gz (-ux wx + uy wy + 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) - (Cos[(ψ/2)] ((-gz wy + 
                                    gy wz) Cos[(ξ/2)] - (-gy uy wx - gz uz wx + gy ux wy + gz ux wz + 
                                    gx (ux wx + uy wy + 
                                        uz wz)) Sin[(ξ/2)]) + ((fy gy wx + fz gz wx - fy gx wy - fz gx wz - 
                                    fx (gx wx + gy wy + 
                                        gz wz)) Cos[(ξ/2)] + (fx (gz uy wx - gy uz wx + gz ux wy + gx uz wy - gy ux wz - 
                                        gx uy wz) - 
                                    fy (gz ux wx + gx uz wx - 
                                        gz uy wy + gy uz wy - 
                                        gx ux wz + gy uy wz + 
                                        gz uz wz) + 
                                    fz (gx uy wx - gx ux wy + 
                                        gz uz wy + gz uy wz + 
                                        gy (ux wx + uy wy - 
                                        uz wz))) Sin[(ξ/2)]) Sin[(ψ/2)]) (Cos[(ξ/2)] (wy Cos[(ψ/2)] + (-fz wx + 
                                    fx wz) Sin[(ψ/2)]) - 
                        Sin[(ξ/2)] ((uz wx - 
                                    ux wz) Cos[(ψ/2)] + (fx uy wx + fx ux wy + fz uz wy + fz uy wz - 
                                    fy (ux wx - uy wy + 
                                        uz wz)) Sin[(ψ/2)])))}}

We find the vector expression of the matrix. The vectors will be the rows of the matrix.

f = {fx, fy, fz}
g = {gx, gy, gz}
w = {wx, wy, wz}
u = {ux, uy, uz}

I AM HERE

Simplify[
      invmatN[[1]] - 
        Flatten[{ 
            Cos[(ψ/2)] g\[Cross](g\[Cross]w) +  
              Sin[(ψ/2)] (g\[Cross](g\[Cross]w))\[Cross]
                  f, -Sin[(ψ/2)] (g\[Cross](g\[Cross]w)) . 
                f}]]

Simplify[
      invmatN2[[2]] - 
        Flatten[{ 
            Cos[(ψ/2)] ((g\[Cross]w)\[Cross]w) +  
              Sin[(ψ/2)] ((g\[Cross]w)\[Cross]w)\[Cross]
                  f, -Sin[(ψ/2)] ((g\[Cross]w)\[Cross]w) . 
                f}]]

Simplify[
      invmatN2[[3]] - 
        Flatten[{ -Cos[(ψ/2)] (g\[Cross]w) -  
              Sin[(ψ/2)] (g\[Cross]w)\[Cross]f, 
            Sin[(ψ/2)] (g\[Cross]w) . f}]]

Simplify[
      invmatN2[[4]] - 
        Flatten[{ -Cos[(ψ/2)] (g . w) (g\[Cross]w) -  
              Sin[(ψ/2)] ((g . w) (g\[Cross]w)\[Cross]f + 
                    f (g\[Cross]w) . (g\[Cross]w)), (- 
                    Cos[(ψ/2)] (g\[Cross]w) . (g\[Cross]w) +  
                Sin[(ψ/2)] (g . 
                      w) (f . (g\[Cross]w)))}]]

invmatS = 
      invmatN2/(-(gz^2) ((wx^2) + (wy^2)) + 2 gx gz wx wz + 
            2 gy wy (gx wx + gz wz) - (gy^2) ((wx^2) + (wz^2)) - 
            (gx^2) ((wy^2) + (wz^2)))

sols = Simplify[invmatS . {px, py, pz, pw}]

ctsp = sols[[1]]
stcp = sols[[2]]
stsp = sols[[3]]
ctcp = sols[[4]]

Extra condition for the values to be angles: equating tangents

extraeq = 
      Simplify[Numerator[stsp]*Numerator[ctcp] - 
          Numerator[stcp]*Numerator[ctsp]]

extraeq2 = 
      Collect[extraeq, {Cos[(ψ/2)]^2, Sin[(ψ/2)]^2, 
          Cos[(ψ/2)] Sin[(ψ/2)]}, Simplify]

Simplify[
      Coefficient[extraeq2, 
          Cos[(ψ/2)]^2] - (g\[Cross]w) . (g\[Cross]w) 
          p . (pw g\[Cross]w - g\[Cross](w\[Cross]p))]

Simplify[
      Simplify[Coefficient[extraeq2, Cos[(ψ/2)] Sin[(ψ/2)]] /. 
            pw^2 -> 
              1 - px^2 - py^2 - pz^2] + (g\[Cross]w) . (g\[Cross]
                w) (f . (g\[Cross]w) - 
              2  pw (g\[Cross]p) . (f\[Cross]w) - 
              2 ((g\[Cross]p)\[Cross]p) . (f\[Cross]
                      w))]

Simplify[
      Simplify[Coefficient[extraeq2, Cos[(ψ/2)] Sin[(ψ/2)]] /. 
            pw^2 -> 
              1 - px^2 - py^2 - pz^2] + (g\[Cross]w) . (g\[Cross]
                w) (f . (g\[Cross]w) - 
              2  pw (g\[Cross](w\[Cross]f)) . p + 
              2 (p . (g\[Cross]((w\[Cross]f)\[Cross]
                            p))))]

Simplify[
      Coefficient[extraeq2, Sin[(ψ/2)]^2] /. 
        pw^2 -> 1 - px^2 - py^2 - pz^2]

FullSimplify[
      Simplify[Coefficient[extraeq2, Sin[(ψ/2)]^2] /. 
            pw^2 -> 
              1 - px^2 - py^2 - 
                pz^2] - ((g\[Cross]w) . (g\[Cross]
                  w)) ((f\[Cross]g) . (f\[Cross]w) + 
              pw (-f . w f . (g\[Cross]p) + (f\[Cross]
                          w) . (f\[Cross](g\[Cross]p))) - 
              2 ((g\[Cross]p)\[Cross]p) . ((w\[Cross]f)\[Cross]
                      f) + f . (g\[Cross]p) f . (w\[Cross]p) - 
              f . g f . (p\[Cross](p\[Cross]w)) + 
              f . p f . (g\[Cross](p\[Cross]w)))]

Simplify[(extraeq2 /. 
            pw^2 -> 
              1 - px^2 - py^2 - 
                pz^2) - (((g\[Cross]w) . (g\[Cross]
                    w)) (Cos[(ψ/2)]^2 (pw (g\[Cross]
                              w) . p - (g\[Cross](w\[Cross]p)) . 
                        p) + 
                Cos[(ψ/2)] 
                  Sin[(ψ/2)] (-(g\[Cross]w) . f + 
                      2 pw (g\[Cross](w\[Cross]f)) . p - 
                      2 (g\[Cross]((w\[Cross]f)\[Cross]p)) . 
                          p) + Sin[(ψ/2)]^2 ((f\[Cross]
                            g) . (f\[Cross]w) + 
                      pw ((f . w) (g\[Cross]f) . 
                                p + (f\[Cross]
                                  w) . (f\[Cross](g\[Cross]p))) - 
                      2 ((g\[Cross]p)\[Cross]
                              p) . ((w\[Cross]f)\[Cross]
                              f) + ((f\[Cross]g) . 
                            p) ((f\[Cross]w) . p) + (f . 
                            g) (w\[Cross]p) . (f\[Cross]p) + (f . 
                            p ) (g\[Cross]f) . (w\[Cross]
                              p))))]

extraeqNEW = ((g\[Cross]w) . (g\[Cross]
                w)) (Cos[(ψ/2)]^2 (pw (g\[Cross]w) . 
                      p - (g\[Cross](w\[Cross]p)) . p) + 
            Cos[(ψ/2)] 
              Sin[(ψ/2)] (-(g\[Cross]w) . f + 
                  2 pw (g\[Cross](w\[Cross]f)) . p - 
                  2 (g\[Cross]((w\[Cross]f)\[Cross]p)) . 
                      p) + Sin[(ψ/2)]^2 ((f\[Cross]
                        g) . (f\[Cross]w) + 
                  pw ((f . w) (g\[Cross]f) . 
                            p + (f\[Cross]
                              w) . (f\[Cross](g\[Cross]p))) - 
                  2 ((g\[Cross]p)\[Cross]
                          p) . ((w\[Cross]f)\[Cross]
                          f) + ((f\[Cross]g) . 
                        p) ((f\[Cross]w) . p) + (f . 
                        g) (w\[Cross]p) . (f\[Cross]p) + (f . 
                        p ) (g\[Cross]f) . (w\[Cross]
                          p)))

The advantage of extraeqNEW is that we can eliminate the (gxw).(gxw) factor that we had also in the denominator.

extraeqFinal = Cos[(ψ/2)]^2 (pw (g\[Cross]w) . 
                  p - (g\[Cross](w\[Cross]p)) . p) + 
        Cos[(ψ/2)] 
          Sin[(ψ/2)] (-(g\[Cross]w) . f + 
              2 pw (g\[Cross](w\[Cross]f)) . p - 
              2 (g\[Cross]((w\[Cross]f)\[Cross]p)) . 
                  p) + Sin[(ψ/2)]^2 ((f\[Cross]
                    g) . (f\[Cross]w) + 
              pw ((f . w) (g\[Cross]f) . 
                        p + (f\[Cross]
                          w) . (f\[Cross](g\[Cross]p))) - 
              2 ((g\[Cross]p)\[Cross]p) . ((w\[Cross]f)\[Cross]
                      f) + ((f\[Cross]g) . p) ((f\[Cross]w) . 
                    p) + (f . 
                    g) (w\[Cross]p) . (f\[Cross]p) + (f . 
                    p ) (g\[Cross]f) . (w\[Cross]p))

This is useless, I don't know why it is here. Anyway, we can compare the above condition with the one we get when we have only two revolute joint variables.

extraeqOLD = 
      gz ((-pw py + px pz) wx + (pw px + 
                    py pz) wy + (-(px^2) - (py^2)) wz) + 
        gx ((-(py^2) - (pz^2)) wx + (px py - 
                    pw pz) wy + (pw py + px pz) wz) + 
        gy ((px py + pw pz) wx + (-(px^2) - 
                    (pz^2)) wy + (-pw px + 
                    py pz) wz)

4. Eliminating ψ from the translations

secRRR = 
Collect[{Qrrrsep2[[5]], 
            Qrrrsep2[[6]], Qrrrsep2[[7]], Qrrrsep2[[8]]} - {px0, 
            py0, pz0, pw0}, {Cos[(θ/2)] Sin[(ϕ/2)], 
          Cos[(ϕ/2)] Sin[(θ/2)], Sin[(θ/2)] Sin[(ϕ/2)], 
          Cos[(θ/2)] Cos[(ϕ/2)]}]

secRRR2 = 
      secRRR /. {Cos[(θ/2)] Sin[(ϕ/2)] -> ctsp,  
          Cos[(ϕ/2)] Sin[(θ/2)] -> stcp, 
          Sin[(θ/2)] Sin[(ϕ/2)] -> stsp,  
          Cos[(θ/2)] Cos[(ϕ/2)] -> ctcp}

Here we clear denominators, what means that we multiply p0 times (gxw).(gxw). Here we can solve for the half angles or for the angle:

1. Solving for the values of Cos[ψ], Sin[ψ]

secRRR3 = 
        Simplify[
          Together[
              secRRR2]*((gz^2) ((wx^2) + (wy^2)) - 2 gx gz wx wz - 
                2 gy wy (gx wx + gz wz) + 
                (gy^2) ((wx^2) + (wz^2)) + 
                (gx^2) ((wy^2) + (wz^2)))];

secRRR4 = 
        Collect[secRRR3, {Cos[ψ], Sin[ψ], px0, py0, 
            pz0}];

Eliminate the last useless equation,

secRRR5 = {secRRR4[[1]], secRRR4[[2]], 
        secRRR4[[3]]}

Add the extra equation:

extraeq2 = Simplify[Expand[extraeq]]

We can see that there is a (gxw).(gxw) mulitplying the whole equation. Eliminating it,

g = {gx, gy, gz}
w = {wx, wy, wz}

Simplify[((gz^2) ((wx^2) + (wy^2)) - 2 gx gz wx wz - 
            2 gy wy (gx wx + gz wz) + (gy^2) ((wx^2) + (wz^2)) + 
            (gx^2) ((wy^2) + (wz^2))) - (g\[Cross]w) . (g\[Cross]
              w)]

extraeq3 = 
      Simplify[extraeq2/((1/2) ((gz^2) ((wx^2) + (wy^2)) - 
                  2 gx gz wx wz - 2 gy wy (gx wx + gz wz) + 
                  (gy^2) ((wx^2) + (wz^2)) + 
                  (gx^2) ((wy^2) + (wz^2))) )]

secRRR6 = Append[secRRR5, extraeq3];

All equations are of the form A cos(psi) +B sin(psi) +C

Solving for psi

col1 = Coefficient[secRRR6, Cos[ψ]];

col2 = Coefficient[secRRR6, Sin[ψ]];

col3 = 
        Coefficient[secRRR6*Cos[ψ]*Sin[ψ], 
          Cos[ψ]*Sin[ψ]];

matpsi = Transpose[{col1, col2, col3}];

This is a 4x3 matrix in{ cos(psi), sin(psi), 1}

symatpsi = {{A1, B1, C1}, {A2, B2, C2}, {A3, B3, C3}, {A4, B4, 
          C4}}

The 4 minors have to be equal to zero,

minspsi = Minors[symatpsi, 3]

But the fact is that I need only three equations! Which one is useless?

2. Solving for the half angle Cos[ψ/2]^2, Sin[ψ/2]^2, Cos[ψ/2]Sin[ψ/2]

secRRR3b = 
        Collect[Together[
              secRRR2]*((gz^2) ((wx^2) + (wy^2)) - 2 gx gz wx wz - 
                2 gy wy (gx wx + gz wz) + 
                (gy^2) ((wx^2) + (wz^2)) + 
                (gx^2) ((wy^2) + (wz^2))), {Cos[ψ]^2, Sin[ψ]^2, 
            Cos[ψ]*Sin[ψ], px0, py0, pz0}];

Eliminate the last useless equation,

secRRR5b = {secRRR3b[[1]], secRRR3b[[2]], 
          secRRR3b[[3]]};

Add the extra equation:

secRRR6b = Append[secRRR5b, extraeqFinal];

All equations are of the form A cos(psi) +B sin(psi) +C

Solving for psi

col1 = Simplify[Coefficient[secRRR6b, Cos[ψ/2]^2]]

col2 = Simplify[Coefficient[secRRR6b, Sin[ψ/2]^2]];

col3 = 
        Simplify[
          Coefficient[secRRR6b, Cos[ψ/2] Sin[ψ/2]]];

col4 = 
      Simplify[Coefficient[
          secRRR6b*Cos[ψ/2]^2*Sin[ψ/2]^2*Cos[ψ/2] 
            Sin[ψ/2], 
          Cos[ψ/2]^2*Sin[ψ/2]^2*Cos[ψ/2] 
            Sin[ψ/2]]]

matpsib = Transpose[{col1, col2, col3, col4}];

This is a 4x3 matrix in{ cos(psi), sin(psi), 1}

symatpsib = {{A1, B1, C1, D1}, {A2, B2, C2, D2}, {A3, B3, C3, D3}, {A4, 
          B4, C4, 0}}

The determinant must be zero

detpsib = Det[symatpsib]

In addition, we eliminate the last row and solve for the three psi variables

symatpsib3 = Drop[symatpsib, -1]

colb1 = 
      Flatten[Take[Transpose[symatpsib3], 1]]
colb2 = 
      Flatten[Take[Transpose[symatpsib3], {2, 2}]]
colb3 = 
      Flatten[Take[Transpose[symatpsib3], {3, 3}]]
indep = Flatten[Take[Transpose[symatpsib3], {4, 4}]]

cos2psi = 
      Det[{-indep, colb2, colb3}]/Det[{colb1, colb2, colb3}]

sin2psi = 
      Det[{colb1, -indep, colb3}]/Det[{colb1, colb2, colb3}]

cossinpsi = 
      Det[{colb1, colb2, -indep}]/Det[{colb1, colb2, colb3}]

The relation among them gives us the last two design equations,

relpsi1 = 
      Numerator[cos2psi] + Numerator[sin2psi] - 
        Denominator[cos2psi]

relpsi2 = 
      Numerator[cos2psi]*
          Numerator[sin2psi] - (Numerator[cossinpsi])^2

The set of design equations

designsP = {detpsib, relpsi1, relpsi2}

constraints = {gx^2 + gy^2 + gz^2 - 1, gx gx0 + gy gy0 + gz gz0, 
        wx^2 + wy^2 + wz^2 - 1, wx wx0 + wy wy0 + wz wz0, 
        fx^2 + fy^2 + fz^2 - 1, fx fx0 + fy fy0 + fz fz0}

Where  the ABCD coefficients are:

Remember to set the ABCD to nonassigned before doing further calculations above!

A1 =. ;  A2 =. ;  A3 =. ;  A4 =. ; B1 =. ;  B2 =. ;  B3 =. ;  
    B4 =. ;
C1 =. ;  C2 =. ;  C3 =. ;  C4 =. ; D1 =. ; D2 =. ;  D3 =. ;  
    D4 =. ;

Aes = col1;
Bes = col2;
Ces = col3;
Des = col4;

A1 = Aes[[1]]
A2 = Aes[[2]]
A3 = Aes[[3]]
A4 = Aes[[4]]

B1 = Bes[[1]]
B2 = Bes[[2]]
B3 = Bes[[3]]
B4 = Bes[[4]]

C1 = Ces[[1]]
C2 = Ces[[2]]
C3 = Ces[[3]]
C4 = Ces[[4]]

D1 = Des[[1]]
D2 = Des[[2]]
D3 = Des[[3]]
D4 = Des[[4]]

Create the data file

puf = OpenWrite["RRR.dat", PageWidth -> 200]

This writes out an expression, then closes the stream.

\("********The ABCD coefficients for quaternion P*******" >>> "RRR.dat"\)

"********{A1,A2,A3,A4} coefficients*******" >>> "RRR.dat"

{A1, A2, A3, A4} >>> "RRR.dat"

"********{B1,B2,B3,B4} coefficients*******" >>> "RRR.dat"

{B1, B2, B3, B4} >>> "RRR.dat"

"********{C1,C2,C3,C4} coefficients*******" >>> "RRR.dat"

{C1, C2, C3, C4} >>> "RRR.dat"

"********{D1,D2,D3,D4} coefficients*******" >>> "RRR.dat"

{D1, D2, D3, D4} >>> "RRR.dat"

A1 =. ;  A2 =. ;  A3 =. ;  A4 =. ; B1 =. ;  B2 =. ;  B3 =. ;  
    B4 =. ;
C1 =. ;  C2 =. ;  C3 =. ;  C4 =. ; D1 =. ; D2 =. ;  D3 =. ;  
    D4 =. ;

"********Design equations per dual quaternion*******" >> "RRR.dat\
"

designsP[[1]] >>> "RRR.dat"

"***************" >> "RRR.dat"

designsP[[2]] >>> "RRR.dat"

"***************" >> "RRR.dat"

designsP[[3]] >>> "RRR.dat"

Close[puf]

2. Analysis

1. The goal relative positions in quaternion form

qS12 = {{0.7715500492004921, -0.5789752698012732,  
          0.11937933787826784,  0.2350462342285181 }, {1.212670959545422 , 
          1.5419206725328352,  
          2.193473119045673, -1.2965864810179322}}
qS13 = {{-0.8618183013324527,  0.1609974994520099,  
          0.23907761282003118,  
          0.41736185224205735}, {-0.7029087510648233, -0.9464044712951724, -0.8121132551694119, -0.621170269174268}}
qS14 = {{0.11565089621249165,  
          0.13118300091238583, -0.15381589878830845 , 
          0.9725001592579887}, {-0.6260162759078233, -0.5552729243327017 , 0.38365035396299163,  
          0.21002899995699775}}
qS15 = {{0.49663418478173815,  
          0.3247529754196911, -0.8001059749341218,  
          0.08786592250101499}, {-1.483877888770251, -1.6344082244303662, -1.6073021186135898, -0.20814223042566574}}
qS16 = {{-0.6855583163603962, -0.6883409138591217 , 
          0.20377277954552686,  
          0.12113313127703194}, {-0.4368010925800464,  
          0.6666342614335875,  
          1.5452302035229777, -1.283355237821838}}
qS17 = {{-0.36484932165685185, -0.025692428710243558 , 
          0.9253737571645607 , 
          0.09954034932954853}, {-0.020295628254378938, -0.6752367107184916,  
          0.13346334658310335, -1.4894140589602725}}
qS18 = {{0.21965525972796746, -0.21603165047330938,  
          0.2921698301940352,  
          0.9053831692670268}, {-0.3953600920034315, -0.2057223290230593, -1.2307076742566916,  
          0.443984442463763}}
\(\(\(qS19 = {{\(-0.6337061170068199\), \ 0.1993190283532657, \ 
          0.7330754780625363\ , 
          0.1459069075398173}, {\(-0.27956290318767\)\ , 
          2.4976713517202773, \(-0.6525305958151515\), \(-1.3477084371594918\)}}\)\(
\)\(
\)
    \)\)

The screw axes and the rotation and translation:

b12 = {qS12[[1, 1]], qS12[[1, 2]], 
        qS12[[1, 3]]}
point12 = (b12/
            Sqrt[b12 . b12])\[Cross]({qS12[[2, 1]], qS12[[2, 2]], 
              qS12[[2, 3]]}/Sqrt[b12 . b12])
s12 = {b12/Sqrt[b12 . b12], 
        point12\[Cross]b12/Sqrt[b12 . b12]}
psi12 = 
      2 ArcTan[qS12[[1, 4]], Sqrt[b12 . b12]]
t12 = -2 qS12[[2, 4]]/Sqrt[b12 . b12]

b13 = {qS13[[1, 1]], qS13[[1, 2]], 
        qS13[[1, 3]]}
point13 = (b13/
            Sqrt[b13 . b13])\[Cross]({qS13[[2, 1]], qS13[[2, 2]], 
              qS13[[2, 3]]}/Sqrt[b13 . b13])
s13 = {b13/Sqrt[b13 . b13], 
        point13\[Cross]b13/Sqrt[b13 . b13]}
psi13 = 
      2 ArcTan[qS13[[1, 4]], Sqrt[b13 . b13]]
t13 = -2 qS13[[2, 4]]/Sqrt[b13 . b13]

b14 = {qS14[[1, 1]], qS14[[1, 2]], 
        qS14[[1, 3]]}
point14 = (b14/
            Sqrt[b14 . b14])\[Cross]({qS14[[2, 1]], qS14[[2, 2]], 
              qS14[[2, 3]]}/Sqrt[b14 . b14])
s14 = {b14/Sqrt[b14 . b14], 
        point14\[Cross]b14/Sqrt[b14 . b14]}
psi14 = 
      2 ArcTan[qS14[[1, 4]], Sqrt[b14 . b14]]
t14 = -2 qS14[[2, 4]]/Sqrt[b14 . b14]

b15 = {qS15[[1, 1]], qS15[[1, 2]], 
        qS15[[1, 3]]}
point15 = (b15/
            Sqrt[b15 . b15])\[Cross]({qS15[[2, 1]], qS15[[2, 2]], 
              qS15[[2, 3]]}/Sqrt[b15 . b15])
s15 = {b15/Sqrt[b15 . b15], 
        point15\[Cross]b15/Sqrt[b15 . b15]}
psi15 = 
      2 ArcTan[qS15[[1, 4]], Sqrt[b15 . b15]]
t15 = -2 qS15[[2, 4]]/Sqrt[b15 . b15]

b16 = {qS16[[1, 1]], qS16[[1, 2]], 
        qS16[[1, 3]]}
point16 = (b16/
            Sqrt[b16 . b16])\[Cross]({qS16[[2, 1]], qS16[[2, 2]], 
              qS16[[2, 3]]}/Sqrt[b16 . b16])
s16 = {b16/Sqrt[b16 . b16], 
        point16\[Cross]b16/Sqrt[b16 . b16]}
psi16 = 
      2 ArcTan[qS16[[1, 4]], Sqrt[b16 . b16]]
t16 = -2 qS16[[2, 4]]/Sqrt[b16 . b16]

b17 = {qS17[[1, 1]], qS17[[1, 2]], 
        qS17[[1, 3]]}
point17 = (b17/
            Sqrt[b17 . b17])\[Cross]({qS17[[2, 1]], qS17[[2, 2]], 
              qS17[[2, 3]]}/Sqrt[b17 . b17])
s17 = {b17/Sqrt[b17 . b17], 
        point17\[Cross]b17/Sqrt[b17 . b17]}
psi17 = 
      2 ArcTan[qS17[[1, 4]], Sqrt[b17 . b17]]
t17 = -2 qS17[[2, 4]]/Sqrt[b17 . b17]

b18 = {qS18[[1, 1]], qS18[[1, 2]], 
        qS18[[1, 3]]}
point18 = (b18/
            Sqrt[b18 . b18])\[Cross]({qS18[[2, 1]], qS18[[2, 2]], 
              qS18[[2, 3]]}/Sqrt[b18 . b18])
s18 = {b18/Sqrt[b18 . b18], 
        point18\[Cross]b18/Sqrt[b18 . b18]}
psi18 = 
      2 ArcTan[qS18[[1, 4]], Sqrt[b18 . b18]]
t18 = -2 qS18[[2, 4]]/Sqrt[b18 . b18]

b19 = {qS19[[1, 1]], qS19[[1, 2]], 
        qS19[[1, 3]]}
point19 = (b19/
            Sqrt[b19 . b19])\[Cross]({qS19[[2, 1]], qS19[[2, 2]], 
              qS19[[2, 3]]}/Sqrt[b19 . b19])
s19 = {b19/Sqrt[b19 . b19], 
        point19\[Cross]b19/Sqrt[b19 . b19]}
psi19 = 
      2 ArcTan[qS19[[1, 4]], Sqrt[b19 . b19]]
t19 = -2 qS19[[2, 4]]/Sqrt[b19 . b19]

quat12 = 
      N[{Flatten[{Sin[psi12/2]*s12[[1]], Cos[psi12/2]}], 
          Flatten[{Sin[psi12/2]*s12[[2]], -t12/2 Sin[psi12/2]}] + 
            Flatten[{t12/2 Cos[psi12/2]*s12[[1]], 0}]}]

quat13 = 
      N[{Flatten[{Sin[psi13/2]*s13[[1]], Cos[psi13/2]}], 
          Flatten[{Sin[psi13/2]*s13[[2]], -t13/2 Sin[psi13/2]}] + 
            Flatten[{t13/2 Cos[psi13/2]*s13[[1]], 0}]}]

quat14 = 
      N[{Flatten[{Sin[psi14/2]*s14[[1]], Cos[psi14/2]}], 
          Flatten[{Sin[psi14/2]*s14[[2]], -t14/2 Sin[psi14/2]}] + 
            Flatten[{t14/2 Cos[psi14/2]*s14[[1]], 0}]}]

quat15 = 
      N[{Flatten[{Sin[psi15/2]*s15[[1]], Cos[psi15/2]}], 
          Flatten[{Sin[psi15/2]*s15[[2]], -t15/2 Sin[psi15/2]}] + 
            Flatten[{t15/2 Cos[psi15/2]*s15[[1]], 0}]}]

quat16 = 
      N[{Flatten[{Sin[psi16/2]*s16[[1]], Cos[psi16/2]}], 
          Flatten[{Sin[psi16/2]*s16[[2]], -t16/2 Sin[psi16/2]}] + 
            Flatten[{t16/2 Cos[psi16/2]*s16[[1]], 0}]}]

qS16

quat17 = 
      N[{Flatten[{Sin[psi17/2]*s17[[1]], Cos[psi17/2]}], 
          Flatten[{Sin[psi17/2]*s17[[2]], -t17/2 Sin[psi17/2]}] + 
            Flatten[{t17/2 Cos[psi17/2]*s17[[1]], 0}]}]

qS17

quat18 = 
      N[{Flatten[{Sin[psi18/2]*s18[[1]], Cos[psi18/2]}], 
          Flatten[{Sin[psi18/2]*s18[[2]], -t18/2 Sin[psi18/2]}] + 
            Flatten[{t18/2 Cos[psi18/2]*s18[[1]], 0}]}]

qS18

quat19 = 
      N[{Flatten[{Sin[psi19/2]*s19[[1]], Cos[psi19/2]}], 
          Flatten[{Sin[psi19/2]*s19[[2]], -t19/2 Sin[psi19/2]}] + 
            Flatten[{t19/2 Cos[psi19/2]*s19[[1]], 0}]}]

qS19

so the rotation axis and angles are correct for the four position.

2. Transform to four homogeneous matrices

Pos1 = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 
          1}}

The relative rotations are found by using Cayley's formula

Id = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}

makeB[b_, psi_] := 
      Tan[psi/2]*{{0, -b[[3]], b[[2]]}, {b[[3]], 
            0, -b[[1]]}, {-b[[2]], b[[1]], 0}}

makeA[Id_, B_] := (Inverse[Id - B]) . (Id + B)

B12 = N[makeB[s12[[1]], psi12]]
A12 = makeA[Id, B12]
A2 = A12

The absolute rotations will be given by A12=A2A1^(-1), but as A1 is the identity, A12=A2 and so on.

B13 = N[makeB[s13[[1]], psi13]]
A13 = makeA[Id, B13]
A3 = A13

B14 = N[makeB[s14[[1]], psi14]]
A14 = makeA[Id, B14]
A4 = A14

B15 = N[makeB[s15[[1]], psi15]]
A15 = makeA[Id, B15]
A5 = A15

B16 = N[makeB[s16[[1]], psi16]]
A16 = makeA[Id, B16]
A6 = A16

B17 = N[makeB[s17[[1]], psi17]]
A17 = makeA[Id, B17]
A7 = A17

B18 = N[makeB[s18[[1]], psi18]]
A18 = makeA[Id, B18]
A8 = A18

B19 = N[makeB[s19[[1]], psi19]]
A19 = makeA[Id, B19]
A9 = A19

For the translation, let me try to use the following formula: (I-A)v=dxs. This will give the component perp. to s and then we need to add t.s.

d1 = {0, 0, 0}
dsol = {dx, dy, dz}

eq12 = (Id - A12) . s12[[2]] - 
        dsol\[Cross]s12[[1]]
sols = 
      Solve[{eq12[[1]] == 0, eq12[[2]] == 0, 
          dsol . s12[[1]] == 0}, {dx, dy, 
          dz}]
d12 = (dsol /. {sols[[1, 1]], sols[[1, 2]], 
              sols[[1, 3]]}) + 
        t12*s12[[1]]
d2 = d12 - A12 . d1
Pos2 = 
      Append[Transpose[Append[Transpose[A2], d2]], {0, 0, 0, 1}]

MatrixForm[Pos2]

eq13 = (Id - A13) . s13[[2]] - 
        dsol\[Cross]s13[[1]]
sols = 
      Solve[{eq13[[1]] == 0, eq13[[2]] == 0, 
          dsol . s13[[1]] == 0}, {dx, dy, 
          dz}]
d13 = (dsol /. {sols[[1, 1]], sols[[1, 2]], 
              sols[[1, 3]]}) + 
        t13*s13[[1]]
d3 = d13 - A13 . d1
Pos3 = 
      Append[Transpose[Append[Transpose[A3], d3]], {0, 0, 0, 1}]

eq14 = (Id - A14) . s14[[2]] - 
        dsol\[Cross]s14[[1]]
sols = 
      Solve[{eq14[[1]] == 0, eq14[[2]] == 0, 
          dsol . s14[[1]] == 0}, {dx, dy, 
          dz}]
d14 = (dsol /. {sols[[1, 1]], sols[[1, 2]], 
              sols[[1, 3]]}) + 
        t14*s14[[1]]
d4 = d14 - A14 . d1
Pos4 = 
      Append[Transpose[Append[Transpose[A4], d4]], {0, 0, 0, 1}]

eq15 = (Id - A15) . s15[[2]] - 
        dsol\[Cross]s15[[1]]
sols = 
      Solve[{eq15[[1]] == 0, eq15[[2]] == 0, 
          dsol . s15[[1]] == 0}, {dx, dy, 
          dz}]
d15 = (dsol /. {sols[[1, 1]], sols[[1, 2]], 
              sols[[1, 3]]}) + 
        t15*s15[[1]]
d5 = d15 - A15 . d1
Pos5 = 
      Append[Transpose[Append[Transpose[A5], d5]], {0, 0, 0, 1}]

eq16 = (Id - A16) . s16[[2]] - 
        dsol\[Cross]s16[[1]]
sols = 
      Solve[{eq16[[1]] == 0, eq16[[2]] == 0, 
          dsol . s16[[1]] == 0}, {dx, dy, 
          dz}]
d16 = (dsol /. {sols[[1, 1]], sols[[1, 2]], 
              sols[[1, 3]]}) + 
        t16*s16[[1]]
d6 = d16 - A16 . d1
Pos6 = 
      Append[Transpose[Append[Transpose[A6], d6]], {0, 0, 0, 1}]

MatrixForm[Pos6]

eq17 = (Id - A17) . s17[[2]] - 
        dsol\[Cross]s17[[1]]
sols = 
      Solve[{eq17[[1]] == 0, eq17[[2]] == 0, 
          dsol . s17[[1]] == 0}, {dx, dy, 
          dz}]
d17 = (dsol /. {sols[[1, 1]], sols[[1, 2]], 
              sols[[1, 3]]}) + 
        t17*s17[[1]]
d7 = d17 - A17 . d1
Pos7 = 
      Append[Transpose[Append[Transpose[A7], d7]], {0, 0, 0, 1}]

MatrixForm[Pos7]

eq18 = (Id - A18) . s18[[2]] - 
        dsol\[Cross]s18[[1]]
sols = 
      Solve[{eq18[[1]] == 0, eq18[[2]] == 0, 
          dsol . s18[[1]] == 0}, {dx, dy, 
          dz}]
d18 = (dsol /. {sols[[1, 1]], sols[[1, 2]], 
              sols[[1, 3]]}) + 
        t18*s18[[1]]
d8 = d18 - A18 . d1
Pos8 = 
      Append[Transpose[Append[Transpose[A8], d8]], {0, 0, 0, 1}]

MatrixForm[Pos8]

eq19 = (Id - A19) . s19[[2]] - 
        dsol\[Cross]s19[[1]]
sols = 
      Solve[{eq19[[1]] == 0, eq19[[2]] == 0, 
          dsol . s19[[1]] == 0}, {dx, dy, 
          dz}]
d19 = (dsol /. {sols[[1, 1]], sols[[1, 2]], 
              sols[[1, 3]]}) + 
        t19*s19[[1]]
d9 = d19 - A19 . d1
Pos9 = 
      Append[Transpose[Append[Transpose[A9], d9]], {0, 0, 0, 1}]

MatrixForm[Pos9]

check the inverse transform:

checkA = Pos4[[{1, 2, 3}, {1, 2, 3}]]

checkB = (checkA - Id) . Inverse[checkA + Id]

checkb = {checkB[[3, 2]], checkB[[1, 3]], 
        checkB[[2, 1]]}

N[Tan[psi14/2]*s14]

checks = checkb/Sqrt[checkb . checkb]

N[s14]

checkd = Flatten[Pos4[[{1, 2, 3}, {4}]] ]

checkv = 
      N[-1/(2*Tan[psi14/2]^3)*
          checkB . checkB . (Id - checkB) . checkd]

so the transform is fine!

3. Solutions solsRRRR: 
good ones are 2, 4, 5, 19, 21 from set 3
1, 3, 4, 7, 15, 21  from set last

Preliminaries

Ztrans[θ_, 
        d_] := {{Cos[θ], -Sin[θ], 0, 0}, {Sin[θ], 
          Cos[θ], 0, 0}, {0, 0, 1, d}, {0, 0, 0, 1}}

Xtrans[α_, 
        a_] := {{1, 0, 0, a}, {0, Cos[α], -Sin[α], 0}, {0, 
          Sin[α], Cos[α], 0}, {0, 0, 0, 1}}

ComputeNum[A_, hnum_, g1num_, g2num_, g3num_, pnum_] := 
      Simplify[A /. {g1x -> g1num[[1, 1]], 
            g1y -> g1num[[1, 2]], g1z -> g1num[[1, 3]], 
            g1x0 -> g1num[[2, 1]], g1y0 -> g1num[[2, 2]], 
            g1z0 -> g1num[[2, 3]], g2x -> g2num[[1, 1]], 
            g2y -> g2num[[1, 2]], g2z -> g2num[[1, 3]], 
            g2x0 -> g2num[[2, 1]], g2y0 -> g2num[[2, 2]], 
            g2z0 -> g2num[[2, 3]], g3x -> g3num[[1, 1]], 
            g3y -> g3num[[1, 2]], g3z -> g3num[[1, 3]], 
            g3x0 -> g3num[[2, 1]], g3y0 -> g3num[[2, 2]], 
            g3z0 -> g3num[[2, 3]], hx -> hnum[[1]], 
            hy -> hnum[[2]], hz -> hnum[[3]], 
            px -> pnum[[1, 1]], py -> pnum[[1, 2]], 
            pz -> pnum[[1, 3]], pw -> pnum[[1, 4]], 
            px0 -> pnum[[2, 1]], py0 -> pnum[[2, 2]], 
            pz0 -> pnum[[2, 3]], 
            pw0 -> pnum[[2, 4]]}]

dataRRRR3 = Import["/Users/alba/solsRRRR3", "Table"]

dataRRRR = 
      Import["/Users/alba/solsRRRRlast", "Table"]

The joint axes G,F,W,U expressed as lines

Length[dataRRRR]

nsol = 1

G = {dataRRRR[[nsol, {1, 2, 3}]], 
        dataRRRR[[nsol, {4, 5, 6}]]}

W = {dataRRRR[[nsol, {7, 8, 9}]], 
        dataRRRR[[nsol, {10, 11, 12}]]}

F = {dataRRRR[[nsol, {13, 14, 15}]], 
        dataRRRR[[nsol, {16, 17, 18}]]}

U = {dataRRRR[[nsol, {19, 20, 21}]], 
        dataRRRR[[nsol, {22, 23, 24}]]}

Make them unit and clean

Gline = {G[[1]]/Sqrt[G[[1]] . G[[1]]], 
        G[[2]]/Sqrt[G[[1]] . G[[1]]]}

Wline = {W[[1]]/Sqrt[W[[1]] . W[[1]]], 
        W[[2]]/Sqrt[W[[1]] . W[[1]]]}

Fline = {F[[1]]/Sqrt[F[[1]] . F[[1]]], 
        F[[2]]/Sqrt[F[[1]] . F[[1]]]}

Uline = {U[[1]]/Sqrt[U[[1]] . U[[1]]], 
        U[[2]]/Sqrt[U[[1]] . U[[1]]]}

Find transformation [G]

common normal to G and W

Ngw = ({Gline[[1]]\[Cross]Wline[[1]], 
            Gline[[1]]\[Cross]Wline[[2]] + 
              Gline[[2]]\[Cross]Wline[[1]]})/
        Sqrt[(Gline[[1]]\[Cross]
                Wline[[1]]) . (Gline[[1]]\[Cross]
                Wline[[1]])]

pointNgw = Ngw[[1]]\[Cross]Ngw[[2]]

Ngwline = {Ngw[[1]], pointNgw\[Cross]Ngw[[1]]}

Intersection point with G: Easier way to compute the intersection point: (pointN+t.Ndir)×Gdir=G0 (they must have same moment)

part = (Gline[[2, 1]] - pointNgw[[2]] Gline[[1, 3]] + 
            pointNgw[[3]] Gline[[1, 2]])/(Ngwline[[1, 
                  2]] Gline[[1, 3]] - 
            Ngwline[[1, 3]] Gline[[1, 2]])

pointInt = pointNgw + part Ngwline[[1]]

transformation [G] goes from origin to the frame {Ngp, G×Ngp,G} at the perpendicular point with G

yaxis = Gline[[1]]\[Cross]Ngw[[1]]

Rgtrans = {Ngw[[1]], yaxis, Gline[[1]], pointInt}

Ghom = Transpose[Rgtrans]
MatrixForm[Ghom]

Gtot = Append[Ghom, {0, 0, 0, 1}]

MatrixForm[Gtot]

Check Gtot:

rotmat = Gtot[[{1, 2, 3}, {1, 2, 3}]]

zeromat = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}

dispmat = Flatten[Gtot[[{1, 2, 3}, {4}]]]

dmat = makeB[dispmat, Pi/2]

bigmat = {Join[rotmat[[1]], {0, 0, 0}], 
        Join[rotmat[[2]], {0, 0, 0}], Join[rotmat[[3]], {0, 0, 0}], 
        Join[(dmat . rotmat)[[1]], rotmat[[1]]], 
        Join[(dmat . rotmat)[[2]], rotmat[[2]]], 
        Join[(dmat . rotmat)[[3]], rotmat[[3]]]}

MatrixForm[bigmat]

bigmat . {0, 0, 1, 0, 0, 0}

Gline

bigmat . {1, 0, 0, 0, 0, 0}

Ngwline

both G and P are as expected, so [G] is right.

Transformation from G to W along X

Length and angle between the first R-joint (G) and the P-joint (P)

cosanglegw = Gline[[1]] . Wline[[1]]

minusasinanglegw = 
      Gline[[1]] . Wline[[2]] + 
        Gline[[2]] . Wline[[1]]

angle1 = ArcCos[cosanglegw]

angle1*180/Pi

length1 = -minusasinanglegw/Sqrt[1 - cosanglegw^2]

If[(Drop[((Gtot . Xtrans[angle1, length1]) . {0, 0, 1, 
                    0}), {4}]) . Wline[[1]] == 1,  , 
      angle1 = -angle1; length1 = -length1]

check that we get the direction P, when the angle θ is zero:

(Gtot . Ztrans[0, 0] . Xtrans[angle1, length1]) . {0, 0, 1, 
        0}

Wline[[1]]

check that we get the point:

(Gtot . Ztrans[0, 0] . Xtrans[angle1, length1]) . {0, 0, 0, 
        1}

part = (Wline[[2, 1]] - pointNgw[[2]] Wline[[1, 3]] + 
            pointNgw[[3]] Wline[[1, 2]])/(Ngwline[[1, 
                  2]] Wline[[1, 3]] - 
            Ngwline[[1, 3]] Wline[[1, 2]])

pointInt2 = pointNgw + part Ngwline[[1]]

The parameters (α1,a1) are right.

Transformation from W origin to common normal with F (along Z)

Length and angle along the P joint to get to the common normal between P and W:

Nwf = {Wline[[1]]\[Cross]Fline[[1]], 
        Wline[[1]]\[Cross]Fline[[2]] + 
          Wline[[2]]\[Cross]Fline[[1]]}

Nwfunit = Nwf/Sqrt[Nwf[[1]] . Nwf[[1]]]

pointNwf = Nwfunit[[1]]\[Cross]Nwfunit[[2]]

Nwfline = {Nwfunit[[1]], 
        pointNwf\[Cross]Nwfunit[[1]]}

The way this is computed, the twist angle between lines will be always between 0 and 180. Otherwise.....

cosangle2 = Ngwline[[1]] . Nwfline[[1]]

minusasinangle2 = 
      Ngwline[[1]] . Nwfline[[2]] + 
        Ngwline[[2]] . Nwfline[[1]]

attention to this minus sign

angle2 = ArcCos[cosangle2]

angle2*180/Pi

length2 = (-minusasinangle2/Sqrt[1 - cosangle2^2])

If[(Drop[((Gtot . Xtrans[angle1, length1] . 
                      Ztrans[angle2, length2]) . {1, 0, 0, 0}), {4}]) . 
          Nwfline[[1]] == 1,  , angle2 = -angle2; 
      length2 = -length2]

check that this transformation is right:

(Gtot . Ztrans[0, 0] . Xtrans[angle1, length1] . Ztrans[0, 0] . 
          Ztrans[angle2, length2]) . {0, 0, 1, 0}

Wline[[1]]

(Gtot . Ztrans[0, 0] . Xtrans[angle1, length1] . Ztrans[0, 0] . 
          Ztrans[angle2, length2]) . {1, 0, 0, 0}

Nwfline[[1]]

checking the moment:

rotmat = (Gtot . Ztrans[0, 0] . Xtrans[angle1, length1] . 
            Ztrans[0, 0] . Ztrans[angle2, length2])[[{1, 2, 3}, {1, 2, 
            3}]]

zeromat = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}

dispmat = 
      Flatten[(Gtot . Ztrans[0, 0] . Xtrans[angle1, length1] . 
              Ztrans[0, 0] . Ztrans[angle2, length2])[[{1, 2, 
              3}, {4}]]]

dmat = makeB[dispmat, Pi/2]

bigmat = {Join[rotmat[[1]], {0, 0, 0}], 
        Join[rotmat[[2]], {0, 0, 0}], Join[rotmat[[3]], {0, 0, 0}], 
        Join[(dmat . rotmat)[[1]], rotmat[[1]]], 
        Join[(dmat . rotmat)[[2]], rotmat[[2]]], 
        Join[(dmat . rotmat)[[3]], rotmat[[3]]]}

MatrixForm[bigmat]

bigmat . {0, 0, 1, 0, 0, 0}

Wline

bigmat . {1, 0, 0, 0, 0, 0}

Nwfline

Transformation between W and F along X

Length and angle between P-joint and second R-joint (W):
need to find common normal line.

cosanglewf = Wline[[1]] . Fline[[1]]

minusasinanglewf = 
      Wline[[1]] . Fline[[2]] + 
        Wline[[2]] . Fline[[1]]

angle3 = ArcCos[cosanglewf]

angle3*180/Pi

length3 = -minusasinanglewf/Sqrt[1 - cosanglewf^2]

If[(Drop[((Gtot . Xtrans[angle1, length1] . 
                      Ztrans[angle2, length2] . 
                      Xtrans[angle3, length3]) . {0, 0, 1, 0}), {4}]) . 
          Fline[[1]] == 1,  , angle3 = -angle3; 
      length3 = -length3]

Check that this arrives to W:

(Gtot . Ztrans[0, 0] . Xtrans[angle1, length1] . Ztrans[0, 0] . 
          Ztrans[angle2, length2] . Xtrans[angle3, length3]) . {0, 0, 1, 
        0}

Fline[[1]]

Transformation from F origin to common normal with U (along Z)

Length and angle along the P joint to get to the common normal between P and W:

Nfu = {Fline[[1]]\[Cross]Uline[[1]], 
        Fline[[1]]\[Cross]Uline[[2]] + 
          Fline[[2]]\[Cross]Uline[[1]]}

Nfuunit = Nfu/Sqrt[Nfu[[1]] . Nfu[[1]]]

pointNfu = Nfuunit[[1]]\[Cross]Nfuunit[[2]]

Nfuline = {Nfuunit[[1]], 
        pointNfu\[Cross]Nfuunit[[1]]}

The way this is computed, the twist angle between lines will be always between 0 and 180. Otherwise.....

cosangle4 = Nwfline[[1]] . Nfuline[[1]]

minusasinangle4 = 
      Nwfline[[1]] . Nfuline[[2]] + 
        Nwfline[[2]] . Nfuline[[1]]

attention to this minus sign

angle4 = ArcCos[cosangle4]

angle4*180/Pi

length4 = (-minusasinangle4/Sqrt[1 - cosangle4^2])

If[(Drop[((Gtot . Xtrans[angle1, length1] . 
                      Ztrans[angle2, length2] . Xtrans[angle3, length3] . 
                      Ztrans[angle4, length4]) . {1, 0, 0, 0}), {4}]) . 
          Nfuline[[1]] == 1,  , angle4 = -angle4; 
      length4 = -length4]

check that this transformation is right:

(Gtot . Ztrans[0, 0] . Xtrans[angle1, length1] . Ztrans[0, 0] . 
          Ztrans[angle2, length2] . Xtrans[angle3, length3] . Ztrans[0, 0] . 
          Ztrans[angle4, length4]) . {0, 0, 1, 0}

Fline[[1]]

(Gtot . Ztrans[0, 0] . Xtrans[angle1, length1] . Ztrans[0, 0] . 
          Ztrans[angle2, length2] . Xtrans[angle3, length3] . 
          Ztrans[angle4, length4]) . {1, 0, 0, 0}

Nfuline[[1]]

checking the moment:

rotmat = (Gtot . Ztrans[0, 0] . Xtrans[angle1, length1] . 
            Ztrans[0, 0] . Ztrans[angle2, length2] . 
            Xtrans[angle3, length3] . Ztrans[0, 0] . 
            Ztrans[angle4, length4])[[{1, 2, 3}, {1, 2, 3}]]

zeromat = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}

dispmat = 
      Flatten[(Gtot . Ztrans[0, 0] . Xtrans[angle1, length1] . 
              Ztrans[0, 0] . Ztrans[angle2, length2] . 
              Xtrans[angle3, length3] . Ztrans[0, 0] . 
              Ztrans[angle4, length4])[[{1, 2, 3}, {4}]]]

dmat = makeB[dispmat, Pi/2]

bigmat = {Join[rotmat[[1]], {0, 0, 0}], 
        Join[rotmat[[2]], {0, 0, 0}], Join[rotmat[[3]], {0, 0, 0}], 
        Join[(dmat . rotmat)[[1]], rotmat[[1]]], 
        Join[(dmat . rotmat)[[2]], rotmat[[2]]], 
        Join[(dmat . rotmat)[[3]], rotmat[[3]]]}

MatrixForm[bigmat]

bigmat . {0, 0, 1, 0, 0, 0}

Fline

bigmat . {1, 0, 0, 0, 0, 0}

Nfuline

Transformation between F and U along X

Length and angle between P-joint and second R-joint (W):
need to find common normal line.

cosanglefu = Fline[[1]] . Uline[[1]]

minusasinanglefu = 
      Fline[[1]] . Uline[[2]] + 
        Fline[[2]] . Uline[[1]]

angle5 = ArcCos[cosanglefu]

angle5*180/Pi

length5 = -minusasinanglefu/Sqrt[1 - cosanglefu^2]

If[(Drop[((Gtot . Ztrans[0, 0] . Xtrans[angle1, length1] . 
                      Ztrans[0, 0] . Ztrans[angle2, length2] . 
                      Xtrans[angle3, length3] . Ztrans[0, 0] . 
                      Ztrans[angle4, length4] . 
                      Xtrans[angle5, length5]) . {0, 0, 1, 0}), {4}]) . 
          Uline[[1]] == 1,  , angle5 = -angle5; 
      length5 = -length5]

Check that this arrives to W:

(Gtot . Ztrans[0, 0] . Xtrans[angle1, length1] . Ztrans[0, 0] . 
          Ztrans[angle2, length2] . Xtrans[angle3, length3] . Ztrans[0, 0] . 
          Ztrans[angle4, length4] . Xtrans[angle5, length5] . 
          Ztrans[0, 0]) . {0, 0, 1, 0}

Uline[[1]]

Relative transformation [H] from W to the gripper.

Relative transformation from last R-joint to the gripper:  calculated as position 1. We need the angles theta and phi to go there to the first position, which are zero.....
Just in case, we already have the procedure to find them.

thetas = 
      dataRRRR[[nsol, {25, 26, 27, 28, 29, 30, 31, 32}]]

phis = 
      dataRRRR[[nsol, {33, 34, 35, 36, 37, 38, 39, 40}]]

psis = 
      dataRRRR[[nsol, {41, 42, 43, 44, 45, 46, 47, 48}]]

xis = dataRRRR[[nsol, {49, 50, 51, 52, 53, 54, 55, 56}]]

this is the solution for the angles to reach the second position wrt the angles for the first,

theta10 = 0
phi10 = 0
psi10 = 0
xi10 = 0

now we compose the kinematic equations

KE10 = 
      Gtot . Ztrans[0, 0] . Xtrans[angle1, length1] . Ztrans[0, 0] . 
        Ztrans[angle2, length2] . Xtrans[angle3, length3] . Ztrans[0, 0] . 
        Ztrans[angle4, length4] . Xtrans[angle5, length5] . 
        Ztrans[0, 0]

The homogeneous transform for the first position P1: it is a rotation about the z-axis only,

Htot = Inverse[KE10] . Pos1

KE10 . Htot

Check if the chain hits the rest of positions

Position 2: Find values of the joint variables

Assemble kinematic equations and compare with absolute transformation

KE12 = 
      Gtot . Ztrans[thetas[[1]], 0] . Xtrans[angle1, length1] . 
        Ztrans[phis[[1]], 0] . Ztrans[angle2, length2] . 
        Xtrans[angle3, length3] . Ztrans[psis[[1]], 0] . 
        Ztrans[angle4, length4] . Xtrans[angle5, length5] . 
        Ztrans[xis[[1]], 0] . Htot

Pos2

Now the reaching of the second position is right!!! :))))

Position 3: Find values of the joint variables

KE13 = 
      Gtot . Ztrans[thetas[[2]], 0] . Xtrans[angle1, length1] . 
        Ztrans[phis[[2]], 0] . Ztrans[angle2, length2] . 
        Xtrans[angle3, length3] . Ztrans[psis[[2]], 0] . 
        Ztrans[angle4, length4] . Xtrans[angle5, length5] . 
        Ztrans[xis[[2]], 0] . Htot

Pos3

?????

Position 4: Find values of the joint variables

KE14 = 
      Gtot . Ztrans[thetas[[3]], 0] . Xtrans[angle1, length1] . 
        Ztrans[phis[[3]], 0] . Ztrans[angle2, length2] . 
        Xtrans[angle3, length3] . Ztrans[psis[[3]], 0] . 
        Ztrans[angle4, length4] . Xtrans[angle5, length5] . 
        Ztrans[xis[[3]], 0] . Htot

Pos4

Position 5: Find values of the joint variables

KE15 = 
      Gtot . Ztrans[thetas[[4]], 0] . Xtrans[angle1, length1] . 
        Ztrans[phis[[4]], 0] . Ztrans[angle2, length2] . 
        Xtrans[angle3, length3] . Ztrans[psis[[4]], 0] . 
        Ztrans[angle4, length4] . Xtrans[angle5, length5] . 
        Ztrans[xis[[4]], 0] . Htot

Pos5

Position 6: Find values of the joint variables

KE16 = 
      Gtot . Ztrans[thetas[[5]], 0] . Xtrans[angle1, length1] . 
        Ztrans[phis[[5]], 0] . Ztrans[angle2, length2] . 
        Xtrans[angle3, length3] . Ztrans[psis[[5]], 0] . 
        Ztrans[angle4, length4] . Xtrans[angle5, length5] . 
        Ztrans[xis[[5]], 0] . Htot

Pos6

Position 7: Find values of the joint variables

KE17 = 
      Gtot . Ztrans[thetas[[6]], 0] . Xtrans[angle1, length1] . 
        Ztrans[phis[[6]], 0] . Ztrans[angle2, length2] . 
        Xtrans[angle3, length3] . Ztrans[psis[[6]], 0] . 
        Ztrans[angle4, length4] . Xtrans[angle5, length5] . 
        Ztrans[xis[[6]], 0] . Htot

Pos7

Position 8: Find values of the joint variables

KE18 = 
      Gtot . Ztrans[thetas[[7]], 0] . Xtrans[angle1, length1] . 
        Ztrans[phis[[7]], 0] . Ztrans[angle2, length2] . 
        Xtrans[angle3, length3] . Ztrans[psis[[7]], 0] . 
        Ztrans[angle4, length4] . Xtrans[angle5, length5] . 
        Ztrans[xis[[7]], 0] . Htot

Pos8

Position 9: Find values of the joint variables

KE19 = 
      Gtot . Ztrans[thetas[[8]], 0] . Xtrans[angle1, length1] . 
        Ztrans[phis[[8]], 0] . Ztrans[angle2, length2] . 
        Xtrans[angle3, length3] . Ztrans[psis[[8]], 0] . 
        Ztrans[angle4, length4] . Xtrans[angle5, length5] . 
        Ztrans[xis[[8]], 0] . Htot

Pos9

Chop[KE12 - Pos2]
Chop[KE13 - Pos3]
Chop[KE14 - Pos4]
Chop[KE15 - Pos5]
Chop[KE16 - Pos6]
Chop[KE17 - Pos7]
Chop[KE18 - Pos8]
Chop[KE19 - Pos9]

OK!

final check

Chop[KE12 - Pos2]
Chop[KE13 - Pos3]
Chop[KE14 - Pos4]
Chop[KE15 - Pos5]
Chop[KE16 - Pos6]
Chop[KE17 - Pos7]
Chop[KE18 - Pos8]
Chop[KE19 - Pos9]

dataAxes3 = 
        Table[dataRRRR3[[i, Range[1, 24]]], {i, 1, 
            Length[dataRRRR3]}];

dataAxes = 
        Table[dataRRRR[[i, Range[1, 24]]], {i, 1, 
            Length[dataRRRR]}];

RRRRgood = 
      Flatten[{dataAxes3[[{2, 4, 5, 19, 21}]], 
          dataAxes[[{1, 3, 4, 7, 15, 21}]]}, 1]

RRRRgood3 = SetPrecision[RRRRgood, 3]

Export["RRRRthesis", RRRRgood3, "Lines"]

dataRRRR[[2]]

Data for Synthetica

Gtot

Htot

angle1*180/Pi
length1

angle2*180/Pi
length2

angle3*180/Pi
length3

angle4*180/Pi
length4

angle5*180/Pi
length5

thetaS = 
      Prepend[thetas, 0]*180/Pi + {0, 0, -2*Pi, 0, 0, 0, 0, 0, 0}*180/
            Pi

phiS = (Prepend[phis, 0] + {angle2, angle2, angle2, angle2, angle2, 
                angle2, angle2, angle2, angle2})*180/Pi + {0, 0, 0, 0, 0, 
            0, 0, 0, 0}*180/Pi

psiS = (Prepend[psis, 0] + {angle4, angle4, angle4, angle4, angle4, 
                angle4, angle4, angle4, angle4})*180/Pi + {0, 0, 0, 0, 0, 
            0, 0, -2*Pi, 0}*180/Pi

xiS = (Prepend[xis, 0])*180/Pi + {0, 0, 0, 0, 2*Pi, 0, 0, 0, 0}*180/
            Pi

4. Checking the analysis

The dual quaternion

quat4r = {Cos[(θ/2)] (Cos[(ϕ/2)] (ux Cos[(ψ/2)] Sin[(ξ/2)] + (fx Cos[(ξ/2)] + (-fz uy + 
                                  fy uz) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] (Cos[(ψ/2)] (wx Cos[(ξ/2)] + (uz wy - uy wz) Sin[(ξ/2)]) + ((fz wy - 
                                  fy wz) Cos[(ξ/2)] + (-fx ux wx - fy uy wx - fz uz wx - fy ux wy + fx uy wy - fz ux wz + 
                                  fx uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
          Sin[(θ/2)] (Cos[(ϕ/2)] (Cos[(ψ/2)] (gx Cos[(ξ/2)] + (-gz uy + gy uz) Sin[(ξ/2)]) + ((fz gy - 
                                  fy gz) Cos[(ξ/2)] + (-fx gx ux - fy gy ux - fz gz ux - fy gx uy + fx gy uy - fz gx uz + 
                                  fx gz uz) Sin[(ξ/2)]) Sin[(ψ/2)]) + Sin[(ϕ/2)] (Cos[(ψ/2)] ((-gz wy + 
                                  gy wz) Cos[(ξ/2)] + (-gx ux wx + gy uy wx + gz uz wx - gy ux wy - gx uy wy - gz ux wz - 
                                  gx uz wz) Sin[(ξ/2)]) + ((-fx gx wx + fy gy wx + fz gz wx - fy gx wy - fx gy wy - 
                                  fz gx wz - 
                                  fx gz wz) Cos[(ξ/2)] + (fz gy ux wx - fy gz ux wx + fz gx uy wx + fx gz uy wx - fy gx uz wx - 
                                  fx gy uz wx - fz gx ux wy + 
                                  fx gz ux wy + fz gy uy wy + 
                                  fy gz uy wy + fx gx uz wy - 
                                  fy gy uz wy + fz gz uz wy + 
                                  fy gx ux wz - fx gy ux wz - 
                                  fx gx uy wz - fy gy uy wz + 
                                  fz gz uy wz - fz gy uz wz - 
                                  fy gz uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])), 
        Cos[(θ/2)] (Cos[(ϕ/2)] (uy Cos[(ψ/2)] Sin[(ξ/2)] + (fy Cos[(ξ/2)] + (fz ux - 
                                  fx uz) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] (Cos[(ψ/2)] (wy Cos[(ξ/2)] + (-uz wx + ux wz) Sin[(ξ/2)]) + ((-fz wx + 
                                  fx wz) Cos[(ξ/2)] + (fy ux wx - 
                                  fx uy wx - fx ux wy - fy uy wy - 
                                  fz uz wy - fz uy wz + 
                                  
                                  fy uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
          Sin[(θ/2)] (Cos[(ϕ/2)] (Cos[(ψ/2)] (gy Cos[(ξ/2)] + (gz ux - gx uz) Sin[(ξ/2)]) + ((-fz gx + 
                                  fx gz) Cos[(ξ/2)] + (fy gx ux - 
                                  fx gy ux - fx gx uy - fy gy uy - 
                                  fz gz uy - fz gy uz + 
                                  fy gz uz) Sin[(ξ/2)]) Sin[(ψ/2)]) + Sin[(ϕ/2)] (Cos[(ψ/2)] ((gz wx - 
                                  gx wz) Cos[(ξ/2)] + (-gy ux wx - gx uy wx + gx ux wy - gy uy wy + gz uz wy - gz uy wz - 
                                  gy uz wz) Sin[(ξ/2)]) + ((-fy gx wx - fx gy wx + fx gx wy - fy gy wy + fz gz wy - 
                                  fz gy wz - 
                                  fy gz wz) Cos[(ξ/2)] + (-fz gx ux wx - fx gz ux wx + fz gy uy wx - fy gz uy wx + 
                                  fx gx uz wx - fy gy uz wx - 
                                  fz gz uz wx - fz gy ux wy - 
                                  fy gz ux wy - fz gx uy wy + 
                                  fx gz uy wy + fy gx uz wy + 
                                  fx gy uz wy + fx gx ux wz + 
                                  fy gy ux wz - fz gz ux wz + 
                                  fy gx uy wz - fx gy uy wz + 
                                  fz gx uz wz + 
                                  fx gz uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])), 
        Cos[(θ/2)] (Cos[(ϕ/2)] (uz Cos[(ψ/2)] Sin[(ξ/2)] + (fz Cos[(ξ/2)] + (-fy ux + 
                                  fx uy) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] (Cos[(ψ/2)] (wz Cos[(ξ/2)] + (uy wx - ux wy) Sin[(ξ/2)]) + ((fy wx - 
                                  fx wy) Cos[(ξ/2)] + (fz ux wx - 
                                  fx uz wx + fz uy wy - fy uz wy - 
                                  fx ux wz - fy uy wz - 
                                  fz uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
          Sin[(θ/2)] (Cos[(ϕ/2)] (Cos[(ψ/2)] (gz Cos[(ξ/2)] + (-gy ux + gx uy) Sin[(ξ/2)]) + ((fy gx - 
                                  fx gy) Cos[(ξ/2)] + (fz gx ux - 
                                  fx gz ux + fz gy uy - fy gz uy - 
                                  fx gx uz - fy gy uz - 
                                  fz gz uz) Sin[(ξ/2)]) Sin[(ψ/2)]) + Sin[(ϕ/2)] (Cos[(ψ/2)] ((-gy wx + 
                                  
                                  gx wy) Cos[(ξ/2)] + (-gz ux wx - gx uz wx - gz uy wy - gy uz wy + gx ux wz + gy uy wz - 
                                  gz uz wz) Sin[(ξ/2)]) + ((-fz gx wx - fx gz wx - fz gy wy - fy gz wy + fx gx wz + 
                                  fy gy wz - 
                                  fz gz wz) Cos[(ξ/2)] + (fy gx ux wx + fx gy ux wx - fx gx uy wx + fy gy uy wx + fz gz uy wx + 
                                  fz gy uz wx - fy gz uz wx - 
                                  fx gx ux wy + fy gy ux wy - 
                                  fz gz ux wy - fy gx uy wy - 
                                  fx gy uy wy - fz gx uz wy + 
                                  fx gz uz wy + fz gy ux wz + 
                                  fy gz ux wz - fz gx uy wz - 
                                  fx gz uy wz + fy gx uz wz - 
                                  fx gy uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])), 
        Cos[(θ/2)] (Cos[(ϕ/2)] (Cos[(ξ/2)] Cos[(ψ/2)] + (-fx ux - fy uy - fz uz) Sin[(ξ/2)] Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] ((-ux wx - uy wy - 
                            uz wz) Cos[(ψ/2)] Sin[(ξ/2)] + ((-fx wx - fy wy - fz wz) Cos[(ξ/2)] + (fz uy wx - 
                                  fy uz wx - fz ux wy + fx uz wy + 
                                  fy ux wz - 
                                  fx uy wz) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
          Sin[(θ/2)] (Cos[(ϕ/2)] ((-gx ux - gy uy - 
                            gz uz) Cos[(ψ/2)] Sin[(ξ/2)] + ((-fx gx - fy gy - fz gz) Cos[(ξ/2)] + (-fz gy ux + 
                                  fy gz ux + fz gx uy - fx gz uy - 
                                  fy gx uz + 
                                  fx gy uz) Sin[(ξ/2)]) Sin[(ψ/2)]) + Sin[(ϕ/2)] (Cos[(ψ/2)] ((-gx wx - gy wy - 
                                  gz wz) Cos[(ξ/2)] + (-gz uy wx + gy uz wx + gz ux wy - gx uz wy - gy ux wz + 
                                  gx uy wz) Sin[(ξ/2)]) + ((fz gy wx - fy gz wx - fz gx wy + fx gz wy + fy gx wz - 
                                  fx gy wz) Cos[(ξ/2)] + (fx gx ux wx - fy gy ux wx - fz gz ux wx + fy gx uy wx + fx gy uy wx + 
                                  fz gx uz wx + fx gz uz wx + 
                                  fy gx ux wy + fx gy ux wy - 
                                  fx gx uy wy + fy gy uy wy - 
                                  fz gz uy wy + fz gy uz wy + 
                                  fy gz uz wy + fz gx ux wz + 
                                  fx gz ux wz + fz gy uy wz + 
                                  fy gz uy wz - fx gx uz wz - 
                                  fy gy uz wz + 
                                  fz gz uz wz) Sin[(ξ/2)]) Sin[(ψ/2)])), 
        Cos[(θ/2)] (Cos[(ϕ/2)] (ux0 Cos[(ψ/2)] Sin[(ξ/2)] + (fx0 Cos[(ξ/2)] + (-fz0 uy - fz uy0 + fy0 uz + 
                                  fy uz0) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] (Cos[(ψ/2)] (wx0 Cos[(ξ/2)] + (uz0 wy + uz wy0 - uy0 wz - 
                                  uy wz0) Sin[(ξ/2)]) + ((fz0 wy + fz wy0 - fy0 wz - fy wz0) Cos[(ξ/2)] + (-fx0 ux wx - 
                                  fx ux0 wx - fy0 uy wx - fy uy0 wx - 
                                  fz0 uz wx - fz uz0 wx - fx ux wx0 - 
                                  fy uy wx0 - fz uz wx0 - fy0 ux wy - 
                                  fy ux0 wy + fx0 uy wy + fx uy0 wy - 
                                  fy ux wy0 + fx uy wy0 - fz0 ux wz - 
                                  fz ux0 wz + fx0 uz wz + fx uz0 wz - 
                                  fz ux wz0 + 
                                  fx uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
          Sin[(θ/2)] (Cos[(ϕ/2)] (Cos[(ψ/2)] (gx0 Cos[(ξ/2)] + (-gz0 uy - gz uy0 + gy0 uz + 
                                  gy uz0) Sin[(ξ/2)]) + ((fz0 gy + fz gy0 - fy0 gz - fy gz0) Cos[(ξ/2)] + (-fx0 gx ux - 
                                  fx gx0 ux - fy0 gy ux - fy gy0 ux - 
                                  fz0 gz ux - fz gz0 ux - fx gx ux0 - 
                                  fy gy ux0 - fz gz ux0 - fy0 gx uy - 
                                  fy gx0 uy + fx0 gy uy + fx gy0 uy - 
                                  fy gx uy0 + fx gy uy0 - fz0 gx uz - 
                                  fz gx0 uz + fx0 gz uz + fx gz0 uz - 
                                  fz gx uz0 + 
                                  fx gz uz0) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] (Cos[(ψ/2)] ((-gz0 wy - 
                                  gz wy0 + gy0 wz + 
                                  gy wz0) Cos[(ξ/2)] + (-gx0 ux wx - gx ux0 wx + gy0 uy wx + gy uy0 wx + gz0 uz wx + gz uz0 wx - 
                                  gx ux wx0 + gy uy wx0 + gz uz wx0 - 
                                  gy0 ux wy - gy ux0 wy - gx0 uy wy - 
                                  gx uy0 wy - gy ux wy0 - gx uy wy0 - 
                                  gz0 ux wz - gz ux0 wz - gx0 uz wz - 
                                  gx uz0 wz - gz ux wz0 - 
                                  gx uz wz0) Sin[(ξ/2)]) + ((-fx0 gx wx - fx gx0 wx + fy0 gy wx + fy gy0 wx + 
                                  fz0 gz wx + fz gz0 wx - fx gx wx0 + 
                                  fy gy wx0 + fz gz wx0 - fy0 gx wy - 
                                  fy gx0 wy - fx0 gy wy - fx gy0 wy - 
                                  fy gx wy0 - fx gy wy0 - fz0 gx wz - 
                                  fz gx0 wz - fx0 gz wz - fx gz0 wz - 
                                  fz gx wz0 - 
                                  fx gz wz0) Cos[(ξ/2)] + (fz0 gy ux wx + fz gy0 ux wx - fy0 gz ux wx - fy gz0 ux wx + 
                                  fz gy ux0 wx - fy gz ux0 wx + 
                                  fz0 gx uy wx + fz gx0 uy wx + 
                                  fx0 gz uy wx + fx gz0 uy wx + 
                                  fz gx uy0 wx + fx gz uy0 wx - 
                                  fy0 gx uz wx - fy gx0 uz wx - 
                                  fx0 gy uz wx - fx gy0 uz wx - 
                                  fy gx uz0 wx - fx gy uz0 wx + 
                                  fz gy ux wx0 - fy gz ux wx0 + 
                                  fz gx uy wx0 + fx gz uy wx0 - 
                                  fy gx uz wx0 - fx gy uz wx0 - 
                                  fz0 gx ux wy - fz gx0 ux wy + 
                                  fx0 gz ux wy + fx gz0 ux wy - 
                                  fz gx ux0 wy + fx gz ux0 wy + 
                                  fz0 gy uy wy + fz gy0 uy wy + 
                                  fy0 gz uy wy + fy gz0 uy wy + 
                                  fz gy uy0 wy + fy gz uy0 wy + 
                                  fx0 gx uz wy + fx gx0 uz wy - 
                                  fy0 gy uz wy - fy gy0 uz wy + 
                                  fz0 gz uz wy + fz gz0 uz wy + 
                                  fx gx uz0 wy - fy gy uz0 wy + 
                                  fz gz uz0 wy - fz gx ux wy0 + 
                                  fx gz ux wy0 + fz gy uy wy0 + 
                                  fy gz uy wy0 + fx gx uz wy0 - 
                                  fy gy uz wy0 + fz gz uz wy0 + 
                                  fy0 gx ux wz + fy gx0 ux wz - 
                                  fx0 gy ux wz - fx gy0 ux wz + 
                                  fy gx ux0 wz - fx gy ux0 wz - 
                                  fx0 gx uy wz - fx gx0 uy wz - 
                                  fy0 gy uy wz - fy gy0 uy wz + 
                                  fz0 gz uy wz + fz gz0 uy wz - 
                                  fx gx uy0 wz - fy gy uy0 wz + 
                                  fz gz uy0 wz - fz0 gy uz wz - 
                                  fz gy0 uz wz - fy0 gz uz wz - 
                                  fy gz0 uz wz - fz gy uz0 wz - 
                                  fy gz uz0 wz + fy gx ux wz0 - 
                                  fx gy ux wz0 - fx gx uy wz0 - 
                                  fy gy uy wz0 + fz gz uy wz0 - 
                                  fz gy uz wz0 - 
                                  fy gz uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)])), 
        Cos[(θ/2)] (Cos[(ϕ/2)] (uy0 Cos[(ψ/2)] Sin[(ξ/2)] + (fy0 Cos[(ξ/2)] + (fz0 ux + fz ux0 - fx0 uz - 
                                  fx uz0) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] (Cos[(ψ/2)] (wy0 Cos[(ξ/2)] + (-uz0 wx - uz wx0 + ux0 wz + 
                                  ux wz0) Sin[(ξ/2)]) + ((-fz0 wx - fz wx0 + fx0 wz + fx wz0) Cos[(ξ/2)] + (fy0 ux wx + 
                                  fy ux0 wx - fx0 uy wx - fx uy0 wx + 
                                  fy ux wx0 - fx uy wx0 - fx0 ux wy - 
                                  fx ux0 wy - fy0 uy wy - fy uy0 wy - 
                                  fz0 uz wy - fz uz0 wy - fx ux wy0 - 
                                  fy uy wy0 - fz uz wy0 - fz0 uy wz - 
                                  fz uy0 wz + fy0 uz wz + fy uz0 wz - 
                                  fz uy wz0 + 
                                  fy uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
          Sin[(θ/2)] (Cos[(ϕ/2)] (Cos[(ψ/2)] (gy0 Cos[(ξ/2)] + (gz0 ux + gz ux0 - gx0 uz - 
                                  gx uz0) Sin[(ξ/2)]) + ((-fz0 gx - fz gx0 + fx0 gz + fx gz0) Cos[(ξ/2)] + (fy0 gx ux + 
                                  fy gx0 ux - fx0 gy ux - fx gy0 ux + 
                                  fy gx ux0 - fx gy ux0 - fx0 gx uy - 
                                  fx gx0 uy - fy0 gy uy - fy gy0 uy - 
                                  fz0 gz uy - fz gz0 uy - fx gx uy0 - 
                                  fy gy uy0 - fz gz uy0 - fz0 gy uz - 
                                  fz gy0 uz + fy0 gz uz + fy gz0 uz - 
                                  fz gy uz0 + 
                                  fy gz uz0) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] (Cos[(ψ/2)] ((gz0 wx + gz wx0 - 
                                  gx0 wz - 
                                  gx wz0) Cos[(ξ/2)] + (-gy0 ux wx - gy ux0 wx - gx0 uy wx - gx uy0 wx - gy ux wx0 - gx uy wx0 + 
                                  gx0 ux wy + gx ux0 wy - gy0 uy wy - 
                                  gy uy0 wy + gz0 uz wy + gz uz0 wy + 
                                  gx ux wy0 - gy uy wy0 + gz uz wy0 - 
                                  gz0 uy wz - gz uy0 wz - gy0 uz wz - 
                                  gy uz0 wz - gz uy wz0 - 
                                  gy uz wz0) Sin[(ξ/2)]) + ((-fy0 gx wx - fy gx0 wx - fx0 gy wx - fx gy0 wx - 
                                  fy gx wx0 - fx gy wx0 + fx0 gx wy + 
                                  fx gx0 wy - fy0 gy wy - fy gy0 wy + 
                                  fz0 gz wy + fz gz0 wy + fx gx wy0 - 
                                  fy gy wy0 + fz gz wy0 - fz0 gy wz - 
                                  fz gy0 wz - fy0 gz wz - fy gz0 wz - 
                                  fz gy wz0 - 
                                  fy gz wz0) Cos[(ξ/2)] + (-fz0 gx ux wx - fz gx0 ux wx - fx0 gz ux wx - fx gz0 ux wx - 
                                  fz gx ux0 wx - fx gz ux0 wx + 
                                  fz0 gy uy wx + fz gy0 uy wx - 
                                  fy0 gz uy wx - fy gz0 uy wx + 
                                  fz gy uy0 wx - fy gz uy0 wx + 
                                  fx0 gx uz wx + fx gx0 uz wx - 
                                  fy0 gy uz wx - fy gy0 uz wx - 
                                  fz0 gz uz wx - fz gz0 uz wx + 
                                  fx gx uz0 wx - fy gy uz0 wx - 
                                  fz gz uz0 wx - fz gx ux wx0 - 
                                  fx gz ux wx0 + fz gy uy wx0 - 
                                  fy gz uy wx0 + fx gx uz wx0 - 
                                  fy gy uz wx0 - fz gz uz wx0 - 
                                  fz0 gy ux wy - fz gy0 ux wy - 
                                  fy0 gz ux wy - fy gz0 ux wy - 
                                  fz gy ux0 wy - fy gz ux0 wy - 
                                  fz0 gx uy wy - fz gx0 uy wy + 
                                  fx0 gz uy wy + fx gz0 uy wy - 
                                  fz gx uy0 wy + fx gz uy0 wy + 
                                  fy0 gx uz wy + fy gx0 uz wy + 
                                  fx0 gy uz wy + fx gy0 uz wy + 
                                  fy gx uz0 wy + fx gy uz0 wy - 
                                  fz gy ux wy0 - fy gz ux wy0 - 
                                  fz gx uy wy0 + fx gz uy wy0 + 
                                  fy gx uz wy0 + fx gy uz wy0 + 
                                  fx0 gx ux wz + fx gx0 ux wz + 
                                  fy0 gy ux wz + fy gy0 ux wz - 
                                  fz0 gz ux wz - fz gz0 ux wz + 
                                  fx gx ux0 wz + fy gy ux0 wz - 
                                  fz gz ux0 wz + fy0 gx uy wz + 
                                  fy gx0 uy wz - fx0 gy uy wz - 
                                  fx gy0 uy wz + fy gx uy0 wz - 
                                  fx gy uy0 wz + fz0 gx uz wz + 
                                  fz gx0 uz wz + fx0 gz uz wz + 
                                  fx gz0 uz wz + fz gx uz0 wz + 
                                  fx gz uz0 wz + fx gx ux wz0 + 
                                  fy gy ux wz0 - fz gz ux wz0 + 
                                  fy gx uy wz0 - fx gy uy wz0 + 
                                  fz gx uz wz0 + 
                                  fx gz uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)])), 
        Cos[(θ/2)] (Cos[(ϕ/2)] (uz0 Cos[(ψ/2)] Sin[(ξ/2)] + (fz0 Cos[(ξ/2)] + (-fy0 ux - fy ux0 + fx0 uy + 
                                  fx uy0) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] (Cos[(ψ/2)] (wz0 Cos[(ξ/2)] + (uy0 wx + uy wx0 - ux0 wy - 
                                  ux wy0) Sin[(ξ/2)]) + ((fy0 wx + fy wx0 - fx0 wy - fx wy0) Cos[(ξ/2)] + (fz0 ux wx + 
                                  fz ux0 wx - fx0 uz wx - fx uz0 wx + 
                                  fz ux wx0 - fx uz wx0 + fz0 uy wy + 
                                  fz uy0 wy - fy0 uz wy - fy uz0 wy + 
                                  fz uy wy0 - fy uz wy0 - fx0 ux wz - 
                                  fx ux0 wz - fy0 uy wz - fy uy0 wz - 
                                  fz0 uz wz - fz uz0 wz - fx ux wz0 - 
                                  fy uy wz0 - 
                                  fz uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
          Sin[(θ/2)] (Cos[(ϕ/2)] (Cos[(ψ/2)] (gz0 Cos[(ξ/2)] + (-gy0 ux - gy ux0 + gx0 uy + 
                                  gx uy0) Sin[(ξ/2)]) + ((fy0 gx + fy gx0 - fx0 gy - fx gy0) Cos[(ξ/2)] + (fz0 gx ux + 
                                  fz gx0 ux - fx0 gz ux - fx gz0 ux + 
                                  fz gx ux0 - fx gz ux0 + fz0 gy uy + 
                                  fz gy0 uy - fy0 gz uy - fy gz0 uy + 
                                  fz gy uy0 - fy gz uy0 - fx0 gx uz - 
                                  fx gx0 uz - fy0 gy uz - fy gy0 uz - 
                                  fz0 gz uz - fz gz0 uz - fx gx uz0 - 
                                  fy gy uz0 - 
                                  fz gz uz0) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] (Cos[(ψ/2)] ((-gy0 wx - 
                                  gy wx0 + gx0 wy + 
                                  gx wy0) Cos[(ξ/2)] + (-gz0 ux wx - gz ux0 wx - gx0 uz wx - gx uz0 wx - gz ux wx0 - gx uz wx0 - 
                                  gz0 uy wy - gz uy0 wy - gy0 uz wy - 
                                  gy uz0 wy - gz uy wy0 - gy uz wy0 + 
                                  gx0 ux wz + gx ux0 wz + gy0 uy wz + 
                                  gy uy0 wz - gz0 uz wz - gz uz0 wz + 
                                  gx ux wz0 + gy uy wz0 - 
                                  gz uz wz0) Sin[(ξ/2)]) + ((-fz0 gx wx - fz gx0 wx - fx0 gz wx - fx gz0 wx - 
                                  fz gx wx0 - fx gz wx0 - fz0 gy wy - 
                                  fz gy0 wy - fy0 gz wy - fy gz0 wy - 
                                  fz gy wy0 - fy gz wy0 + fx0 gx wz + 
                                  fx gx0 wz + fy0 gy wz + fy gy0 wz - 
                                  fz0 gz wz - fz gz0 wz + fx gx wz0 + 
                                  fy gy wz0 - 
                                  fz gz wz0) Cos[(ξ/2)] + (fy0 gx ux wx + fy gx0 ux wx + fx0 gy ux wx + fx gy0 ux wx + 
                                  fy gx ux0 wx + fx gy ux0 wx - 
                                  fx0 gx uy wx - fx gx0 uy wx + 
                                  fy0 gy uy wx + fy gy0 uy wx + 
                                  fz0 gz uy wx + fz gz0 uy wx - 
                                  fx gx uy0 wx + fy gy uy0 wx + 
                                  fz gz uy0 wx + fz0 gy uz wx + 
                                  fz gy0 uz wx - fy0 gz uz wx - 
                                  fy gz0 uz wx + fz gy uz0 wx - 
                                  fy gz uz0 wx + fy gx ux wx0 + 
                                  fx gy ux wx0 - fx gx uy wx0 + 
                                  fy gy uy wx0 + fz gz uy wx0 + 
                                  fz gy uz wx0 - fy gz uz wx0 - 
                                  fx0 gx ux wy - fx gx0 ux wy + 
                                  fy0 gy ux wy + fy gy0 ux wy - 
                                  fz0 gz ux wy - fz gz0 ux wy - 
                                  fx gx ux0 wy + fy gy ux0 wy - 
                                  fz gz ux0 wy - fy0 gx uy wy - 
                                  fy gx0 uy wy - fx0 gy uy wy - 
                                  fx gy0 uy wy - fy gx uy0 wy - 
                                  fx gy uy0 wy - fz0 gx uz wy - 
                                  fz gx0 uz wy + fx0 gz uz wy + 
                                  fx gz0 uz wy - fz gx uz0 wy + 
                                  fx gz uz0 wy - fx gx ux wy0 + 
                                  fy gy ux wy0 - fz gz ux wy0 - 
                                  fy gx uy wy0 - fx gy uy wy0 - 
                                  fz gx uz wy0 + fx gz uz wy0 + 
                                  fz0 gy ux wz + fz gy0 ux wz + 
                                  fy0 gz ux wz + fy gz0 ux wz + 
                                  fz gy ux0 wz + fy gz ux0 wz - 
                                  fz0 gx uy wz - fz gx0 uy wz - 
                                  fx0 gz uy wz - fx gz0 uy wz - 
                                  fz gx uy0 wz - fx gz uy0 wz + 
                                  fy0 gx uz wz + fy gx0 uz wz - 
                                  fx0 gy uz wz - fx gy0 uz wz + 
                                  fy gx uz0 wz - fx gy uz0 wz + 
                                  fz gy ux wz0 + fy gz ux wz0 - 
                                  fz gx uy wz0 - fx gz uy wz0 + 
                                  fy gx uz wz0 - 
                                  fx gy uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)])), 
        Cos[(θ/2)] ((-fx0 ux - fx ux0 - fy0 uy - fy uy0 - 
                      fz0 uz - 
                      fz uz0) Cos[(ϕ/2)] Sin[(ξ/2)] Sin[(ψ/2)] + 
                Sin[(ϕ/2)] ((-ux0 wx - ux wx0 - uy0 wy - 
                            uy wy0 - uz0 wz - 
                            uz wz0) Cos[(ψ/2)] Sin[(ξ/2)] + ((-fx0 wx - fx wx0 - fy0 wy - fy wy0 - fz0 wz - 
                                  fz wz0) Cos[(ξ/2)] + (fz0 uy wx + fz uy0 wx - fy0 uz wx - fy uz0 wx + fz uy wx0 - fy uz wx0 - 
                                  fz0 ux wy - fz ux0 wy + fx0 uz wy + 
                                  fx uz0 wy - fz ux wy0 + fx uz wy0 + 
                                  fy0 ux wz + fy ux0 wz - fx0 uy wz - 
                                  fx uy0 wz + fy ux wz0 - 
                                  
                                  fx uy wz0) Sin[(ξ/2)]) Sin[(ψ/2)])) + 
          Sin[(θ/2)] (Cos[(ϕ/2)] ((-gx0 ux - gx ux0 - 
                            gy0 uy - gy uy0 - gz0 uz - 
                            gz uz0) Cos[(ψ/2)] Sin[(ξ/2)] + ((-fx0 gx - fx gx0 - fy0 gy - fy gy0 - fz0 gz - 
                                  fz gz0) Cos[(ξ/2)] + (-fz0 gy ux - fz gy0 ux + fy0 gz ux + fy gz0 ux - fz gy ux0 + fy gz ux0 + 
                                  fz0 gx uy + fz gx0 uy - fx0 gz uy - 
                                  fx gz0 uy + fz gx uy0 - fx gz uy0 - 
                                  fy0 gx uz - fy gx0 uz + fx0 gy uz + 
                                  fx gy0 uz - fy gx uz0 + 
                                  fx gy uz0) Sin[(ξ/2)]) Sin[(ψ/2)]) + 
                Sin[(ϕ/2)] (Cos[(ψ/2)] ((-gx0 wx - 
                                  gx wx0 - gy0 wy - gy wy0 - gz0 wz - 
                                  gz wz0) Cos[(ξ/2)] + (-gz0 uy wx - gz uy0 wx + gy0 uz wx + gy uz0 wx - gz uy wx0 + gy uz wx0 + 
                                  gz0 ux wy + gz ux0 wy - gx0 uz wy - 
                                  gx uz0 wy + gz ux wy0 - gx uz wy0 - 
                                  gy0 ux wz - gy ux0 wz + gx0 uy wz + 
                                  gx uy0 wz - gy ux wz0 + 
                                  gx uy wz0) Sin[(ξ/2)]) + ((fz0 gy wx + fz gy0 wx - fy0 gz wx - fy gz0 wx + fz gy wx0 - 
                                  fy gz wx0 - fz0 gx wy - fz gx0 wy + 
                                  fx0 gz wy + fx gz0 wy - fz gx wy0 + 
                                  fx gz wy0 + fy0 gx wz + fy gx0 wz - 
                                  fx0 gy wz - fx gy0 wz + fy gx wz0 - 
                                  fx gy wz0) Cos[(ξ/2)] + (fx0 gx ux wx + fx gx0 ux wx - fy0 gy ux wx - fy gy0 ux wx - 
                                  fz0 gz ux wx - fz gz0 ux wx + 
                                  fx gx ux0 wx - fy gy ux0 wx - 
                                  fz gz ux0 wx + fy0 gx uy wx + 
                                  fy gx0 uy wx + fx0 gy uy wx + 
                                  fx gy0 uy wx + fy gx uy0 wx + 
                                  fx gy uy0 wx + fz0 gx uz wx + 
                                  fz gx0 uz wx + fx0 gz uz wx + 
                                  fx gz0 uz wx + fz gx uz0 wx + 
                                  fx gz uz0 wx + fx gx ux wx0 - 
                                  fy gy ux wx0 - fz gz ux wx0 + 
                                  fy gx uy wx0 + fx gy uy wx0 + 
                                  fz gx uz wx0 + fx gz uz wx0 + 
                                  fy0 gx ux wy + fy gx0 ux wy + 
                                  fx0 gy ux wy + fx gy0 ux wy + 
                                  fy gx ux0 wy + fx gy ux0 wy - 
                                  fx0 gx uy wy - fx gx0 uy wy + 
                                  fy0 gy uy wy + fy gy0 uy wy - 
                                  fz0 gz uy wy - fz gz0 uy wy - 
                                  fx gx uy0 wy + fy gy uy0 wy - 
                                  fz gz uy0 wy + fz0 gy uz wy + 
                                  fz gy0 uz wy + fy0 gz uz wy + 
                                  fy gz0 uz wy + fz gy uz0 wy + 
                                  fy gz uz0 wy + fy gx ux wy0 + 
                                  fx gy ux wy0 - fx gx uy wy0 + 
                                  fy gy uy wy0 - fz gz uy wy0 + 
                                  fz gy uz wy0 + fy gz uz wy0 + 
                                  fz0 gx ux wz + fz gx0 ux wz + 
                                  fx0 gz ux wz + fx gz0 ux wz + 
                                  fz gx ux0 wz + fx gz ux0 wz + 
                                  fz0 gy uy wz + fz gy0 uy wz + 
                                  fy0 gz uy wz + fy gz0 uy wz + 
                                  fz gy uy0 wz + fy gz uy0 wz - 
                                  fx0 gx uz wz - fx gx0 uz wz - 
                                  fy0 gy uz wz - fy gy0 uz wz + 
                                  fz0 gz uz wz + fz gz0 uz wz - 
                                  fx gx uz0 wz - fy gy uz0 wz + 
                                  fz gz uz0 wz + fz gx ux wz0 + 
                                  fx gz ux wz0 + fz gy uy wz0 + 
                                  fy gz uy wz0 - fx gx uz wz0 - 
                                  fy gy uz wz0 + 
                                  fz gz uz wz0) Sin[(ξ/2)]) Sin[(ψ/2)]))}

Substitute data

ComputeNum[A_, g1num_, g2num_, g3num_, g4num_, pnum_, theta_, phi_, 
        psi_, xi_] := 
      Simplify[A /. {gx -> g1num[[1, 1]], 
            gy -> g1num[[1, 2]], gz -> g1num[[1, 3]], 
            gx0 -> g1num[[2, 1]], gy0 -> g1num[[2, 2]], 
            gz0 -> g1num[[2, 3]], wx -> g2num[[1, 1]], 
            wy -> g2num[[1, 2]], wz -> g2num[[1, 3]], 
            wx0 -> g2num[[2, 1]], wy0 -> g2num[[2, 2]], 
            wz0 -> g2num[[2, 3]], fx -> g3num[[1, 1]], 
            fy -> g3num[[1, 2]], fz -> g3num[[1, 3]], 
            fx0 -> g3num[[2, 1]], fy0 -> g3num[[2, 2]], 
            fz0 -> g3num[[2, 3]], ux -> g4num[[1, 1]], 
            uy -> g4num[[1, 2]], uz -> g4num[[1, 3]], 
            ux0 -> g4num[[2, 1]], uy0 -> g4num[[2, 2]], 
            uz0 -> g4num[[2, 3]], px -> pnum[[1, 1]], 
            py -> pnum[[1, 2]], pz -> pnum[[1, 3]], 
            pw -> pnum[[1, 4]], px0 -> pnum[[2, 1]], 
            py0 -> pnum[[2, 2]], pz0 -> pnum[[2, 3]], 
            pw0 -> pnum[[2, 4]], θ -> 
              theta, ϕ -> phi, ψ -> psi, ξ -> 
              xi}]

ComputeNum[quat4r, Gline, Wline, Fline, Uline, qS12, thetas[[1]], 
      phis[[1]], psis[[1]], xis[[1]]]
qS12

ComputeNum[quat4r, Gline, Wline, Fline, Uline, qS13, thetas[[2]], 
      phis[[2]], psis[[2]], xis[[2]]]
qS13

ComputeNum[quat4r, Gline, Wline, Fline, Uline, qS14, thetas[[3]], 
      phis[[3]], psis[[3]], xis[[3]]]
qS14

ComputeNum[quat4r, Gline, Wline, Fline, Uline, qS15, thetas[[4]], 
      phis[[4]], psis[[4]], xis[[4]]]
qS15

ComputeNum[quat4r, Gline, Wline, Fline, Uline, qS16, thetas[[5]], 
      phis[[5]], psis[[5]], xis[[5]]]
qS16

ComputeNum[quat4r, Gline, Wline, Fline, Uline, qS17, thetas[[6]], 
      phis[[6]], psis[[6]], xis[[6]]]
qS17

ComputeNum[quat4r, Gline, Wline, Fline, Uline, qS18, thetas[[7]], 
      phis[[7]], psis[[7]], xis[[7]]]
qS18

ComputeNum[quat4r, Gline, Wline, Fline, Uline, qS19, thetas[[8]], 
      phis[[8]], psis[[8]], xis[[8]]]
qS19