################################################################################
#                                                                              #
#                    Package sui quaternioni                                   #
#                    -----------------------                                   #
#                                                                              #
# Autore: Claudio Marsan                                                       #
# Ultima revisione: 16 marzo 2003                                              #
# Versione: Maple V Release 6.02 for Windows 2000                              #
#                                                                              #
################################################################################

# Base per i quaternioni
alias( q1 = quat(1, 0, 0, 0),
       qi = quat(0, 1, 0, 0),
       qj = quat(0, 0, 1, 0),
       qk = quat(0, 0, 0, 1)):

# "Pretty print" per i quaternioni
`print/quat` := proc(a0, a1, a2, a3)
  a0 + a1*'i' + a2*'j' + a3*'k';
end:

# Estrazione dei quattro coefficienti dalla struttura di dati di un quaternione
quaternion[qparts] := proc(a, a0, a1, a2, a3)
  a0 := op(1, a);
  a1 := op(2, a);
  a2 := op(3, a);
  a3 := op(4, a);
  NULL;
end:

# Traccia di un quaternione
quaternion[qtrace] := proc(a)
  local a0, a1, a2, a3;
  quaternion[qparts](a, a0, a1, a2, a3);
  2*a0;
end:

# Norma di un quaternione
quaternion[qnorm] := proc(a)
  local a0, a1, a2, a3;
  quaternion[qparts](a, a0, a1, a2, a3);
  a0^2 + a1^2 + a2^2 + a3^2;
end:

# Valore assoluto di un quaternione
quaternion[qabs] := proc(a)
  local a0, a1, a2, a3;
  quaternion[qparts](a, a0, a1, a2, a3);
  sqrt(a0^2 + a1^2 + a2^2 + a3^2);
end:

# Coniugato di un quaternione
quaternion[qconj] := proc(a)
  local a0, a1, a2, a3;
  quaternion[qparts](a, a0, a1, a2, a3);
  quat(a0, -a1, -a2, -a3);
end:

# Prodotto di uno scalare con un quaternione:
quaternion[qscal] := proc(k, a)
  local a0, a1, a2, a3;
  quaternion[qparts](a, a0, a1, a2, a3);
  quat(k*a0, k*a1, k*a2, k*a3);
end:

# Somma di quaternioni
quaternion[qadd] := proc()
  local n, lista, nuovalista, i, a0, a1, a2, a3, b0, b1, b2, b3 ;
  n := nargs;
  lista := [args];
  if n = 2
    then
      quaternion[qparts](lista[1], a0, a1, a2, a3);
      quaternion[qparts](lista[2], b0, b1, b2, b3);
      quat(a0 + b0, a1 + b1, a2 + b2, a3 + b3);
    else
      quaternion[qparts](lista[1], a0, a1, a2, a3);
      quaternion[qparts](lista[2], b0, b1, b2, b3);
      nuovalista := quat(a0 + b0, a1 + b1, a2 + b2, a3 + b3),
		    seq(lista[i], i = 3..n);
      quaternion[qadd](nuovalista);
  fi;
end:

# Differenza di due quaternioni
quaternion[qsub] := proc(a, b)
  local a0, a1, a2, a3, b0, b1, b2, b3;
  quaternion[qparts](a, a0, a1, a2, a3);
  quaternion[qparts](b, b0, b1, b2, b3);
  quat(a0 - b0, a1 - b1, a2 - b2, a3 - b3);
end:

# Prodotto di quaternioni
quaternion[qmul] := proc()
  local n, lista, nuovalista, i, a0, a1, a2, a3, b0, b1, b2, b3 ;
  n := nargs;
  lista := [args];
  if n = 2
    then
      quaternion[qparts](lista[1], a0, a1, a2, a3);
      quaternion[qparts](lista[2], b0, b1, b2, b3);
      quat(a0*b0 - a1*b1 - a2*b2 - a3*b3,
	   a0*b1 + a1*b0 + a2*b3 - a3*b2,
	   a0*b2 - a1*b3 + a2*b0 + a3*b1,
	   a0*b3 + a1*b2 - a2*b1 + a3*b0);
    else
      quaternion[qparts](lista[1], a0, a1, a2, a3);
      quaternion[qparts](lista[2], b0, b1, b2, b3);
      nuovalista := quat(a0*b0 - a1*b1 - a2*b2 - a3*b3,
			 a0*b1 + a1*b0 + a2*b3 - a3*b2,
			 a0*b2 - a1*b3 + a2*b0 + a3*b1,
			 a0*b3 + a1*b2 - a2*b1 + a3*b0),
		    seq(lista[i], i = 3..n);
      quaternion[qmul](nuovalista);
  fi;
end:

# Inversa moltiplicativa di un quaternione non nullo
quaternion[qinv] := proc(a)
  local a0, a1, a2, a3, q;
  quaternion[qparts](a, a0, a1, a2, a3);
  q := quaternion[qnorm](a);
  if q = 0
    then ERROR(`Il quaternione 0 non e' invertibile!`);
    else quat(a0/q, -a1/q, -a2/q, -a3/q);
  fi;
