/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 0.8.5 ] */

/* [wxMaxima: title   start ]
Libreria tensoriale di
Giampaolo Bottoni
(crealibtensori v.1)
   [wxMaxima: title   end   ] */

/* [wxMaxima: title   start ]
Scrive su disco la libreria)
   [wxMaxima: title   end   ] */

/* [wxMaxima: comment start ]
Versione ( in via di sviluppo ) del 22 giugno 2010.
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Premessa
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
Un tensore è una lista di oggetti vari (atomici, 
liste o matrici).
L'ultimo elemento della lista deve essere una lista.
La lunghezza della lista dell'ultimo elemento
è l'ordine del tensore + 2.
Il penultimo elemento della lista dell'ultimo
elemento è un commento di qualsiasi tipo,
un numero o una stringa o una lista etc.
L'ultimo elemento della lista dell'ultimo elemento
è la dimensione dello spazio per il tensore
considerato.
I precedenti elementi indicano la natura degli
indici del tensore ossia se sono controvarianti (1),
covarianti (2), derivativi ordinari(0), derivati
covarianti covarianti (8), derivati covarianti
controvarianti (9), personalizzati ( altre cifre
escluse 0,1,2,8,9 ). 
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Con questa definizione è possibile trattare come
tensore anche uno scalare che ha questa struttura:

[valore,[commento,dimensione_spazio]]
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
La dimensione_spazio minima ammessa è 2 mentre
per descrivere lo spazio tridimensionale deve essere:

dimensione_spazio:3

Per descrivere gli spazi relativistici pseudoeuclidei
ossia con la prima dimensione temporale e le restanti
tre spaziali deve essere:

dimensione_spazio:4
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
I tensori, date le scelte di programmazione ora
fatte, si distinguono tra tensori con un numero
dispari di indici e tensori con un numero pari
di indici.
Se il tensore ha numero di indici pari viene usato
solo il primo elemento del tensore mentre se 
il numero è dispari, vengono usati i primi elementi 
da 1 a dimensione_spazio.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Un vettore è un tensore di ordine 1 di questa
struttura, nel caso in cui la dimensione
dello spazio sia 2:

[v1,v2,[i1,commento,2]]

Se la dimensione dello spazio vale 3 il vettore
ha questa struttura:

[v1,v2,v3,[i1,commento,3]]

Se la dimensione dello spazio vale 4 il vettore
ha questa struttura:

[v1,v2,v3,v4,[i1,commento,4]]

Il valore di i1 è 1 se il vettore è controvariante
ed è 2 se è covariante etc.
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Operazioni di salvataggio su file
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Specifico il file in cui scrivo la nuova libreria:
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
salvoqui:"c:/xmaxima/libtensori.mc";
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
if atom(path_iniziale) then writefile(salvoqui);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
print(["Ora comincio a lavorare su ",salvoqui])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
cartellalavoro:"C:/xmaxima/###.{mc,mac}";
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
E' consigliabile, prima di modificarlo, di salvare il valore 
del path di default di Maxima.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
( if atom(path_iniziale) then 
   (path_iniziale: file_search_maxima) )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
file_search_maxima: cons(cartellalavoro,path_iniziale);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Amplio il path iniziale aggiungendogli una cartella. Con questo 
trucco posso ricaricare varie volte questo documento senza il 
problema di modificare ogni volta il path di ricerca. 
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
file_search_maxima: cons(cartellalavoro,path_iniziale);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
print(["Il path usato è il seguente",cartellalavoro]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Funzioni algebriche
   [wxMaxima: section end   ] */

