/*
Uno studio per la Befana 2015 e rifatto nel 2016...
( mancano 84 anni al 2100 )
*/

/*
Lorentz in modo semplicissimo
*/

/*
Prima variante:
*/

/*
Questa funzione crea la matrice di Lorentz facendo 
pero' uso della funzione sqrt(...) quindi con un piccolo 
rischio di complicazioni numeriche ossia che Maxima si 
incasini e non riesca a semplificare le espressioni.
---
Questa funzione ha bisogno di due dati:
1) Il vettore che specifica la direzione della trasformazione
ossia un vettore di 4 elementi di cui almeno uno dei primi 
tre elementi non deve essere nullo e l'ultimo, che controlla 
la formula globale, DEVE ESSERE OBBLIGATORIAMENTE POSITIVO,
MAI NULLO.
2) La velocita' della luce ossia un numero positivo e mai nullo.
---
Questa funzione produce tre risultati ossia :
1) La trasformazione di Lorentz.
2) Il quadrivettore che viene "fermato" dalla trasformazione
di Lorentz ossia la trasformazione di Lorentz moltiplicata
per questo vettore da' un vettore con le componenti spaziali
tutte nulle.
3) Il tensore metrico covariante che si usa e che dipende
da quale e' considerata la velocita' della luce nel vuoto.
---
*/

lsemplice(br,lux):=block([bd,bx,by,bz,cm,sh,lz,uu,gm],
bd:sqrt(br[1]*br[1]+br[2]*br[2]+br[3]*br[3]),
bx:br[1]/bd,
by:br[2]/bd,
bz:br[3]/bd,
cm:(br[4]+1/br[4])/2-1,
sh:(br[4]-1/br[4])/2,
lz:matrix([1+cm,-sh*bx/lux,-sh*by/lux,-sh*bz/lux],
[-lux*sh*bx,1+cm*bx*bx,cm*bx*by,cm*bx*bz],
[-lux*sh*by,cm*bx*by,1+cm*by*by,cm*by*bz],
[-lux*sh*bz,cm*bx*bz,cm*by*bz,1+cm*bz*bz]),
uu:[1+cm,lux*sh*bx,lux*sh*by,lux*sh*bz],
gm:matrix([lux*lux,0,0,0],
[0,-1,0,0],[0,0,-1,0],[0,0,0,-1]),
return ([lz,uu,gm]))$

/*
Verifico che la funzione data funzioni bene scegliendo 
bene i dati in modo da non ottenere radici quadrate.
Quindi Maxima non deve fare praticamente nessuna fatica 
a semplificare i risultati...
Definisco tre vettori il cui significato sara' 
spiegato poi...
*/

cluce:7;
besea:[3/13,4/13,12/13,9];
beseaq:[3/13,4/13,12/13,81];
beseasq:[3/13,4/13,12/13,3];

/*
Ecco il calcolo applicato al primo vettore....
*/

lzz:lsemplice(besea,cluce);

/*
Verifico che la trasformazione di Lorentz, applicata al 
quadrivettore calcolato, lo fermi veramente ed inoltre la 
trasformazione di Lorentz, applicata al tensore metrico 
covariante, non lo modifica ossia lo lascia immutato...
---
Attenzione. Se la velocita' della luce non vale 1, 
la matrice della trasformazione di Lorentz 
NON E' SIMMETRICA per cui va usata anche la sua 
trasposta ottenuta con la transpose(...) di Maxima.
*/

print("La trasformata di Lorentz e' lzz[1] ",lzz[1])$
print("lzz[2] e' il vettore usato= ",lzz[2])$
print("La norma di lzz[2]= ",lzz[2].lzz[3].lzz[2])$
print("Usando lzz[1] il vettore dovrebbe essere ",
  "FERMO ossia zero tutte le sue componenti spaziali")$
unuovo:lzz[1].lzz[2]$
print("unuovo = ",unuovo)$
print("La norma di unuovo deve essere invariata ",
  "ossia unuovo.lzz[3].unuovo = ",
unuovo.lzz[3].unuovo)$
print("Come e' il tensore metrico all'inizio =",
   lzz[3])$
print("Dopo la trasformazione il tens.metrico",
      " deve essere identico")$
transpose(lzz[1]).lzz[3].lzz[1];
print("Ho scritto: transpose(lzz[1]).lzz[3].lzz[1] ")$

/*
Qui mostro come fare il quadrato della trasformazione di
Lorentz ossia cosa viene, applicata due volte di seguito...
*/

lzzq:lsemplice(beseaq,cluce);

/*
Dimostro che effettivamente lzzq[1] e' il quadrato di lzz[1].
*/

lzzq[1]-lzz[1].lzz[1];

/*
Dimostro come si fa la radice quadrata della matrice che 
rappresenta la trasformazione di Lorentz ossia lzzsq[1] 
moltiplicata per se stessa deve dare lzz[1]...
*/

lzzsq:lsemplice(beseasq,cluce);