end:

# Testo di aiuto per il package dei quaternioni
`help/text/quaternion` := TEXT(
` `,
`HELP FOR:  quaternion`,
` `,
`          Introduzione al package dei quaternioni`,
`          ---------------------------------------`,
` `,
`Il seguente package permette di operare con i quaternioni classici, ossia`,
`della forma  a0 + a1*i + a2*j + a3*k (a0, a1, a2, a3 sono scalari reali) e `,
`con le usuali regole per i prodotti tra i, j, k.`,
`Un quaternione e' definito mediante una struttura di dati della forma`,
`quat(a0, a1, a2, a3), dove a0, a1, a2, a3 sono scalari reali rappresentanti`,
`i coefficienti di, rispettivamente, 1, i, j, k.`,
`I quaternioni speciali 1, i, j, k possono essere richiamati mediante, `,
`rispettivamente, q1, qi, qj, qk.`,
`I risultati vengono restituiti nella forma a0 + a1*i + a2*j + a3*k.`,
` `,
`REMARK:`,
`Come punto di partenza per la programmazione di questo package e' stato`,
`utilizzato il package di Greg Fee dell'aprile 1992, pubblicato su "First`,
`Leaves: A Tutorial Introduction to Maple V" che fa parte della documentazione`,
`ufficiale di Maple V. Questo package contiene molte piu' funzioni di quello`,
`originale che conteneva solo qparts, qadd, qmul e inoltre qadd e qmul sono`,
`state estese per poter supportare la somma e il prodotto di piu' di due`,
`quaternioni con una sola chiamata di funzione.`,
` `,
`SEE ALSO:`,
`    qparts, qtrace, qnorm, qabs, qconj, qscal, qadd, qsubs, qmul, qinv`):

# Testo di aiuto per la funzione qparts
`help/text/qparts` := TEXT(
` `,
`HELP FOR: qparts`,
` `,
`CALLING SEQUENCE:`,
`    qparts(a, a0, a1, a2, a3)`,
` `,
`PARAMETERS:`,
`    - a e' un quaternione`,
`    - a0, a1, a2, a3 sono le variabili alle quali verranno assegnati i `,
`      coefficienti di a`,
` `,
`SYNONPSIS:`,
`qparts estrae i coefficienti del quaternione a e li assegna alle variabili`,
`a0, a1, a2, a3. Non ritorna alcun valore sul terminale.`,
` `,
`EXAMPLE:`,
`> with(quaternion):`,
`> a := quat(3, -1, 4, 2);`,
`                     3 - i + 4j + 2k`,
`> qparts(a, a0, a1, a2, a3);`,
` `,
`> a0; a1; a2; a3;`,
`                     3`,
`                     -1`,
`                     4`,
`                     2`,
` `,
`SEE ALSO:`,
`    quaternion`):

# Testo di aiuto per la funzione qtrace
`help/text/qtrace` := TEXT(
` `,
`HELP FOR: qtrace`,
` `,
`CALLING SEQUENCE:`,
`    qtrace(a)`,
` `,
`PARAMETERS:`,
`    a e' un quaternione`,
` `,
`SYNONPSIS:`,
`qtrace calcola la traccia del quaternione a. Il risultato e' un reale.`,
` `,
`EXAMPLE:`,
`> with(quaternion):`,
`> a := quat(3, -1, 4, 2);`,
`                     3 - i + 4j + 2k`,
`> qtrace(a);`,
`                     6`,
` `,
`SEE ALSO:`,
`    quaternion`):

# Testo di aiuto per la funzione qnorm
`help/text/qnorm` := TEXT(
` `,
`HELP FOR: qnorm`,
` `,
`CALLING SEQUENCE:`,
`    qnorm(a)`,
` `,
`PARAMETERS:`,
`    a e' un quaternione`,
` `,
`SYNONPSIS:`,
`qnorm calcola la norma del quaternione a. Il risultato e' un reale.`,
` `,
`EXAMPLE:`,
`> with(quaternion):`,
`> a := quat(3, -1, 4, 2);`,
`                     3 - i + 4j + 2k`,
`> qnorm(a);`,
`                     30`,
` `,
`SEE ALSO:`,
`    quaternion`):

# Testo di aiuto per la funzione qabs
`help/text/qabs` := TEXT(
` `,
`HELP FOR: qabs`,
` `,
`CALLING SEQUENCE:`,
`    qabs(a)`,
` `,
`PARAMETERS:`,
`    a e' un quaternione`,
` `,
`SYNONPSIS:`,
`qabs calcola il valore assoluto del quaternione a. Il risultato e' un reale.`,
` `,
`EXAMPLE:`,
`> with(quaternion):`,
`> a := quat(3, -1, 4, 2);`,
`                     3 - i + 4j + 2k`,
`> qabs(a);`,
`                     sqrt(30)`,
` `,
`SEE ALSO:`,
`    quaternion`):

