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

/* [wxMaxima: title   start ]
Esempio n.2 d'uso della libreria
libtensori.mc
   [wxMaxima: title   end   ] */

/* [wxMaxima: comment start ]
Divergenze ed altri operatori differenziali...
Vedere:
http://www.elegio.it/mc2/maxwell-generale.html
http://www.elegio.it/mc2/Ricci-Riemann.html
   [wxMaxima: comment end   ] */

/* [wxMaxima: subsect start ]
Premessa: carico la libreria
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Utilizza la libreria tensoriale. Deve dunque caricarla...
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
load_questo:"libtensori.mc";
/* [wxMaxima: input   end   ] */

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

/* [wxMaxima: comment start ]
Dato che nel file che voglio rileggere con una load(...) ci
può essere una closefile(), per non avere segnalazioni di
errore apro un file di dribbling il cui contenuto
non mi interessa e dunque non importa se lo sovrascriverò.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
if atom(path_iniziale) then writefile("eliminando.mc");
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Ora amplia il path
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
altracartella:"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: comment start ]
Ora amplio il path iniziale aggiungendogli la 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(altracartella,path_iniziale);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Dopo avere ampliato il path, cerco di caricare la libreria tensoriale.
   [wxMaxima: comment end   ] */

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

/* [wxMaxima: comment start ]
Ora salva su un nuovo file ossia salvoqui...
   [wxMaxima: comment end   ] */

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