/*
Dimostro che realmente lzzsq[1] e' la radice quadrata di lzz[1]...
*/

lzz[1]-lzzsq[1].lzzsq[1];

/*
Seconda variante: 
W la razionalita' !!!
*/

/*
Questa funzione e' simile a quella precedente ma ha 
il GRANDE PREGIO DI NON FARE USO DELLA sqrt(...).
Per evitare problemi impedisco di usare come primo
argomento il valore zero che darebbe divisione per zero
e dunque dati infiniti ossia errore di calcolo....
Il secondo ed il terzo argomento possono essere dati 
qualsiasi ed il quarto e' la velocita' della luce nel 
vuoto che deve ovviamente non essere mai uno zero..
Se si assegna zero viene usato il valore 1.
*/

lrazionale(rd,f,h,xlux):=block([r,lux,bx,by,bz,cm,sh,
    lz,uu,gm],
if rd=0 then r:1 else r:rd,
if xlux=0 then lux:1 else lux:abs(xlux),
bx:2*h*(1-f*f)/((1+f*f)*(1+h*h)),
by:4*h*f/((1+f*f)*(1+h*h)),
bz:(1-h*h)/(1+h*h),
cm:(r+1/r)/2-1,
sh:(r-1/r)/2,
lz:matrix([1+cm,-sh*bx/lux,-sh*by/lux,-sh*bz/lux],
[-lux*sh*bx,1+cm*bx*bx,cm*bx*by,cm*bx*bz],
[-lux*sh*by,cm*bx*by,1+cm*by*by,cm*by*bz],
[-lux*sh*bz,cm*bx*bz,cm*by*bz,1+cm*bz*bz]),
uu:[1+cm,lux*sh*bx,lux*sh*by,lux*sh*bz],
gm:matrix([lux*lux,0,0,0],
[0,-1,0,0],[0,0,-1,0],[0,0,0,-1]),
return ([lz,uu,gm]))$

/*
Ecco l'applicazione della funzione che produce 
sempre numeri razionali se gli argomenti di ingresso 
sono anche loro razionali...
*/

lrr:lrazionale(1+random(7)/11,random(19)/17,
    random(23)/19,cluce);
lrr[2].lrr[3].lrr[2];
urnuovo:lrr[1].lrr[2];
urnuovo.lrr[3].urnuovo;
print("Come e' il tensore metrico all'inizio")$
lrr[3];
print("Dopo la trasformazione ",
   " il tens.metrico deve essere identico")$
transpose(lrr[1]).lrr[3].lrr[1];

lrra:lrazionale(25,4,7,cluce);
lrraq:lrazionale(25*25,4,7,cluce);
lrrasq:lrazionale(5,4,7,cluce);

/*
Bene se ottengo solo zeri...
*/

lrra[1].lrra[1]-lrraq[1];
lrrasq[1].lrrasq[1]-lrra[1];

/*
Altre funzioni di servizio, utili per gestire i vettori
con le convenzioni usate qui...
*/

/*
Questa piccola funzione normalizza ad 1 la sola parte 
spaziale di un quadrivettore a meno che non abbia tutte 
le componenti spaziali nulle nel qual caso da' 
un vettore di 3 dimensioni di zeri. 
*/

v3norma(vv):=block([vsq,vn],
  vsq:sqrt(vv[2]^2+vv[3]^2+vv[4]^2),
  if vsq=0 then return([0,0,0]),
  vn:ratsimp([vv[2]/vsq,vv[3]/vsq,vv[4]/vsq]),
  return (vn))$

/*
In coordinate sferiche indicate da ra, f ed h puo' 
servire passare alle tradizionali cartesiane...
*/

incart(ra,f,h):=block([bx,by,bz],
  bx:2*h*(1-f*f)/((1+f*f)*(1+h*h)),
  by:4*h*f/((1+f*f)*(1+h*h)),
  bz:(1-h*h)/(1+h*h),
  return( [ra*bx,ra*by,ra*bz] ) )$

/*
Conclusione
*/

/*
Le due funzioni qui definite e sperimentate:
lsemplice(br,lux) e lrazionale(rd,f,h,xlux)
chiariscono la struttura della trasformazione di Lorentz 
e, la seconda, a mio parere, e' utilissima per lavorare 
solo con numeri razionali e quindi operando con la massima 
precisione possibile su qualsiasi calcolatore.
*/

print("Qui sono definite due funzioni:")$
print("lsemplice(br,lux)")$
print("lrazionale(rd,f,h,xlux)")$
print("E varie variabili, ad esempio: besea =",besea,
 " e cluce=",cluce)$
print("Matrici lzz[1]=",lzz[1]," e lrra[1]=",lrra[1])$

lzz[1].lzz[2];

/*
Qui sono definite anche le due funzioni ausiliarie
v3norma(vv) e incart(ra,f,h).
*/

/*
Le formule qui applicate sembrano banali, data la
loro semplicita' ma... a me sembra che siano
raramente usate in letteratura e nell'insegnamento...
*/