# Testo di aiuto per la funzione qconj
`help/text/qconj` := TEXT(
` `,
`HELP FOR: qconj`,
` `,
`CALLING SEQUENCE:`,
`    qconj(a)`,
` `,
`PARAMETERS:`,
`    a e' un quaternione`,
` `,
`SYNONPSIS:`,
`qconj calcola il quaternione coniugato di a. Il risultato e' un quaternione.`,
` `,
`EXAMPLE:`,
`> with(quaternion):`,
`> a := quat(3, -1, 4, 2);`,
`                     3 - i + 4j + 2k`,
`> qconj(a);`,
`                     3 + i - 4j - 2k`,
` `,
`SEE ALSO:`,
`    quaternion`):

# Testo di aiuto per la funzione qscal
`help/text/qscal` := TEXT(
` `,
`HELP FOR: qscal`,
` `,
`CALLING SEQUENCE:`,
`    qscal(k, a)`,
` `,
`PARAMETERS:`,
`    - k e' uno scalare (un numero reale)`,
`    - a e' un quaternione`,
` `,
`SYNONPSIS:`,
`qscal calcola il prodotto tra lo scalare k e il quaternione a. Il risultato`,
`e' un quaternione`,
` `,
`EXAMPLE:`,
`> with(quaternion):`,
`> a := quat(3, -1, 4, 2);`,
`                     3 - i + 4j + 2k`,
`> qscal(-5, a);`,
`                     -15 + 5i - 20j - 10k`,
` `,
`SEE ALSO:`,
`    quaternion`):

# Testo di aiuto per la funzione qadd
`help/text/qadd` := TEXT(
` `,
`HELP FOR: qadd`,
` `,
`CALLING SEQUENCE:`,
`    qadd(args)`,
` `,
`PARAMETERS:`,
`    args e' una successione di quaternioni`,
` `,
`SYNONPSIS:`,
`qadd somma i quaternioni che appaiono nella successione degli argomenti. Il`,
`risultato e' un quaternione`,
` `,
`EXAMPLE:`,
`> with(quaternion):`,
`> a := quat(2, 3, 1, 5);`,
`                       2 + 3i + j + 5k`,
`> b := quat(9, 4, 7, 6);`,
`                       9 + 4i + 7j + 6k`,
`> c := quat(0, 2, -1, 2);`,
`                       2i - j + 2k`,
`> qadd(a, b, c);`,
`                       11 + 9i + 7j + 13k`,
` `,
`SEE ALSO:`,
`    quaternion`):

# Testo di aiuto per la funzione qsub
`help/text/qsub` := TEXT(
` `,
`HELP FOR: qsub`,
` `,
`CALLING SEQUENCE:`,
`    qsub(a, b)`,
` `,
`PARAMETERS:`,
`    a, b sono due quaternioni`,
` `,
`SYNONPSIS:`,
`qsub calcola a-b. Il risultato e' un quaternione`,
` `,
`EXAMPLE:`,
`> with(quaternion):`,
`> a := quat(2, 3, 1, 5);`,
`                       2 + 3i + j + 5k`,
`> b := quat(9, 4, 7, 6);`,
`                       9 + 4i + 7j + 6k`,
`> qsub(a, b);`,
`                       -7 - 5i - 6j - k`,
` `,
`SEE ALSO:`,
`    quaternion`):

# Testo di aiuto per la funzione qmul
`help/text/qmul` := TEXT(
` `,
`HELP FOR: qmul`,
` `,
`CALLING SEQUENCE:`,
`    qmul(args)`,
` `,
`PARAMETERS:`,
`    args e' una successione di quaternioni`,
` `,
`SYNONPSIS:`,
`qmul moltiplica, da sinistra a destra, i quaternioni che appaiono nella`,
`successione degli argomenti. Il risultato e' un quaternione.`,
` `,
`EXAMPLE:`,
`> with(quaternion):`,
`> a := quat(2, 3, 1, 5);`,
`                       2 + 3i + j + 5k`,
`> b := quat(9, 4, 7, 6);`,
`                       9 + 4i + 7j + 6k`,
`> c := quat(0, 2, -1, 2);`,
`                       2i - j + 2k`,
`> qmul(a, b, c);`,
`                       -135 + 62i + 167j - 118k`,
` `,
`SEE ALSO:`,
`    quaternion`):

# Testo di aiuto per la funzione qinv
`help/text/qinv` := TEXT(
` `,
`HELP FOR: qinv`,
` `,
`CALLING SEQUENCE:`,
`    qinv(a)`,
` `,
`PARAMETERS:`,
`    a e' un quaternione`,
` `,
`SYNONPSIS:`,
`qconj calcola il quaternione inverso di a, se a e' diverso da zero. Il`,
`risultato e' un quaternione.`,
` `,
`EXAMPLE:`,
`> with(quaternion):`,
`> a := quat(3, -1, 4, 2);`,
`                     3 - i + 4j + 2k`,
`> qinv(a);`,
`                     1/10 + i/30 - 2j/15 - k/15`,
` `,
`SEE ALSO:`,
`    quaternion`):