/* [wxMaxima: subsect start ]
Funzioni di informazione
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
L'ordine di un tensore è ottenibile con la seguente
funzione che, per ragioni di velocità, non effettua
nessuna verifica della correttezza del tensore tn
passato in argomento.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
ordine(tn):=block([],
    length(tn[length(tn)])-2
    )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:["ordine(tn)"];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Il commento che ogni tensore contiene e che può essere
una stringa ma anche un numero o una lista o un 
qualsiasi oggetto, è estraibile in questo modo:
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
commento(tn):=block([u],
    u:length(tn),
    tn[u][length(tn[u])-1])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("commento(tn)",libmia);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Per inserire un commento nel tensore, uso la commentalo:
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
commentalo(tn,comme):=block([u],
    u:length(tn),
    tn[u][length(tn[u])-1]:comme,tn[u] )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("commentalo(tn,comme)",libmia);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Ogni tensore contiene, come si è detto, la dimensione
dello spazio a cui si riferisce. Per ottenere questo
dato si può usare questa funzione:
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
dimensione(tn):=block([],
    tn[length(tn)][length(tn[length(tn)])])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("dimensione(tn)",libmia);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Per estrarre dal tensore la lista dei tipi (covarianti,
controvarianti, etc.) dei suoi indici si può usare
questa funzione:
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
indici(tn):=block([],
    rest(tn[length(tn)],-2))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("indici(tn)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Indichiamo con tn un generico tensore ed esaminamo
i test che ogni funzione, quando non ci sono 
esigenze di prestazioni, potrebbe fare per 
controllare che la struttura del tensore è corretta.
Raccogliamo i test nella funzione tensorp(tn)
che dà risultato true se tn è un tensore
ben fatto e false se tn non è un tensore o è un
tensore mal fatto.
La tensorp utilizza in modo ricorsivo la
matrixxp che opera su matrici di matrici.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
matrixxp(tn,o,n):=block([basta,smetti],
   basta:false,smetti:true,
   if o>1 then (
      smetti:false,
      if not(matrixp(tn)) then basta:true,
      if length(tn[1])#n then basta:true
      ),
   if basta then return(false), 
   if smetti then return(true),
   matrixxp(tn[1,1],o-2,n)
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
tensorp(tn):=block([u,o,n,pd],
  if not(listp(tn)) then return(false),
  u:length(tn),
  if not(listp(tn[u])) then return(false),
  o:length(tn[u])-2,
  if 0>o then return(false),
  n:floor(tn[u][o+2]),
  if 2>n then return(false),
  pd:mod(o,2),
  if pd=0 then ( if u>n then return(false))
  else if n>u-1 then return(false),
  if matrixxp(tn[1],o,n) then return(true),
  false )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tensorp(tn)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
zerotensor: funzione di inizializzazione
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
La funzione zerotensor genera un tensore con elementi 
tutti nulli, con indici di valore prescelto,
con commento arbitrario e di dimensione assegnata.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
zerotensor(tipi,comme,ndim):=block([tn,j,o,u,nd],
   if not(listp(tipi)) then 
      return("Errore, non lista tipi"),
   o:length(tipi), nd:floor(ndim),
   if 2>nd then 
      return("Errore, dimensione inferiore a 2"),
   if o=0 then return([0,[comme,nd]]),
   if mod(o,2)=0 then tn:makelist(0,j,1,2)
   else tn:makelist(0,j,1,nd+1),
   u:length(tn),
   tn[u]:append(makelist(tipi[j],j,1,o),[comme,nd]),
   for j:1 thru u-1 do tn[j]:zeromatmat(o,nd),
   tn
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
zeromatmat(o,nd):=block([tr,omm,j,k],
   if 2>o then return(0),
   tr:zeromatrix(nd,nd),
   omm:o-2,
   for j:1 thru nd do
   for k:1 thru nd do
       tr[j,k]:zeromatmat(omm,nd),
   tr )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("zerotensor(tipi,comme,ndim)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Funzioni per generare tensori da scalari,
liste o matrici
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Per metrica si intende una lista di tensori utili
per gestire una data metrica. 
La lista viene generata da una apposita funzione
fa_metrica(matrice,listavariabili) definita
in qui in seguito.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Crea un tensore scalare 
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
tsca(sc,metrica):=block([],
    [sc,["(scal)",dimensione(metrica[2])]])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tsca(sc,metrica)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Crea un vettore covariante, data una lista e
la metrica.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
tcov(lv,metrica):=block([nl,nd],
    if not(listp(lv)) then
       return("Non lista"),
    nl:length(lv),
    nd:dimensione(metrica[2]),
    [ makelist(lv[1+mod(j-1,nl)],j,1,nd),
      [2,"(v2)",nd]])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tcov(lv,metrica)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Crea un vettore controvariante, data la lista e
la metrica.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
tcontrov(lv,metrica):=block([nl,nd],
    if not(listp(lv)) then
       return("Non lista"),
    nl:length(lv),
    nd:dimensione(metrica[2]),
    [ makelist(lv[1+mod(j-1,nl)],j,1,nd),
      [1,"(v1)",nd]])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tcontrov(lv,metrica)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Crea un tensore di ordine 2, totalmente controvariante,
data una matrice e la metrica.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
tmat11(mat,metrica):=block([],
    if not(matrixp(mat)) then
       return("Non matrice"),
    [mat,[1,1,"(tm11)",dimensione(metrica[2])]])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tmat11(mat,metrica)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Crea un tensore di ordine 2, misto controvariante,
covariante, data una matrice e la metrica.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
tmat12(mat,metrica):=block([],
    if not(matrixp(mat)) then
       return("Non matrice"),
    [mat,[1,2,"(tm12)",dimensione(metrica[2])]])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tmat12(mat,metrica)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Crea un tensore di ordine 2, misto covariante,
controvariante, data una matrice e la metrica.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
tmat21(mat,metrica):=block([],
    if not(matrixp(mat)) then
       return("Non matrice"),
    [mat,[2,1,"(tm11)",dimensione(metrica[2])]])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tmat21(mat,metrica)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Crea un tensore di ordine 2, totalmente covariante,
data una matrice e la metrica.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
tmat22(mat,metrica):=block([],
    if not(matrixp(mat)) then
       return("Non matrice"),
    [mat,[2,2,"(tm22)",dimensione(metrica[2])]])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tmat22(mat,metrica)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Operazioni di somma, differenza e
prodotto per scalare, tra tensori. 
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Somma di due tensori con uguale dimensione, ordine
e tipo di indici.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
tsomma(tena,tenb):=block(
   [oa,ob,na,nb,ina,inb,stop],
   if not(tensorp(tena)) then
      return("Errore: primo arg. non tensore"),
   if not(tensorp(tenb)) then
      return("Errore: secondo arg. non tensore"),
   oa:ordine(tena), ob:ordine(tenb),
   if oa#ob then
      return("Errore: tensori di ordine diverso"),
   na:dimensione(tena),  nb:dimensione(tenb),
   if na#nb then
      return("Errore: tensori di dimensione diversa"),
   ina:indici(tena), inb:indici(tenb),
   stop:false,
   for j:1 thru oa do if mod(ina[j]+inb[j],2)=1 
       then stop:true,
   if stop then 
      return("Errore: indici incongruenti"),
   append(ratsimp(rest(tena,-1)+rest(tenb,-1)),
      [append(ina,["(Somma tensori)",na])])
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tsomma(tena,tenb)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Differenza tra due tensori con uguale dimensione, ordine
e tipo di indici.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
tmeno(tena,tenb):=block(
   [oa,ob,na,nb,ina,inb,stop],
   if not(tensorp(tena)) then
      return("Errore: primo arg. non tensore"),
   if not(tensorp(tenb)) then
      return("Errore: secondo arg. non tensore"),
   oa:ordine(tena), ob:ordine(tenb),
   if oa#ob then
      return("Errore: tensori di ordine diverso"),
   na:dimensione(tena),  nb:dimensione(tenb),
   if na#nb then
      return("Errore: tensori di dimensione diversa"),
   ina:indici(tena), inb:indici(tenb),
   stop:false,
   for j:1 thru oa do if mod(ina[j]+inb[j],2)=1 
       then stop:true,
   if stop then 
      return("Errore: indici incongruenti"),
   append(ratsimp(rest(tena,-1)-rest(tenb,-1)),
      [append(ina,["(Differenza tra tensori)",na])])
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tmeno(tena,tenb)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Prodotto di un tensore scalare per un tensore della 
stessa dimensione.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
tprod(tscala,tenb):=block(
   [oa,ob,na,nb,inb],
   if not(tensorp(tscala)) then
      return("Errore: primo arg. non tensore"),
   if not(tensorp(tenb)) then
      return("Errore: secondo arg. non tensore"),
   oa:ordine(tscala), ob:ordine(tenb),
   if oa#0 then
      return("Errore: non scalare"),
   na:dimensione(tscala),  nb:dimensione(tenb),
   if na#nb then
      return("Errore: tensori di dimensione diversa"),
   inb:indici(tenb),
   append(ratsimp(tscala[1]*rest(tenb,-1)),[append(indici(tenb),
      ["(Prodotto scalare*tensore)",nb])])
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tprod(tscala,tenb)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
tassegna e tvale : per assegnare
o  conoscere valori di singoli
elementi.
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Per non dovere specificare gli indici accoppiandoli.
Assegna un dato ad un elemento di matrice di 
matrice oppure trova il valore di quell'elemento.
Sono funzioni ricorsive:
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
assegnam(mt,o,val,p):=block([],
    if o=length(p) then (
       mt[p[o-1],p[o]]:val,
       return(true)),
    assegnam(mt[p[o-1],p[o]],o+2,val,p)
    )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
valem(miot,o,p):=block([],
    if o=length(p) then return(miot[p[o-1],p[o]]),
    valem(miot[p[o-1],p[o]],o+2,p)
    )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Assegna un valore ad un elemento del tensore oppure
trova il valore di quell'elemento.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
tassegna(tn,val,p):=block([o],
    o:length(tn[length(tn)])-2,
    if length(p)#o then 
    return("errore: incongruenza di indici"),
    if o=0 then ( tn[1]:val,return(true)),
    if o=1 then ( tn[p[1]]:val, return(true)),
    if mod(o,2)=0 then assegnam(tn[1],2,val,p)
    else assegnam(tn[p[1]],2,val,rest(p))
    )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
tvale(tn,p):=block([o],
    o:length(tn[length(tn)])-2,
    if length(p)#o then 
    return("errore: incongruenza di indici"),
    if o=0 then return(tn[1]),
    if o=1 then return(tn[p[1]]),
    if mod(o,2)=0 then valem(tn[1],2,p)
    else ( 
        valem(tn[p[1]],2,rest(p)))
    )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tassegna(tn,val,p)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tvale(tn,p)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
permuta : permutazione di indici
di un tensore
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
La funzione che permuta gli indici di un qualsiasi 
tensore creando un altro tensore con gli indici permutati
in base all'ordinamento specificato dalla lista per.
La lista per deve indicare da dove è stato preso il dato
che compare in quella posizione.
Per esempio [3,1,2,4] indica che il dato che compare
come primo indice è quello che prima era situato in 
colonna 3, come secondo dato quello che prima stava in
colonna 1 etc.
La rotazione degli indici sa fa con [2,3,4,1].
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
permuta(tn,per,comm):=block([tr,u,o,a,b,j,nd,va,vb,id],
    if not(tensorp(tn)) 
       then return("Errore, non tensore"),
    if not(listp(per)) then 
       return("Errore, non lista permutazioni"),
    u:length(tn),o:length(tn[u])-2,
    if 2>o then return("Applicabile da ordine 2 in su"),
    if o#length(per) then
       return("Errore, lista permutazioni incongruente"),
    j:1,for a:1 thru o-1 do for b:a+1 thru o do
        if per[a]=per[b] then j:0,
    if j=0 then return("Permutazione errata"),
    nd:tn[u][o+2],
    va:makelist(0,j,1,o),
    vb:makelist(0,j,1,o),
    id:makelist(0,j,1,o),
    for j:1 thru o do id[j]:tn[u][per[j]],
    tr:zerotensor(id,comm,nd),
    funfor(tr,tn,va,vb,per,1,o,nd),
    tr
    )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Faccio un ciclo for in forma ricorsiva...
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
funfor(ta,tb,va,vb,per,n,m,nd):=block([j,np],
   if n=m then (
      for j:1 thru nd do (
         va[n]:j,
         vb[per[n]]:j,
         tassegna(ta,tvale(tb,vb),va) ))
   else (
      np:n+1,
      for j:1 thru nd do (
          va[n]:j,
          vb[per[n]]:j,
          funfor(ta,tb,va,vb,per,np,m,nd))
         )
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("permuta(tn,per,com)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Funzione per trovare come è permutata una lista ossia
data una lista e una permutazione, cosa diventa la
lista permutata.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
ada(la,per):=block([n,j,k,lr],
    n:length(per),
    if n#length(la) then return("Disuguale lunghezza liste"),
    for j:1 thru n-1 do for k:j+1 thru n do if
       per[j]=per[k] then return("Indici duplicati"),
    lr:makelist(0,j,1,n),
    for j:1 thru n do lr[j]:la[per[j]],
    lr)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("ada(la,per)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
tensoredalista : per creare un tensore
usando una lista per inizializzarlo
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Crea un tensore del tipo voluto, prelevando i dati da una lista
che se è troppo breve viene riutilizzata ciclicamente e
naturalmente, se è troppo lunga, viene usata solo parzialmente.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
forlista(tb,vb,n,m,nd,lista,dove):=block([j,np],
   if n=m then (
      for j:1 thru nd do (
         vb[n]:j,
         dove[1]:mod(dove[1],length(lista))+1,
         tassegna(tb,lista[dove[1]],vb) ))
   else (
      np:n+1,
      for j:1 thru nd do (
          vb[n]:j,
          forlista(tb,vb,np,m,nd,lista,dove))
         )
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
tensoredalista(lista,tipi,comme,ndim):=block(
    [tr,dove,m,vb,j],
    tr:zerotensor(tipi,comme,ndim),
    dove[1]:0,
    m:length(tipi),
    vb:makelist(0,j,1,m),
    forlista(tr,vb,1,m,ndim,lista,dove),
    tr)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("tensoredalista(lista,tipi,comme,ndim)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
listadatensore : per estrarre gli elementi
di un tensore mettendoli in una lista
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Estrae gli elementi di un tensore producendo una lista.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
forinlista(lista,tb,vb,n,m,nd,dove):=block([j,np],
   if n=m then (
      for j:1 thru nd do (
         vb[n]:j,
         dove[1]:mod(dove[1],length(lista))+1,
         lista[dove[1]]:tvale(tb,vb)))
   else (
      np:n+1,
      for j:1 thru nd do (
         vb[n]:j,
         forinlista(lista,tb,vb,np,m,nd,dove))
         )
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
listadatensore(tn):=block([lista,u,o,dove,vb,j],
    u:length(tn),o:length(tn[u])-2,
    nd:tn[u][o+2],
    dove[1]:0,
    vb:makelist(0,j,1,o),
    lista:makelist(0,j,1,nd^o),
    forinlista(lista,tn,vb,1,o,nd,dove),
    lista)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("listadatensore(tn)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
traccia : generalizzazione del calcolo
della traccia di una matrice.
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
traccia di un tensore, generalizzazione di quella di una 
matrice...
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
fortraccia(tr,vr,tn,vn,h,k,n,nn,m,nd):=block([j,np,nnp,ss],
   if n=m then (
      ss:0,
      for j:1 thru nd do (
         vn[h]:j,vn[k]:j,
         ss:ss+tvale(tn,vn)),
         tassegna(tr,ss,vr))
   else (
      np:n+1, nnp:nn+1,
      if nnp=h then nnp:nnp+1,
      if nnp=k then nnp:nnp+1,
      for j:1 thru nd do (
         vr[np]:j,
         vn[nnp]:j,
         fortraccia(tr,vr,tn,vn,h,k,np,nnp,m,nd))
         )
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
traccia(tn,ih,ik):=block([u,o,nd,vr,vb,h,k,tipi,i,j],
    if not(tensorp(tn)) then
       return("Errore: non tensore"),
    u:length(tn),o:length(tn[u])-2,nd:tn[u][o+2],
    if 2>o then return("Errore, ordine inferiore a 2"),
    if ih=ik then return("Errore, indici uguali"),
    h:min(ih,ik),k:max(ih,ik),
    if k>o then 
       return(["Errore, indice superiore ad ordine",k,o]),
    vb:makelist(0,j,1,o),vr:[],
    if o>2 then vr:makelist(0,j,1,o-2),
    tipi:indici(tn),
    if tipi[h]=tipi[k] then 
        print("Attenzione: indici dello stesso tipo"),
    tr:zerotensor(vr,"ridotto",nd),
    i:0,for j:1 thru o do (
        if j#h then if j#k then (
           i:i+1,
           tr[u][i]:tipi[j])),
    fortraccia(tr,vr,tn,vb,h,k,0,0,o-2,nd),
    tr)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("traccia(tn,ih,ik)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
scala : generalizzazione del prodotto
scalare tra vettori e tra matrici. 
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Prodotto scalare generalizzato. Attenzione... è posizionale
e per ottenere gli indici nell'ordine voluto bisogna poi
riordinare il risultato con la funzione permuta.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
forscala(tr,vr,ta,va,tb,vb,h,k,n,nn,ma,m,nd):=block(
   [j,np,nnp,ss],
   if n=m then (
      va[h]:1,
      vb[k]:1,
      ss:0,
      for j:1 thru nd do (
         va[h]:j,vb[k]:j,
         ss:ss+tvale(ta,va)*tvale(tb,vb)),
         tassegna(tr,ss,vr))
   else (
      np:n+1, nnp:nn+1,
      if nnp=h then nnp:nnp+1,
      if nnp=k+ma then nnp:nnp+1,
      for j:1 thru nd do (
         vr[np]:j,
         if nnp>ma then vb[nnp-ma]:j
         else va[nnp]:j,
         forscala(tr,vr,ta,va,tb,vb,h,k,np,nnp,ma,m,nd))
         )
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
scala(ta,ia,tb,ib):=block(
    [ua,oa,ub,ob,ur,ot,nd,vr,va,vb,tipa,tipb,i,j],
    if not(tensorp(ta)) then
       return("Errore: ta non tensore"),
    if not(tensorp(tb)) then
       return("Errore: tb non tensore"),
    ua:length(ta),oa:length(ta[ua])-2,nd:ta[ua][oa+2],
    if 1>oa then return("Errore, ordine ta inferiore a 1"),
    ub:length(tb),ob:length(tb[ub])-2,
    if nd#tb[ub][ob+2] then
       return("Errore, dimensioni diverse"),
    if 1>ob then return("Errore, ordine tb inferiore a 1"),
    if ia>oa then 
       return(["Errore, indice superiore ad ordine",ia,oa]),
    if ib>ob then 
       return(["Errore, indice superiore ad ordine",ib,ob]),
    va:makelist(0,j,1,oa),vb:makelist(0,j,1,ob),
    ot:oa+ob-2,
    vr:[],if ot>0 then vr:makelist(0,j,1,ot),
    tipa:indici(ta),tipb:indici(tb),
    if tipa[ia]=tipb[ib] then 
        print("Attenzione: indici dello stesso tipo"),
    tr:zerotensor(vr,"prod-scalare",nd),
    ur:length(tr),
    i:0,for j:1 thru oa+ob do (
       if j>oa then ( if (j-oa)#ib then (
          i:i+1,
          tr[ur][i]:tipb[j-oa]))
       else if j#ia then (
          i:i+1, 
          tr[ur][i]:tipa[j])),
    forscala(tr,vr,ta,va,tb,vb,ia,ib,0,0,oa,ot,nd),
    tr)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("scala(ta,ia,tb,ib)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Funzioni per il calcolo differenziale
   [wxMaxima: section end   ] */