/* [wxMaxima: comment start ]
Guardo i nomi delle funzioni della libreria.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
libmia;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Creiamo una metrica per la sperimentazione
delle formule:
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
Usiamo la più generica metrica diagonale di uno spazio
a quattro dimensioni.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
ct_coords: [t,r,h,p];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
depends([gt,gr,gh,gp],[t,r,h,p]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
ggg22:matrix(
[gt,0,0,0],
[0,gr,0,0],
[0,0,gh,0],
[0,0,0,gp]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
metrica:fa_metrica(ggg22,ct_coords)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Divergenza di un vettore.
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
La divergenza diver di un vettore ( che è uno scalare ) 
è data, in generale, da (diver:A18[1][i,i]) ma può essere 
calcolata senza bisogno di fare derivate covarianti 
ossia detto sqrtabsdet la radice quadrata del valore 
assoluto del determinante del tensore metrico covariante, 
si ottiene con la formula:
(AD1[j]:sqrtabsdet*A1[i], diver:AD10[1][i,i]/sqrtabsdet )
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Scrivo le derivate parziali in forma abbreviata ossia
indicando con un pedice la variabile di derivazione.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
derivabbrev:true;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
metrica[2];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
La radice quadrata del valore assoluto del determinante
vale:
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
rest(metrica[2],-1);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
metrica[3];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
metrica[4];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
metrica[5];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
metrica[6];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
metrica[7];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Data l'importanza dei simboli di Christoffel di seconda specie
ripeto qui la stampa delle quattro matrici che costituiscono
appunto lo pseudotensore del terzo ordine.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
ch122:[0,0,0,0]$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
ch122[1]:matrix([gt[t]/(2*gt),gt[r]/(2*gt),gt[h]/(2*gt),gt[p]/(2*gt)],
[gt[r]/(2*gt),-gr[t]/(2*gt),0,0],[gt[h]/(2*gt),0,-gh[t]/(2*gt),0],
[gt[p]/(2*gt),0,0,-gp[t]/(2*gt)]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
ch122[2]:matrix([-gt[r]/(2*gr),gr[t]/(2*gr),0,0],
[gr[t]/(2*gr),gr[r]/(2*gr),gr[h]/(2*gr),gr[p]/(2*gr)],
[0,gr[h]/(2*gr),-gh[r]/(2*gr),0],
[0,gr[p]/(2*gr),0,-gp[r]/(2*gr)]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
ch122[3]:matrix([-gt[h]/(2*gh),0,gh[t]/(2*gh),0],
[0,-gr[h]/(2*gh),gh[r]/(2*gh),0],
[gh[t]/(2*gh),gh[r]/(2*gh),gh[h]/(2*gh),gh[p]/(2*gh)],
[0,0,gh[p]/(2*gh),-gp[h]/(2*gh)]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
ch122[4]:matrix([-gt[p]/(2*gp),0,0,gp[t]/(2*gp)],
[0,-gr[p]/(2*gp),0,gp[r]/(2*gp)],
[0,0,-gh[p]/(2*gp),gp[h]/(2*gp)],
[gp[t]/(2*gp),gp[r]/(2*gp),gp[h]/(2*gp),gp[p]/(2*gp)]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
La radice quadrata del valore assoluto del determinante del
tensore metrico covariante svolge un ruolo fondamentale
in parecchie formule. Ovviamete, per un tensore metrico
diagonale specificato simbolicamente vale:
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
metrica[11][1];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Crea la funzione divergenza del tensore rispetto ad un dato indice.
Fa uso della derivata covariante per cui questa funzione è
applicabile a qualunque tensore di ordine positivo.

Nota: qui serve per test ma è già inclusa nella libreria...
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
xdivergenza(tn,kkk,metrica):=block([k,tm,j],
    if not(tensorp(tn)) then
       return("Errore: non tensore"),
    if not(numberp(kkk)) then
       return("Errore: non numero intero"),
    k:floor(kkk),
    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,"!")
    else ( tm:scala(metrica[3],1,tn,k),
           tm:diffcov(tm,metrica,"!")),
    traccia(tm,k,ordine(tm))
    )$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Scrivo un vettore controvariante generico:
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
depends([At,Ar,Ah,Ap],ct_coords);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
A1:[At,Ar,Ah,Ap,[1,"Vettore controvariante",4]];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Sfrutto le proprietà della divergenza di un vettore controvariante
per non calcolare la derivata covariante.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
AD1:append(metrica[11][1]*rest(A1,-1),
 [[1,"sqrt(|g22|)*A1",4]]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
vala:traccia(append([radcan(fa_diff(AD1,ct_coords,"gradiente di AD1")[1]/metrica[11][1])],
 [[1,0,"derivata ordinaria",4]]),1,2);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
xssa:xdivergenza(A1,1,metrica)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Controllo se ci sono differenze rispetto alla versione
in libreria...
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
tmeno(xssa,divergenza(A1,1,metrica));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
I due modi di calcolare la divergenza di un vettore controvariante
coincidono. Infatti...
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
ratsimp(rest(ratsimp(xssa),-1)-rest(vala,-1));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Mi aspetto che siano coincidenti anche i due modi di calcolare la 
divergenza di un vettore covariante.
Sia B2 un vettore covariante...
Va prima convertito in controvariante e poi derivato.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
depends([Bt,Br,Bh,Bp],ct_coords);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
B2:[Bt,Br,Bh,Bp,[2,"Vettore covariante",4]];
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Calcolo la sua divergenza in modo standard...
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
xssb:xdivergenza(B2,1,metrica)$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Controllo se ci sono differenze rispetto alla versione in libreria
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
tmeno(xssb,divergenza(B2,1,metrica));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Ora invece sfrutto le proprietà della divergenza di un vettore...
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
B1:scala(metrica[3],2,B2,1);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
BD1:append(metrica[11][1]*rest(B1,-1),
 [[1,"sqrt(|g22|)*B1",4]]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
valb:append([radcan(fa_diff(BD1,ct_coords,
 "gradiente di BD1")[1]/metrica[11][1])],
 [[1,0,"derivata ordinaria",4]]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
scab:traccia(valb,1,2);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
I due metodi, quello standard e quello specifico del caso
debbono dare lo stesso risultato. Infatti...
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
ratsimp(xssb[1]-scab[1]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Ora lavoro su tensori del secondo ordine...
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
depends([a11,a12,a13,a14,a21,a22,a23,a24,
a31,a32,a33,a34,a41,a42,a43,a44],ct_coords);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
mat:matrix(
  [a11,a12,a13,a14],
  [a21,a22,a23,a24],
  [a31,a32,a33,a34],
  [a41,a42,a43,a44]);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
tsmat:tmat11(mat+transpose(mat),metrica);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
tamat:tmat11(mat-transpose(mat),metrica);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Ridefinisco qui varie funzioni già presenti
nella libreria per poterle modificare facendo
delle stampe intermedie.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
xsomma(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(rest(tena,-1)+rest(tenb,-1),
      [append(ina,["(Somma tensori)",na])])
   )$
/* [wxMaxima: input   end   ] */

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

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

/* [wxMaxima: input   start ] */
xtcov(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 ] */
xtcontrov(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 ] */
xtmat11(mat,metrica):=block([],
    if not(matrixp(mat)) then
       return("Non matrice"),
    [mat,[1,1,"(tm11)",dimensione(metrica[2])]])$
/* [wxMaxima: input   end   ] */

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

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

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

/* [wxMaxima: input   start ] */
xtcov([1,2,3],metrica);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
xtcontrov([1,2,3],metrica);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
tab:tsomma(tamat,tsmat);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
tprod(tsca(1/2,metrica),tab);
/* [wxMaxima: input   end   ] */

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