/* [wxMaxima: subsect start ]
fa_diff : derivata ordinaria di un tensore
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Richiede l'uso di una permutazione finale perché spontaneamente
aggiungerebbe l'indice 0 all'inizio e non alla fine.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
fa_diff(tn,variabili,comme):=block([tn0,tnn,indi,o,j],
    tnn:rest(tn,-1),
    o:ordine(tn),
    indi:indici(tn),
    nd:dimensione(tn),
    if mod(o,2)=0 then (
       tn0:makelist(0,j,1,nd+1),
       for j:1 thru nd do 
       tn0[j]:ratsimp(diff(tnn[1],variabili[j])),
       tn0[nd+1]:append(cons(0,indi),[comme,nd])
       )
    else (
       tn0:[0,append(cons(0,indi),[comme,nd])],
       tn0[1]:ratsimp(matrix(diff(tnn,variabili[1]))),
       for j:2 thru nd do
       tn0[1]:addrow(tn0[1],
              ratsimp(matrix(diff(tnn,variabili[j]))))
       ),
    if o>0 then 
       permuta(tn0,makelist(mod(j,o+1)+1,j,1,o+1),comme)
    else tn0
    )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("fa_diff(tn,variabili,comme)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
fa_ch222 : Simboli di Christoffel di prima specie 
   [wxMaxima: subsect end   ] */

/* [wxMaxima: input   start ] */
fa_ch222(g220):=block([g022,g202],
     if controllo then print("Inizia fa_ch222"),
     g022:permuta(g220,[2,3,1],"-"),
     g202:permuta(g022,[2,3,1],"-"),
     append((+rest(g220,-1)+rest(g022,-1)-rest(g202,-1))/2,
     [[2,2,2,"Christoffel prima specie",dimensione(g220)]])
     )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("fa_ch222(g220)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
fa_ch122 : Per fare i simboli di Christoffel di 
seconda specie
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Eleva il primo indice.. e fa il simbolo di Christoffel
di seconda specie. Anche in questo caso il secondo e 
il terzo indice sono scambiabili.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
fa_ch122(g22,listavariabili):=block([g11,g220],
   if controllo then print("Inizia fa_ch222"),
   g11: [ratsimp(invert(g22[1])),
        [1,1,"-",dimensione(g22)]],
   g220: fa_diff(g22,listavariabili,"g220"),
   ratsimp(scala(g11,2,fa_ch222(g220),1))
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("fa_ch122(g22,listavariabili)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
fa_r1222 : genera il tensore di Riemann 
col primo indice controvariante.
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Ora calcolo il tensore di Riemann con la formula dell'HEL 
Hobson,Efstathiou,Lasenby ISBN 9780521829519 a pag. 158.
Ricopio la formula dal libro: 
r1d2a2b2c = ch1d2a2a0b - 
            ch1d2a2b0c + 
            ch1e2a2c*ch1d2e2b - 
            ch1e2a2b*ch1d2e2c
Il nome degli indici però non mi piace e preferisco usare 
come nomi a,b,c,d sostituti di d,a,b,c ed inoltre usare 
il nome j al posto del nome e, per cui la formula diventa:
r1a2b2c2d = ch1a2b2d0c - 
            ch1a2b2c0d + 
            ch1j2b2d*ch1a2j2c - 
            ch1j2b2c*ch1a2j2d
In linguaggio Maximese non interpongo i nomi degli indici 
alle cifre indicatrici di tipo ma scrivo i nomi tra 
parentesi quadre:
r1222[a,b][c,d]:ch1220[a,b][d,c]-
                ch1220[a,b][c,d]+
                ch122[j][b,d]*ch122[a][j,c]-
                ch122[j][b,c]*ch122[a][j,d]
Per fare il calcolo bisogna disporre sia del simbolo di 
Christoffel di seconda specie che dello pseudotensore 
delle sue derivate.
Spezzo il calcolo in due parti in modo che sia possibile
partire o dal tensore metrico covariante o dai simboli
di Christoffel di seconda specie.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
fa_r1222(t22,listavariabili):=block(
   [nd,g22,g11,g220,g022,g202,ch222,ch122],
   if controllo then print("Inizia fa_r1222"),
   if not(tensorp(t22)) then 
      return("Primo argomento non tensore"),
   nd:dimensione(t22),
   g22:[ratsimp((t22[1]+transpose(t22[1]))/2),[2,2,"(g22)",nd]],
   g11:[ratsimp(invert(g22[1])),[1,1,"(g11)",nd]],
   g220:fa_diff(g22,listavariabili,"(g220)"),
   g022:permuta(g220,[2,3,1],"-"),
   g202:permuta(g022,[2,3,1],"-"),
   ch222:append((rest(g220,-1)+rest(g022,-1)-rest(g202,-1))/2,
         [[2,2,2,"Christoffel prima specie",dimensione(g220)]]),
   ch122:scala(g11,2,ch222,1),
   fa_ch_r1222(ch122,listavariabili))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
fa_ch_r1222(ch122,listavariabili):=block(
   [r1222,ch1a2b2c2d,ch1a2b2d2c,xabcd,xacbd,xadbc],
   if controllo then print("Inizia fa_ch_r1222"),
   ch1a2b2c2d:fa_diff(ch122,listavariabili,"ch1a2b2c2d"),
   ch1a2b2d2c:permuta(ch1a2b2c2d,[1,2,4,3],"ch1a2b2d2c"),
   xabcd:scala(ch122,3,ch122,1),
   xacbd:permuta(xabcd,[1,3,2,4],"xacbd"),
   xadbc:permuta(xabcd,[1,4,3,2],"xadbc"),
   r1222:[ch1a2b2d2c[1]-ch1a2b2c2d[1]+ xacbd[1] - xadbc[1],
         [1,2,2,2,"(r1222)",dimensione(ch122)]]
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("fa_r1222(t22,listavariabili)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
VERIFICA:
Le funzioni della libreria alternativa, usabile 
per calcolare il tensore di Riemann r1222 sono queste
che non sono ricorsive e non gestiscono gli indici.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
der2(t2,variabili):=block([t3,nd,j],
   nd:length(variabili),
   t3:makelist(0,j,1,nd),
   for j:1 thru nd do
       t3[j]:diff(t2,variabili[j]),
   t3)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
chris222(g22,variabili):=block([g220,ss,tn],
   g220:r3(der2(g22,variabili)),
   ss:r3(g220),
   tn:g220+ss,
   ss:r3(ss),
   ratsimp((tn-ss)/2)
)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
chris122(g22,variabili):=block([g11],
g11:invert(g22),
ch222:chris222(g22,variabili),
transpose(g11.ch222)[1]
)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
r3(t3):=block([nd,j,k,h,t3n],
 nd:length(t3),
 t3n:makelist(ident(nd),j,1,nd),
 for j:1 thru nd do 
 for k:1 thru nd do
 for h:1 thru nd do
   t3n[j][k,h]:t3[h][j,k],
 t3n
)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
r4(t4):=block([nd,i,j,k,h,t4n],
nd:length(t4),
t4n:genmatrix(lambda([i,j],ident(nd)),nd,nd),
for i:1 thru nd do
for j:1 thru nd do 
for k:1 thru nd do
for h:1 thru nd do
   t4n[i,j][k,h]:t4[h,i][j,k],
t4n
)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
der3(t3,variabili):=block([t4,t4d,nd,j],
   if not(listp(t3)) then (print(["Errore",t3]),
   return("Non lista t3")),
   nd:length(variabili),
   t4d:der2(t3,variabili),
   t4:matrix(t4d[1]),
   for j:2 thru nd do
     t4:addrow(t4,t4d[j]),
   ratsimp(t4))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Bisogna chiamare la riemann_1222 passandogli il tensore metrico 
covariante e lei fa tutto...
Questa funzione è utile perché è più semplice di quella
ricorsiva ed adatta ad essere riempita di stampe di controllo...
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
riemann_1222(t22,variabili):=block(
   [g22,ch122,ch1220,r1222,nd,a,b,c,d,j,sp],
   if not(tensorp(t22)) then 
      return("Primo argomento non tensore"),
   nd:dimensione(t22),
   g22:ratsimp((t22[1]+transpose(t22[1]))/2),
   ch122:chris122(g22,variabili),
   ch1220:r4(der3(ch122,variabili)),
   r1222:genmatrix(lambda([i,j],ident(nd)),nd,nd),
   for a:1 thru nd do
   for b:1 thru nd do
   for c:1 thru nd do
   for d:1 thru nd do (
      sp:ch1220[a,b][d,c]-ch1220[a,b][c,d],
      for j:1 thru nd do sp:sp+
           ch122[j][b,d]*ch122[a][j,c]-
           ch122[j][b,c]*ch122[a][j,d],
      r1222[a,b][c,d]:sp),  
   [ ratsimp(r1222),[1,2,2,2,"(r1222)",nd]]    
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("riemann_1222(t22,variabili)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
fa_r2222 : genera il tensore di Riemann
totalmente covariante.
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Calcolo il tensore di Riemann in forma totalmente 
covariante, in base a quanto sta scritto sull' HEL.

r2a2b2c2d = (g2b2c0d0a-
             g2a2c0d0b+
             g2a2d0b0c-
             g2b2d0a0c )/2 +
             ch1j2b2c*ch2j2a2d - 
             ch1j2b2d*ch2j2a2c 

ovvero in Maximese

r2222[a,b][c,d] = ( g2200[b,c][a,d] - 
                    g2200[a,c][b,d] + 
                    g2200[a,d][b,c] - 
                    g2200[b,d][a,c] )/2 +
                    ch122[j][a,d]*ch222[j][b,c] - 
                    ch122[j][a,c]*ch222[j][b,d]

Il tensore di Riemann in forma totalmente covariante
evidenzia le sue simmetrie. E' antisimmetrico scambiando
gli indici della prima coppia ed e' antisimmetrico
scambiando gli indici della seconda coppia mentre
è simmetrico scambiando le coppie tra loro.
Per queste simmetrie il calcolo del tensore di Riemann
potrebbe essere molto più veloce evitando inutili ricalcoli
ma qui ...non si bada a spese...
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
fa_r2222(t22,listavariabili):=block(
   [nd,g22,xg11,g11,g220,g022,g202,ch222,ch122,g2200],
   if controllo then print("Inizia fa_r2222"),
   if not(tensorp(t22)) then 
      return("Primo argomento non tensore"),
   nd:dimensione(t22),
   g22:[ratsimp((t22[1]+transpose(t22[1]))/2),[2,2,"(g22)",nd]],
   xg11:trigsimp(ratsimp(invert(g22[1]))),
   g11:[xg11,[1,1,"(g11)",nd]],
   g220:fa_diff(g22,listavariabili,"(g220)"),
   g022:permuta(g220,[2,3,1],"-"),
   g202:permuta(g022,[2,3,1],"-"),
   ch222:append(ratsimp((rest(g220,-1)+
         rest(g022,-1)-rest(g202,-1))/2),
         [[2,2,2,"Christoffel prima specie",dimensione(g220)]]),
   ch122:scala(g11,2,ch222,1),
   g2200:fa_diff(g220,listavariabili,"(yabcd) ovvero g2200"),
   fa_poi_r2222(ch222,ch122,g2200)
        )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
fa_poi_r2222(ch222,ch122,g2200):=block(
   [zabcd,zadbc,zacbd,ybcad,yacbd,yadbc,ybdac],
   if controllo then print("Inizia fa_poi_r2222"),
   zabcd:scala(ch122,1,ch222,1),
   zadbc:permuta(zabcd,[1,3,4,2],"(zadbc)"),
   zacbd:permuta(zabcd,[1,3,2,4],"(zacbd)"),
   ybcad:permuta(g2200,[3,1,2,4],"(ybcad)"),
   yacbd:permuta(g2200,[1,3,2,4],"(yacbd)"),
   yadbc:permuta(g2200,[1,3,4,2],"(yadbc)"),
   ybdac:permuta(g2200,[3,1,4,2],"(ybdac)"),
   [ratsimp((ybcad[1]-yacbd[1]+yadbc[1]-ybdac[1])/2 
             +zadbc[1]-zacbd[1]),
     [2,2,2,2,"(r2222)",nd]] )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("fa_r2222(t22,listavariabili)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Ecco la funzione tratta della altra libreria... 
Sembra piu' veloce essendo meno generale dato che non fa uso
di funzioni ricorsive ma utilizza il costrutto for.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
riemann_2222(t22,variabili):=block(
   [ss,tn,nd,g22,ch222,ch122,
   g2200,r2222,a,b,c,d,j],
   if not(tensorp(t22)) then 
      return("Primo argomento non tensore"),
   nd:dimensione(t22),
   g22:ratsimp((t22[1]+transpose(t22[1]))/2),
   g220:r3(der2(g22,variabili)),
   ss:r3(g220),
   tn:g220+ss,
   ss:r3(ss),
   ch222:ratsimp((tn-ss)/2),
   ch122:transpose(invert(g22).ch222)[1],
   g2200:r4(der3(g220,variabili)),
   r2222:genmatrix(lambda([i,j],ident(nd)),nd,nd),
   for a:1 thru nd do
   for b:1 thru nd do
   for c:1 thru nd do
   for d:1 thru nd do (
      ss:(g2200[b,c][a,d] - g2200[a,c][b,d] + 
          g2200[a,d][b,c] - g2200[b,d][a,c] )/2,
      for j:1 thru nd do ss:ss+
         ch122[j][a,d]*ch222[j][b,c] - 
         ch122[j][a,c]*ch222[j][b,d], 
      r2222[a,b][c,d]:ss),  
   [ ratsimp(r2222),[2,2,2,2,"(r2222)",nd]]   
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("riemann_2222(t22,variabili)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
fa_r22 : genera il tensore di Ricci 
totalmente covariante.
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Definisco il tensore di Ricci doppiamente covariante
seguendo ancora la definizione dell' HEL.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
fa_r22(t22,listavariabili):=block([],
    if controllo then print("Inizia fa_r22"),
    traccia(fa_r1222(t22,listavariabili),1,4)
    )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("fa_r22(t22,listavariabili)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Utilizza la libreria alternativa non ricorsiva.
Dal tensore metrico covariante mi fornisce subito
il tensore di Ricci covariante.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
ricci_22(t22,variabili):=block(
   [nd,r1222,ric22,a,b,sp,j],
    if not(tensorp(t22)) then 
       return("Primo argomento non tensore"),
    nd:dimensione(t22),
    r1222:riemann_1222(t22,variabili)[1],
    r22:zeromatrix(nd,nd),
    for a:1 thru nd do
    for b:a thru nd do (
        sp:0,
        for j:1 thru nd do sp:sp+r1222[j,a][b,j],
        r22[a,b]:sp,
        r22[b,a]:sp
        ),
    [ ratsimp(r22),[2,2,"(Ricci 22)",nd ] ]
    )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
fametrica: costruttore di metrica. Genera 
una lista che include tutto quello che 
serve per fare derivate covarianti.
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Il calcolo differenziale tensoriale si basa su una
funzione, fametrica, che genera la metrica usata per
calcolare la derivata covariante.
Questa funzione deve ricevere in argomento una matrice
che, per sicurezza, viene simmetrizzata e usata come
tensore metrico covariante della metrica. 
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
fa_metrica(mg,listavariabili):=block(
   [metrica,m2,m3,m4,m4,m6,m7,m8,m9,ma],
   if not(matrixp(mg)) then 
      return("Errore, primo arg. non matrice"),
   if not(listp(listavariabili)) then
      return("Errore, secondo arg. non lista variabili"),
   if length(listavariabili)#length(mg) then
      return("Errore: incongruenza tra matrice e lista variabili"),
   if controllo then print("Inizia g22"),
   m2: [ratsimp((mg + transpose(mg) )/2),
       [2,2,"(tensore metrico covariante)",length(mg)]],
   m2: trigsimp(m2),
   if controllo then print("Inizia g11"),
   m3: [ratsimp(invert(m2[1])),
       [1,1,"(tensore metrico controvariante)",length(mg)]],
   if controllo then print("Inizia g220"),
   m4: fa_diff(m2,listavariabili,"g220"),
   if controllo then print("Inizia g2200"),
   m5: fa_diff(m4,listavariabili,"g2200"),
   m6: fa_ch222(m4),
   commentalo(m6,"(Christoffel prima specie)"),
   if controllo then print("Inizia ch122"),
   m7: scala(m3,2,m6,1),
   commentalo(m7,"(Christoffel seconda specie)"),
   m8: fa_r1222(m2,listavariabili),
   m9: fa_r2222(m2,listavariabili),
   if controllo then print("Inizia r22"),
   ma: traccia(m8,1,4),
   commentalo(ma,"(Ricci totalmente covariante)"),
   mb:[ratsimp(sqrt(abs(determinant(m2[1])))),
      ["sqrt(|determinant(g22)|)",length(mg)]],
   metrica:[listavariabili,m2,m3,m4,m5,m6,m7,m8,m9,ma,mb],
   infometrica(),
   metrica)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
infometrica():=block([],
   print("Ha creato un vettore di undici elementi:"),
   print("[1] : lista delle variabili"),
   print("[2] : tensore metrico covariante"),
   print("[3] : tensore metrico controvariante"),
   print("[4] : derivate prime del tens. metrico cov."),
   print("[5] : derivate seconde del tens. m. cov."),
   print("[6] : simboli di Christoffel di prima specie"),
   print("[7] : simboli di Christoffel di seconda specie"),
   print("[8] : Riemann col solo primo indice controvariante"),
   print("[9] : Riemann totalmente controvariante"),
   print("[10]: Ricci totalmente covariante"),
   print("[11]: sqrt(abs(determinate del tens.metr.cov.))")
   )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("fa_metrica(mg,listavariabili)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Derivata covariante
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Per convenzione la cifra che segnala un indice ottenuto da
derivazione covariante è 8 e quello trasformato in 
controvariante è 9. Dunque xx28 è un vettore covariante ( xx2 )
diventato tensore del secondo ordine per derivazione covariante
ed analogamente xx18 è un vettore controvariante ( xx1), 
derivato covariantemente. Se alzo l'ultimo indice di xx28 
ottengo xx29, etc...
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Dato che la derivazione covariante produce tensori,
è possibile derivare qualsiasi tensore usando i simboli
di Christoffel in due modi diversi a seconda della
natura covariante o controvariante di un dato
indice.
Per fare la derivazione occorrono tanti simboli di Christoffel
quanti sono gli indici del tensore e dunque per derivare 
uno scalare non occorrono simboli di Christoffel ossia
la derivazione di uno scalare coincide con quella ordinaria.
Esaminiamo i due casi distinti usando le convenzioni della
presente libreria che impone che i tensori di qualunque
ordine siano delle liste di qualcosa, scalari o matrici
e che l'ultimo elemento della lista sia a sua volta una
lista che contiene informazioni sulla natura degli 
indici del tensore.

Dunque:

A[a] è uno scalare o un vettore del primo ordine. Se è
uno scalare esiste solo A[1] mentre se è un vettore
esistono anche le componenti A[2]...A[nd] dove nd
è il numero di dimensioni dello spazio, ovvero 3
in meccanica classica e 4 in meccanica relativistica.
B[a][b,c] è un tensore del secondo o del terzo ordine.
C[a][b,c][d,e] è un tensore del quarto o del quinto ordine.

Per accedere alla lista dei tipi degli indici
bisogna prendere l'ultimo elemento del primo indice ossia
A[length(A)] o B[length(B)] o C[length[C)] è la
lista dei tipi degli indici, qualunque sia l'ordine
del tensore A o B o C.
L'ultima delle componenti, se il tensore ha ordine pari 
è sempre 2 mentre se ha ordine dispari è sempre
data da nd+1.
In pratica però è meglio non fare affidamento su
questa regola e prendere sempre l'ultimo indice
della lista perché in questo modo è possibile
usare anche liste più lunghe del minimo indispensabile
e memorizzare negli elementi sovrabbondanti i
dati che si desidera associare in più al 
tensore.
Notare che in Maxima è ammesso omettere gli indici
tra parentesi quadra per cui se scrivo B[1] ottengo 
l'unica o la prima matrice del tensore B ( o un singolo
dato numerico se B è uno scalare o un vettore); la
matrice ( o scalare) è l'unica se B è un tensore di 
ordine pari ossia ad esempio uno scalare o un tensore 
del secondo ordine ed è la prima se B è un tensore di 
ordine dispari ossia, per esempio un vettore o un 
tensore del terzo ordine.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Per la derivata di un vettore covariante A2, ossia di un 
tensore di ordine 1 dove la cifra 2 ricorda che il suo
primo ed unico indice è covariante, la formula è questa:

   A2a8b = A2a0b - ch1j2a2b*A2j

ovvero in Maximese ( essendo 8 la cifra che indica 
derivata covariante ) :

   A28[1][a,b] = A20[1][a,b] - ch122[j][a,b]*A2[j]

dove la cifra 0 specifica la derivata ordinaria e ch122 
rappresenta il simbolo di Christoffel di seconda specie 
avente la cifra 1 che ricorda che il suo primo indice è
(pseudo)controvariante e la cifra 2 che ricorda che
il secondo e terzo indice sono (pseudo)covarianti.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Per la derivata di un tensore del secondo ordine totalmente
covariante A22, la formula è questa:

   A2a2b8c = A2a2b0c - ch1j2a2c*A2j2b - 
                       ch1j2b2c*A2a2j

ovvero in Maximese, dato il tensore A[1][a,b] ottengo:

   A228[a][b,c] = A220[a][b,c] - 
                      ch122[j][a,c]*A22[1][j,b] -
                      ch122[j][b,c]*A22[1][a,j]

per cui è facile intuire la regola generale da applicarsi
ad un tensore di ordine qualsiasi. 
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
La derivata di un tensore del terzo ordine totalmente
covariante ossia A222, la formula è questa:

   A2a2b2c2d = A2a2b2c0d - ch1j2a2d*A2j2b2c -
                           ch1j2b2d*A2a2j2c -
                           ch1j2c2d*A2a2b2j

ovvero in Maximese, dato il tensore A[a][b,c]:

   A2228[1][a,b][c,d] = A2220[1][a,b][c,d] - 
                        ch122[j][a,d]*A222[j][b,c] -
                        ch122[j][b,d]*A222[a][j,c] -
                        ch122[j][c,d]*A222[a][b,j]

A questo punto è facile generalizzare il metodo
a tensori di qualsiasi ordine.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
La derivata di un vettore controvariante si fa invece
in questo modo ( usando uno dei due indici 
pseudocovarianti del simbolo di Christoffel che è 
uno pseudo tensore) :

   A1a8b = A1a0b + ch1a2j2b*A1j

Notare che il simbolo di Christoffel è simmetrico nei 
suoi due ultimi indici trattati come covarianti, per cui
avrei potuto scrivere anche in questo modo ( è solo una 
questione estetica )

   A1a8b = A1a0b + ch1a2b2j*A1j

ossia in Maximese, detto A1[a] il vettore da derivare
ho la seguente formula

  A18[1][a,b] = A10[1][a,b] + ch122[a][j,b]*A1[j]

   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
La derivata di un tensore del secondo ordine totalmente
controvariante A1a1b si fa, di conseguenza in questo modo:

   A1a1b8c = A1a1b0c + ch1a2j2c*A1j1b
                     + ch1b2j2c*A1a1j

ossia in Maximese, dato A11[1][a,b] si ha:

   A118[a][b,c] = A110[a][b,c] +
                     ch122[a][j,c]*A11[1][j,b] +
                     ch122[b][j,c]*A11[1][a,j]

Attenzione: il primo contributo è dato da A11[1][j,b]
che non è equivalente a A11[1][b,j] a meno che il 
tensore del secondo ordine non sia, in pratica,
una matrice simmetrica.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
diffcov(tn,metrica,nome):=block([indi,o,ch122,tr,lp,j,temp],
    if not(tensorp(tn)) 
       then return("Errore: non tensore"),
    indi:indici(tn),
    o:ordine(tn),
    if dimensione(tn)#dimensione(metrica[2]) then
       return("Errore, dimensioni incongruenti"),
    ch122:metrica[7],
    tr:rest(fa_diff(tn,metrica[1],"-"),-1),
    lp:makelist(mod(j,o+1)+1,j,1,o+1),
    lp[1]:1,lp[o+1]:2,
    for j:1 thru o do (
        if mod(indi[j],2)=0 then (
            tr:tr-rest(permuta(scala(ch122,1,tn,j),lp,o),-1)
            )
        else (
            tr:tr+rest(permuta(scala(ch122,2,tn,j),lp,o),-1)
            ),
        temp:lp[j],lp[j]:lp[j+1],lp[j+1]:temp
        ),
    endcons(append(indi,
       [8,nome,dimensione(tn)]),ratsimp(tr))
    )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("diffcov(tn,metrica,nome)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Per mettere alla dura prova questa libreria proviamo a
verificare alcune identità che debbono sempre sussistere.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
A pag. 344, formula (92,4), il Landau afferma.
La somma ciclica di qualsiasi terna di indici del tensore
di Riemann r2222 deve essere sempre nulla.
In altre parole :

    r2222[1][a,b][c,d] + 
    r2222[1][a,d][b,c] + r2222[1][a,c][d,b] = 0

Supponiamo di aver calcolato r2222....
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Usando la funzione permuta si ottiene in generale un tensore 
del quarto ordine che, a quanto dice il Landau ma non solo
lui, dovrebbe risultare totalmente nullo.
Con le funzioni di questa libreria si dovrebbe fare così:

rest(r2222,-1)+
rest(permuta(r2222,[1,3,4,2],""),-1)+
rest(permuta(r2222,[1,4,2,3],""),-1);


   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Per provare la "identità di Bianchi" ossia la seguente
relazione:

    r12228[a][b,c][d,e] +
    r12228[a][b,e][c,d] +
    r12228[a][b,d][e,c] = 0

bisogna fare la derivata covariante del tensore di Riemann 
col solo primo indice controvariante ossia prelevarlo dalla
metrica.
Con questa libreria si fa così:

r12228:ratsimp(diffcov(r1222,metrica,"r12228"))$

Notare che è un calcolo lungo perché si tratta di calcolare
un tensore del quinto ordine ossia con 4^5 = 1024 componenti!
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Noto r12228, ecco come va calcolata l'identità di Bianchi con 
la libreria ossia permutando gli indici. In quattro dimensioni 
questo tensore del quinto ordine ha 1024 componenti che 
dovrebbero risultare tutte nulle!

idbianchi:rest(r12228,-1)+
          rest(permuta(r12228,[1,2,4,5,3],""),-1)+
          rest(permuta(r12228,[1,2,5,3,4],""),-1)$

Stupefacente! Vengono tutti zeri come deve essere...
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Mi sembra che questo risultato sia una buona dimostrazione
della correttezza della libreria...
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Ora uso l'invariante del tensore di Ricci ossia r1i2i
che ho chiamato rinva.Supponiamo di aver già calcolato
il tensore di Ricci totalmente covariante ossia r22.
Allora, supponendo che g11 sia il tensore metrico 
controvariante:

r12: scala(g11,2,r22,1)$
rinva : traccia(scala(g11,2,r22,1),1,2)$

Derivando covariantemente rinva, usando la metrica 
appropriata, si ottiene un vettore ossia:

rinva8: ratsimp(diffcov(rinva,metrica,"rinva8"))$

Notare che si sarebbe ottenuto lo stesso vettore 
usando la derivata ordinaria ossia usando la
fa_diff(tn,listavariabili,comme) nel seguente modo:

rinva0 : fa_diff(rinva,listavariabili,"rinva0")

infatti in caso di scalari la derivata covariante 
coincide con la derivata ordinaria.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
La derivata covariante del tensore di Ricci in forma 
mista è data da r128 ossia

r128:ratsimp(diffcov(scala(g11,2,r22,1),
             metrica,"r128"))$


   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
L'identità che dobbiamo verificare è la seguente:

2*r128[a][b,a]-rinva8[b] = 0

Con questa libreria si fa così:

ratsimp(2*rest(traccia(r128,1,3),-1)-rest(rinva8,-1));

Notare che con la rest(list,-1) si taglia via dalla
lista il suo ultimo elemento ossia quello che contiene
le info sul tipo di indici del tensore.
Quello che resta è dunque qualcosa di trattabile
algebricamente.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Se quindi viene un vettore di zeri è un trionfo...

E, fatte le verifiche, viene proprio un vettore 
di zeri e dunque, salvo bachi nascosti, tutto,
fino qui... FUNZIONA!
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Divergenza di un tensore
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Crea la funzione divergenza di un tensore rispetto ad un 
dato indice, nota la metrica e dunque la dimensione,
l'ordine del tensore e il tensore metrico controvariante.
Fa uso della derivata covariante per cui questa funzione 
è applicabile a qualunque tensore di ordine positivo.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
divergenza(tn,ki,metrica):=block([k,tm,j],
    if not(tensorp(tn)) then
       return("Errore: non tensore"),
    if not(numberp(ki)) then
       return("Errore: non numero intero"),
    k:floor(ki),
    if 1>k then
       return("Errore: indice non positivo"),
    j:indici(tn),
    if k>length(j) then
       return("Errore: indice troppo elevato"),
    if mod(j[k],2)=1 then tm:diffcov(tn,metrica,"(div)")
    else ( tm:scala(metrica[3],1,tn,k),
           tm:diffcov(tm,metrica,"(div)")),
    traccia(tm,k,ordine(tm))
    )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
libmia:cons("divergenza(tn,ki,metrica)",libmia)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
In complesso le funzioni attualmente definite in questa
libreria sono queste:
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
libmia:sort(libmia);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Conclusione...provvisoria
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Chiunque ami le formule matematiche ricche di simmetrie non può
non ammirare la bellezza della matrice del tensore di Ricci
quando la metrica è quella di Reissner Nordstrom ma non nella
forma originale bensì nella forma proposta da Schild.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
La matrice del tensore metrico, nel caso della variante di
Schild, non è più diagonale bensì piena e questo ha un 
effetto piuttosto... devastante sulla velocità del calcolo
effettuato da Maxima.... ma il risultato è il poter
fare calcoli in coordinate "quasi cartesiane" visto che
al tendere all'infinito della distanza dal buco nero la metrica
ritorna ad essere quella pseudoeuclidea della Relatività
Speciale
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Anche se i normali buchi neri macroscopici ( quello al centro
della via Lattea, ad esempio ) sono neutri ritengo molto
opportuno ipotizzare che il buco nero gigantesco abbia la
carica di almeno un elettrone... per i seguenti motivi.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
In Relatività lo stesso problema fisico puo' essere trattato
con infinite metriche diverse... come nella meccanica classica.
Infatti, anche in meccanica classica si possono usare 
i piu' disparati sistemi di riferimento che producono,
derivando lo Jacobiano, diversissimi tensori metrici.
Cito a memoria... sistema cilindrico, sferico, parabolico,
toroidale, ellissoidale etc... Di tutti questi sistemi
e di molti altri, con il package ctensor di
Maxima è possibile calcolare l'appropriato 
tensore metrico.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Quello pero' che NON CAMBIA è il valore del prodotto
scalare ossia il modulo delle forze in azione.
In Relatività, giustamente, più che di "scalare" si parla
di "invariante" perché appunto, una quantità 
invariante non cambia qualunque sia il sistema di
riferimento e dunque la sua metrica associata. 
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Ad una data distanza dal buco nero esiste una forza 
necessaria per impedire che l'oggetto cada nel
buco nero. Se anche il buco nero ha la carica di un
solo elettrone, dato che la forza elettrostatica 
repulsiva dipende dal prodotto delle due cariche
( dello stesso segno ) dei due oggetti che si
respingono. Dunque esisterà sempre una carica della
particella tale da impedire che la particella
acceleri verso il buco nero ossia cada in esso.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Ma l'entità di questa carica è UN INVARIANTE !
ed è un invariante il modulo della forza
necessaria per impedire alla particella ferma
di cominciare ad accelere
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Numericamente il problema ossia il calcolo 
di questa forza è facile perchè si tratta di
risolvere un problemino di ELETTROSTATICA,
sia pure in uno spazio deformato dalla presenza
del grande ( o piccolo ) buco nero.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Siccome la forza è un vettore il suo modulo
è un invariante e dunque non cambia qualunque
sia il sistema di riferimento ( da cui deriva la
metrica ) che sto usando.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Insomma ..non importa che stia usando, in pratica
la metrica di Schwarzschild o quella di Reissner
Nordstrom o quella di chiunque altro inventi 
una sua metrica... Si pensi a quella del 
reverendo Lemaitre:
http://en.wikipedia.org/wiki/Lemaitre_metric

IL MODULO DELLA FORZA CHE EQUILIBRA L'ACCELERAZIONE
DI GRAVITA' NON CAMBIA ossia LA RELATIVITA' NON E' 
ASSOLUTAMENTE RELATIVA ma INVARIANTE ed ASSOLUTA !
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Chiusura file di salvataggio.
   [wxMaxima: subsect end   ] */

/* [wxMaxima: input   start ] */
print("Fine della libreria tensoriale. Vedere libmia")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
closefile();
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Non resta che utilizzare questa libreria effettuando
un load(salvoqui).
   [wxMaxima: comment end   ] */

/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$