/* FileName: Riemann.mac Contents: maxima program to calculate the curvature tensor of a metric CreatedBy: Kodama Hideo CreationDate: 2022/11/14 < Riemann.mpl 2022/11/21: ver. 1.0 (no help vrsion) 2022/11/25: proc: DV & proc: NPC are modified. 2022/11/27: proc:NPC are modified. Checked by the Kerr metric 2022/11/27: ver. 1.1 (still without procinfo ) 2022/11/28: help proc: procinfo is added. 2022/11/28: ver. 1.2-pub : UtilityFunctions.mac is included. 2023/1/20: ver. 1.21-pub: printRiem/.../printET in procinfo() are corrected. 2023/3/31: ver. 1.22-pub: a bug in Div() and typo in procinfo() are fixed 2023/5/4: utility functions are moved to separate files LastUpdate: 2023/5/4 */ /******>TOP */ stringdisp : true$ simp: true$ /*****>packages */ scalarmatrixp: false$ stringdisp: true$ load("eigen")$ /*****>Path */ /*load("operatingsystem")$*/ loadpathname: load_pathname$ MyLibPath: pathname_directory(loadpathname)$ /*****>loading my lib programs */ /*****>list & lattice functions */ load( concat(MyLibPath,"mylib_list_v1_00_pub.mac"))$ /*****>matrix functions */ load( concat(MyLibPath,"mylib_matrix_v1_10_pub.mac"))$ /*****>dictionary variables and functions */ load( concat(MyLibPath,"mylib_dict_v1_00_pub.mac"))$ /******>Basic procedures/functions */ /******>Global variables */ _u:'_u, _d:'_d, /* index position values*/ /* FB:table(antisymmetric), */ /******>conventions */ /* T :: tensordict = [[ ["value"] = [[ [1,2, ...] = [T[1,2,..]], ... ]], [ "indextype]=[_u,_d,...], ["dim"]=[] ]] */ /******>Data to be given */ /*****>varlist */ /* Coordinates*/ /*phi:'phi, theta:'theta, psi:'psi, t:'t, r:'r, x:'x, y:'y, varlist:[t,r,theta,phi], */ /*****>metric */ /* Parameters */ /*a:'a, ds2:-(1-a/4/r)^2/(1+a/4/r)^2*dt^2+(1+a/4/r)^4*(dr^2+r^2*(dtheta^2+sin(theta)^2*dphi^2)), */ /******>Curvature calculation */ /*****>proc:df */ /* df(eq, varlist) => deq by dvarlist */ df(eq, varlist) := block ( [ dim, dvarlist ], if not listp(varlist) then ( error("Usage, df(eq, varlist::list)") ), /*end if*/ dim : length(varlist), dvarlist : map(lambda([var], concat(d, var)), varlist), return ( sum(diff(eq, varlist[i])*dvarlist[i], i, 1, dim) ) )$ /*end proc*/ /*****>proc:Riemann */ /* Procedure to calculate the connection and the Riemann tensor. The output is used to calculate derived quantities such as the Ricci tensor, covariant derivatives and so on. */ /* Riemann(varlist,ds2[,sw]) => Gdata = geometrical data sw=0 -> simplification, no display sw=1 -> simplification with stage messages sw=2 -> + some results Gdata = dict [["varlist"]=varlist, ["dvarlist"] =dvarlist, ["ds2"]=ds2, ["metric"]=[G::matrix], ["Ch"]=Chdict, ["Riem13"]=Riem13, ["Riem04"]=Riem04, ["Riem22"]=Riem22, ["Ric02"]=[Ric02], ["Ric11"] = [Ric11], ["RS"}=[RS], ["ET02"]=[ET02], ["ET11"]=[ET11], ["Weyl04"]=Weyl04 ] ds2= non-degenerate quadratic equation in dvarlist, Ch = dict[ [i,j,k]=[Ch[i; j,k]], .. ] Riem13, Riem04, Riem22, Weyl04: tensor dict of items of the type [i,j,k,l]=[T[i,j,k,l] Ric02, Ric11, ET02, ET11 : matri]x */ Riemann( varlist, ds2, [opt]) := block( [/*scalar */ sw, dim, RS, /*list */ dvarlist, xx,dxx, /* matrix */ G, GI, Ric02, Ric11, ET02, ET11, /* dict */ Ch, Riem13, Riem04,Riem22, Weyl04, /* misc */ Gdata, val ], local(xx, dxx), /****> read the input data */ sw: 0, if length(opt)>0 then ( sw: opt[1] ), /*end if*/ dim: length(varlist), /* dvarlist */ dvarlist : map(lambda([var], concat(d, var)), varlist), /****> general coordinates */ xx : varlist, dxx : dvarlist, /****>G[i,j] & GI[i,j] */ G: zeromatrix(dim,dim), ds2: expand(ds2), for ii : 1 thru dim do ( G[ii,ii]:factor(coeff(ds2,dxx[ii]^2)) ), /*end do*/ for ii : 1 thru dim do ( for jj : ii+1 thru dim do ( G[ii,jj]: 1/2*factor(coeff(coeff(ds2,dxx[ii]),dxx[jj])), G[jj,ii] : G[ii, jj] ) /*end do*/ ). /*end do*/ GI : G^^(-1), /****> Christoffel symbols */ Ch : newdict(), if sw>0 then ( print("Calculating the Christoffel symbols") ), /*end if*/ for ii : 1 step 1 thru dim do ( for jj1 : 1 step 1 thru dim do ( for jj2 : 1 step 1 thru dim do ( val : sum(1/2*GI[ii,jj3]*(diff(G[jj1,jj3],xx[jj2]) +diff(G[jj2,jj3],xx[jj1]) -diff(G[jj1,jj2],xx[jj3])), jj3, 1, dim), val : factor(trigsimp(xthru(val))), val : factor(trigsimp(xthru(val))), adddict(Ch, [ii, jj1, jj2]=[val] ) ) /*end do*/ ) /*end do*/ ), /*end do*/ /****>Riemann curvature */ if sw>0 then ( print("Calculating Curvatures") ), /*end if*/ /***>Riem13 & Riem04& Riem22 */ Riem13: newdict(), Riem04: newdict(), Riem22: newdict(), if sw>0 then ( print("Riem13[*,*,*,*], Riem04[*,*,*,*], Riem22[*,*,*,*]") ), /*end if*/ for kk : 1 thru dim do ( for ll : kk thru dim do ( for ii : 1 thru dim do ( for jj : 1 thru dim do ( val : diff(getdict(Ch, [ii,jj,ll])[1], xx[kk])-diff(getdict(Ch, [ii,jj,kk])[1], xx[ll]) +sum(getdict(Ch,[ii,kk,mm])[1]*getdict(Ch, [mm,jj,ll])[1] -getdict(Ch, [ii,ll,mm])[1]*getdict(Ch, [mm,jj,kk])[1],mm, 1, dim), val : factor(trigsimp(xthru(val))), adddict(Riem13, [ii,jj,kk,ll]=[val] ), repldict(Riem13, [ii,jj, ll, kk]= [-val] ), if sw=2 and val#0 then ( disp(printf(false, "Riem13[~d,~d,~d,~d]=~a~%", ii, jj, kk, ll, val) ) ) /*end if*/ ) /*end do*/ ), /*end do*/ for ii : 1 thru dim do ( for jj : 1 thru dim do ( val : sum(G[ii,mm]*getdict(Riem13, [mm,jj,kk,ll])[1], mm, 1, dim), val : factor(trigsimp(xthru(val))), adddict(Riem04, [ii,jj,kk,ll]=[val] ), repldict(Riem04, [ii, jj, ll, kk]=[-val]), if sw=2 and val#0 then ( disp(printf(false, "Riem04[~d,~d,~d,~d]=~a~%", ii, jj, kk, ll, val) ) ), /*end if*/ val: sum(GI[jj,mm]*getdict(Riem13,[ii,mm,kk,ll])[1], mm, 1, dim), val : factor(trigsimp(xthru(val))), adddict(Riem22, [ii, jj, kk, ll]=[val]), repldict(Riem22, [ii, jj, ll, kk]=[-val]), if sw=2 and val#0 then ( disp(printf(false,"Riem22[~d,~d,~d,~d]=~a~%", ii, jj, kk, ll, val )) ) /*end if*/ ) /*end do*/ ) /*end do*/ ) /*end do*/ ), /*end do*/ /****>Ricci curvatures */ if sw>0 then ( print("Ricci curvatures, Ric02[*,*], Ric11[*,*], RS") ), /*end if*/ /***>Ric02 */ Ric02 : zeromatrix(dim,dim), for ii : 1 step 1 thru dim do ( for jj : ii step 1 thru dim do ( val: sum(getdict(Riem13,[kk,ii,kk,jj])[1],kk,1, dim), val : factor(trigsimp(xthru(val))), Ric02[ii,jj] : val, if ii#jj then ( Ric02[jj,ii]: val ) /*end if*/ ) /*end do*/ ), /*end do*/ /***>Ric11 */ Ric11: zeromatrix(dim,dim), for ii : 1 step 1 thru dim do ( for jj : 1 step 1 thru dim do ( Ric11[ii,jj]: factor(trigsimp(xthru(sum(GI[ii,kk]*Ric02[kk,jj], kk, 1, dim)))) ) /*end do*/ ), /*end do*/ /***>Rs */ RS: factor(trigsimp(xthru(sum(Ric11[ii,ii], ii, 1, dim)))), if sw>0 then ( print("Einstein tensor E02, E11") ), /*end if*/ /****>Einstein tensor */ if sw=2 then ( print("Einstein tensors, ET02[*,*], ET11[*,*]") ), /*end if*/ ET02: zeromatrix(dim,dim), /* Einstein tensor*/ ET11: zeromatrix(dim, dim), /* ET11[i,j]=ET^i{}_j*/ for ii : 1 thru dim do ( for jj : ii thru dim do ( ET02[ii,jj] : factor(trigsimp(xthru( Ric02[ii,jj]-1/2*RS*G[ii,jj]))), if jj=ii then ( ET11[ii,jj]:factor(trigsimp(xthru( Ric11[ii,jj]-1/2*RS))) ) else ( ET02[jj,ii]: ET02[ii,jj], ET11[ii,jj]: Ric11[ii,jj], ET11[jj,ii]: ET11[ii,jj] ) /*end if*/ ) /*end do*/ ), /*end do*/ /****>Weyl tensor */ if sw>0 then ( print("Weyl tensor, Weyl04[*,*,*,*]") ), /*end if*/ Weyl04 : newdict(), if dim>=4 then ( for ii : 1 thru dim do ( for jj : ii thru dim do ( for kk : 1 thru dim do ( for ll : kk thru dim do ( val: getdict(Riem04,[ii,jj,kk,ll])[1] -1/(dim-2)*(G[ii,kk]*Ric02[jj,ll]-G[ii,ll]*Ric02[jj,kk] +G[jj,ll]*Ric02[ii,kk]-G[jj,kk]*Ric02[ii,ll]) +1/(dim-1)/(dim-2)*(G[ii,kk]*G[jj,ll]-G[ii,ll]*G[jj,kk])*RS, val : factor(trigsimp(xthru(val))), adddict(Weyl04, [ii,jj,kk,ll]=[val]) , repldict(Weyl04, [ii,jj,ll,kk]=[-val]), repldict(Weyl04, [jj,ii,kk,ll]=[-val]), repldict(Weyl04, [jj,ii,ll,kk]=[val]) ) /*end do*/ ) /*end do*/ ) /*end do*/ ) /*end do*/ ), /*end if*/ /****>constructing Gdata dict */ Gdata: newdict( ["varlist"]=varlist, ["dvarlist"]= dvarlist, ["dim"]=[dim] ), adddict(Gdata, ["ds2"]= [ds2] ), adddict(Gdata, ["metric"]=[G]) , adddict(Gdata, ["Ch"]=Ch ), adddict(Gdata, ["Riem13"]=Riem13 ), adddict(Gdata, ["Riem04"]=Riem04), adddict(Gdata, ["Riem22"]=Riem22 ), if dim>2 then ( adddict(Gdata, ["Weyl04"]=Weyl04 ) ), /*end if*/ adddict(Gdata, ["Ric02"]=[Ric02] ), adddict(Gdata, ["Ric11"]=[Ric11] ), adddict(Gdata, ["RS"]=[RS] ), adddict(Gdata, ["ET02"]=[ET02] ), adddict(Gdata, ["ET11"]=[ET11] ), /****>output */ if sw>0 then ( print("Done.") ), /*end if*/ return (copy( Gdata)) )$ /*end proc: Rieman*/ /******>print procedures for Riemann */ /*****>proc:printCh */ /* printCh(Gdata) => printing sum_{j1,j2} Ch[i,j1,j2] dx[j1] dx[j2] */ printCh(Gdata) := block( [/*scalar*/ dim, /*list*/ varlist, dx, /*dict*/ Ch, /*misc*/ tmp ], Ch: getdict(Gdata, ["Ch"]), dim: getdict(Gdata, ["dim"])[1], varlist: getdict(Gdata, ["varlist"]), dx : getdict(Gdata, ["dvarlist"]), for ii : 1 thru dim do ( tmp: sum(sum(getdict(Ch, [ii,jj1,jj2])[1] *dx[jj1] *dx[jj2], jj1, 1, dim),jj2, 1, dim), tmp: apply(collectterms, cons(tmp, makelist( dx[jj1], jj1, 1, dim))) , disp( printf(false, "Gamma^~a_(jk) dx^j dx^k", ii) = tmp) ) /*end do*/ )$ /*end proc: printCh*/ /*****>proc:printRiem */ /* printRiem(Gdata[,Ttype]) => printing Riem04/13/22[*,*,*,*] <>0 Ttype = "Riem04", "Riem13", "Rime22", "Weyl04" */ printRiem(Gdata, [opt]) :=block( [/*scalar */ dim, outcount, /* string */ Ttype, /* dict */ RiemT, /* misc*/ tmp ], dim: getdict(Gdata, ["dim"])[1], Ttype:"Riem13", if length(opt)>0 then ( Ttype: opt[1], if not element(Ttype, ["Riem04", "Riem13", "Rime22", "Weyl04"]) then ( error("Tensor type value should be \"Riem04\" or \"Riem13\" or \"Riem22\" or \"Weyl04\". ") ) /*end if*/ ), /*end if*/ RiemT: getdict( Gdata, [Ttype]), outcount:0, for ii1 : 1 thru dim do ( for ii2 : 1 thru dim do ( for jj1 : 1 thru dim do ( for jj2 : jj1+1 thru dim do ( if (tmp: getdict(RiemT, [ii1,ii2,jj1,jj2]))#[] then ( tmp: ev(tmp[1], simp), if tmp#0 then ( disp(printf(false, "~a[~d,~d,~d,~d]", Ttype, ii1, ii2, jj1, jj2)= tmp), outcount : outcount +1 ) /*end if*/ ) /*end if*/ ) /*end do*/ ) /*end do*/ ) /*end do*/ ), /*end do*/ if outcount=0 then ( disp("All components vanish.") ) /*end if*/ )$ /*end proc: printRiem*/ /*****>proc:printRic */ /* printRic(Gdata [, Ttype]) => print Ric02/11/ET02/ET11[*,*]<>0 Ttyep = "Ric02", "Ric11", "ET02", "ET11" */ printRic(Gdata, [opt]) := block( [ /*scalar */ dim,outcount, /*string*/ Ttype, /*dict*/ RicT, /* misc */ tmp ], dim: getdict(Gdata, ["dim"])[1], Ttype:"Ric11", if length(opt)>0 then ( Ttype: opt[1], if not element(Ttype, ["Ric02", "Ric11", "ET02", "ET11"]) then ( error("Tensor type value should be \"Ric02\" or \"Ric11\" or \"ET02\" or \"ET11\". ") ) /*end if*/ ), /*end if*/ RicT: getdict( Gdata, [Ttype])[1], outcount:0, for ii1 : 1 thru dim do ( for ii2 : 1 thru dim do ( tmp: ev(RicT[ii1,ii2],simp), if tmp # 0 then ( disp(printf(false,"~a[~d,~d]",Ttype, ii1,ii2)= tmp), outcount:outcount+1 ) /*end if*/ ) /*end if*/ ), /*end do*/ if outcount=0 then ( print("All components vanish.") ) /*end if*/ )$ /*end proc: printRic*/ /*****>proc: tensordictp */ /* tensordictp(T[,sw]) => true/false sw=0 (default): no message sw=1 : with message */ tensordictp(T,[opt]) := block( [/*scalar*/ sw, /*misc*/ var, dim, indextype, val, indexrange, result ], sw : 0, if length(opt)>0 then ( sw : opt[1] ), /*end if*/ var : getdict(T, ["dim"]), indextype: getdict(T, ["indextype"]), val : getdict(T,["value"]), if var=[] or not listp(var) or not nonnegintegerp(dim: var[1]) then ( if sw>0 then ( print(" no or improper item on [\"dim\"].") ), /*end if*/ return ( false) ) elseif ( indextype=[] or not listp(indextype) or not subsetp(setify(indextype),{_u, _d}) ) then ( if sw>0 then ( print(" no or improper item on [\"indextype\"].") ), /*end if*/ return (false) ) elseif not dictp(val) then ( if sw>0 then ( print("the data for [\"value\"] is not dict type.") ), /*end if*/ return (false) ) else ( result: true, indset : setify(makelist(i,i, 1, dim)), for x in dicttaglist(val) do ( if not subsetp(setify(x), indset) then ( if sw>0 then ( print("an index key in the data for [\"value\"] is wrong:", x) ), /*end if*/ result : false, break ) /*end if*/ ), /*end do*/ return (result) ) /*end if*/ )$ /*end proc: tensordictp*/ /*****>proc:printTensor */ /* printTensor(T) => printing the nozero components of T T= tensordict. */ printTensor(T) := block( [/*scalar*/ dim,Trank, nonzerocount, /* list */ indextype, indset, index, /* dict */ valdict /*misc*/ ], if not tensordictp(T) then ( error("The argument T must be a tensor-type dict") ), /*end if*/ dim: getdict(T,["dim"])[1], indextype: getdict(T,["indextype"]), Trank: length(indextype), valdict: getdict(T, ["value"]), if Trank=0 then ( disp("Tensor type is scalar."), disp("value"= getdict(valdict, [])[1] ) ) else ( print("Tensor type is tensor of the indextype"=indextype), nonzerocount:0, indset: mklattice(Trank,[1,dim]), for index in indset do ( if (val: getdict(valdict, index))#[] and val[1]#0 then ( if nonzerocount=0 then ( disp("Non-zero components are as follows,") ), /*end if*/ disp(index = val[1] ), nonzerocount: nonzerocount+1 ) /*end if*/ ), /*end do*/ if nonzerocount=0 then ( disp("All components vanish!") ) /*end if*/ ) /*end if*/ )$ /*end proc*/ /******>Kretschmann invariants */ /*****>proc:Kretschmann */ /* Kretschmann(Gdata) -> the Kretschmann invariant */ Kretschmann(Gdata) := block( [ dim, val, /* dict*/ Riem22] , dim: getdict(Gdata, ["dim"])[1], Riem22: getdict(Gdata, ["Riem22"]), val:0, for i : 1 thru dim-1 do ( for j : i+1 thru dim do ( for k : 1 thru dim-1 do ( for l : k+1 thru dim do ( val: val+4*getdict(Riem22,[i,j,k,l])[1]*getdict(Riem22,[k,l,i,j])[1] ) /*end do*/ ) /*end do*/ ) /*end do*/ ), /*end do*/ return (val) )$ /*end proc: Kretschmann*/ /*****>proc:KretschmannS */ /* KretschmannS(Gdata) => Kretschmann1, Kretschmann2 for static metric with xx[1] in the static direction */ KretschmannS(Gdata):= block( [ dim, Kretschmann1, Kretschmann2, /*dict*/ Riem22], dim: getdict(Gdata, ["dim"])[1], Riem22: getdict(Gdata, ["Riem22"]), Kretschmann1:0, for ii : 2 thru dim do ( Kretschmann1: Kretschmann1+4*getdict(Riem22,[1,ii,1,ii])[1]^2, for jj : ii+1 thru dim do ( Kretschmann1: Kretschmann1+8*getdict(Riem22,[1,ii,1,jj])[1]*getdict(Riem22, [1,jj,1,ii])[1] ) /*end do*/ ), /*end do*/ Kretschmann2:Kretschmann(Gdata)-Kretschmann1, return ([Kretschmann1, Kretschmann2]) )$ /*end proc: KretschmannS*/ /******>Tensor algebra */ /*****>proc:islinear */ /* islinear (eq, basis) => true/false (eq=dform?) */ islinear(eq, basis) := block( [ a ], eq: expand(eq), a: map(lambda([x], coeff(eq,x)), basis), if eq-inprod(a,basis)=0 then ( return (true) ) else ( return (false) ) )$ /*end proc*/ /*****>proc:ContractT */ /* ContractT(T::tensordict, a, b, G]) => CT::tensordict a, b: index positions, */ ContractT(T, a, b, G, [opt]) := block( [/*scalar*/ dim,Trank, s1,s2, Q, /*list*/ indextype, indset, index, newind, /*matrix*/ GI, M0, /*dict*/ CTval, Tval, CT, /*misc*/ tmp ], dim: getdict(T, ["dim"])[1], indextype: getdict(T, ["indextype"]), Trank: length(indextype), if Trank<2 then ( error("The rank of the tensor must be greater than 1") ), /*end if*/ [s1,s2] : sort([a,b]), if not ( 1<=s1 and s1proc:IP */ /* Inner product of two tensors/vectors/forms IP(V0,W0, G) => val,scalar V0,W0 :: tensordict/vectors/covectors */ IP( V0, W0, G) := block( [ /*scalar*/ dim,Trank, Q, q, N, /* list */ xx, dxx,indextypeV, indextypeW,indset, ind1,ind2, /* matrix*/ GI, IDM, /* dict */ Vval, Wval, /* misc*/ V, W ], V: ev(copy(V0), simp), W: ev(copy(W0), simp), GI: G^^(-1), dim: matrixtype(G)[1], IDM : ident(dim), if matrixp(V) and matrixp(W) then ( /*print("Case Vector"),*/ if matrixtype(V)[2]=1 and matrixtype(W)[2]=1 then ( /* vectors */ Q: inprod(V, G . W)[1,1] ) elseif matrixtype(V)[1]=1 and matrixtype(W)[1]=1 then ( /*covectors*/ Q: inprod(V, GI . W)[1,1] ) elseif (matrixtype(V)[1]=1 and matrixtype(W[2]=1) ) or (matrixtype(V)[2]=1 and matrixtype(W[1]=1) ) then ( /*vector x covector*/ Q: inprod(V, transpose(W))[1,1] ) /*end if*/ ) elseif tensordictp(V) and tensordictp(W) then ( /* tensors */ /*print("Case Tensor"),*/ indextypeV: getdict(V, ["indextype"]), indextypeW: getdict(W,["indextype"]), Trank: length(indextypeV), if length(indextypeW)#Trank then ( error("Two tensors should have the same rank") ), /*end if*/ Q : 0, indset: mklattice(Trank,[1,dim]), N: ev(dim^Trank), Vval: newdict(), for x in getdict(V, ["value"])[1] do ( if args(x)[2][1]#0 then ( adddict(Vval, x) ) /*end if*/ ), /*end do*/ Wval: newdict(), for x in getdict(W, ["value"])[1] do ( if args(x)[2][1]#0 then ( adddict(Wval, x) ) /*end if*/ ), /*end do*/ for x in Vval[1] do ( ind1: args(x)[1], for y in Wval[1] do ( ind2: args(y)[1], q: args(x)[2][1]*args(y)[2][1], for i :1 thru Trank do ( if indextypeV[i]='_u and indextypeW[i]='_u then ( q: q*G[ind1[i], ind2[i]] ) elseif indextypeV[i]='_d and indextypeW[i]='_d then ( q: q*GI[ind1[i], ind2[i]] ) else ( q: q*IDM[ind1[i], ind2[i]] ) /*end if*/ ), /*end do*/ Q: Q+ q ) /*end do*/ ) /*end do*/ ) else ( error("Usage: IP(V, W, Gdata). V, W must be vectors/covectors/tensor dicts of the same rank.") ), /*end if*/ if Q#'Q then ( return (Q) ) else ( error("Usage: IP(V,W, Gdata). V, W= vectors/covectors/tensor dicts") ) /*end if*/ )$ /*end proc: IP */ /******>General frame & Covariant derivative */ /*****>proc: mkVdict */ /* mkVdict( V, indextype) => a tensor dict for the vector v V= list / vector /covector */ mkVdict( V, indextype) := block( [ dim, Vdict], if matrixp(V) then ( V : listme(V) ), /*end if*/ dim : length(V), Vdict : newdict( ["dim"]=[dim], ["indextype"]=indextype ), adddict( Vdict, ["value"]=[makelist( [i]= [V[i]] , i, 1, dim )] ), return (copy(Vdict)) )$ /*end proc: mkVdict */ /*****>proc: mkT2dict */ /* mkT2dict( T2, indextype) => a tensor dict for the vector v T2 = a matrix or a list of lists representing a 2nd rank tensor */ mkT2dict( T2, indextype) := block( [ dim, m,n, T2dict], if listp(T2) then ( T2 : apply(matrix, T2) ), /*end if*/ if matrixp(T2) then ( [m,n] : matrixtype(T2), if m#n then ( error("The matrix must be a square matrix") ), /*end if*/ dim : m ) else ( error("Usage: mkT2dict(T2: square matrix, indextype)") ), /*end if*/ T2dict: newdict(["dim"]=[dim], ["indextype"]=indextype), T2val: newdict(), for i :1 thru dim do ( for j:1 thru dim do ( adddict(T2val, [i,j]= [T2[i,j]] ) ) /*end do*/ ), /*end do*/ adddict(T2dict, ["value"]= T2val ), return (copy(T2dict)) )$ /*end proc: mkT2dict */ /*****>proc:Covder0 */ /* Covder0(T, i, [j,k,, ], Gdata) => val::scalar = (nabla_i T)[j, k, ..] T = a tensor dict i , j, k = 1,.., dim */ Covder0(T, derind, indlist, Gdata) := block( [ /*scalar */ dim, Trank, /* list */ xx, indextype, Tindex, /* matrix */ /* dict */ Ch, /*misc*/ tmpindex, val ], Ch: [], /* set parameters */ Trank:0, if dictp(T) then ( indextype: getdict(T, ["indextype"]), /* =[_u, _d, , ] */ if length(indextype)=0 then ( Tindex : [], Trank: 0, T : getdict( getdict(T, ["value"]),[])[1] /* converted to a scalar function */ ) else ( Tindex: copy(indlist), Trank: length(indextype) ) /*end if*/ ) else ( Tindex: [], Trank: 0 /* treated as a scalar function */ ), /*end if*/ dim: getdict(Gdata,["dim"])[1], if dictp(T) and getdict(T, ["dim"])[1] # dim then ( error("The dimensions of the tensor and space do not match") ), /*end if*/ Ch: getdict(Gdata,["Ch"]), xx: getdict(Gdata,["varlist"]), /* T=scalar case */ if Trank=0 then ( return (diff(T,xx[derind])) ) elseif length(Tindex) # Trank then ( print("Trank"=Trank,"Tindex"=Tindex), error("Type of Tensor does not match the index list.") ) else ( /* T=non-scalar tensor case */ Tval : getdict(T, ["value"]), val: diff(getdict(Tval, Tindex)[1], xx[derind]), for ii : 1 thru Trank do ( tmpindex: copy(Tindex), for jj : 1 thru dim do ( tmpindex[ii]: jj, if indextype[ii]='_d then ( val: val-getdict(Ch, [jj,derind,Tindex[ii]])[1]*getdict(Tval, tmpindex)[1] ) elseif indextype[ii]='_u then ( val: val + getdict(Ch, [Tindex[ii], derind, jj])[1]*getdict(Tval, tmpindex)[1] ) else ( error("Invalid indextype of T.") ) /*end if*/ ) /*end do*/ ), /*end do*/ return (val) ) /*end if*/ )$ /*end proc: Covder0*/ /*****>proc: CovDT */ /* CovDT(T, Gdata) -> the tensordict DT */ CovDT(T, Gdata) := block( [ /*scalar */ dim, Trank, /* list */ indextype, indextype1, index, indexset1, /* dict */ DTval, DTdict ], if not tensordictp(T,1) then ( error("Usage: DT ( T:: tensordict, Gdata).") ), /*end if*/ dim : getdict(T, ["dim"])[1], if getdict(Gdata,["dim"])[1] # dim then ( error("The dimensions of the tensor and Gdata do not match.") ), /*end if*/ indextype : getdict(T, ["indextype"]), Trank : length(indextype), indextype1 : cons('_d, indextype), indexset1 : mklattice( Trank+1, [1, dim]), DTval : newdict(), for index in indexset1 do ( adddict(DTval, index= [ Covder0(T, index[1], slist(index, 2, Trank+1), Gdata) ] ) ), /*end do*/ return ([[ ["dim"]=[dim], ["indextype"]=indextype1, ["value"]=DTval ]]) )$ /*end proc: DT */ /*****>proc: CovDnT */ /* CovDnT(T, n, Gdata) -> tensordict D^nT */ CovDnT(T, n, Gdata) := block( if not nonnegintegerp(n) then ( error("Usage: CovDnT(T:: tesor dict, n::nonnegative integer, Gdata)") ), /*end if*/ if n=0 then ( return (copy(T)) ) else ( return (copy( CovDT( CovDnT(T, n-1, Gdata), Gdata))) ) /*end if*/ )$ /*end proc: CovDnT*/ /*****>proc:Covder */ /* Covder(T, derindex::list, Tindex::list, Gdata)=> D_i1 , . D_in T[j,k,,, .] T= a tensordict */ Covder(T, derind, Tind, Gdata) := block( [/*scalar */ derrank, /* list */ Tind1 ], derrank : length(derind), if derrank=0 then ( return( getdict(getdict(T,["value"]),Tind)[1]) ) elseif derrank= 1 then ( return( Covder0(T, derind[1], Tind, Gdata)) ) else ( Tind1 : append(slist(derind, 2, derrank), Tind), return( Cover0(DnT(T, derrank-1, Gdata), derind[1], Tind1, Gdata)) ) /*end if*/ )$ /*end proc: Covder*/ /*****>proc:LaplacianT */ /* LaplacianT(T,Gdata) -> Laplacian of tensor T T ; tensor dict => (D^2 T) */ LaplacianT(T, Gdata) := block ( [tmp], tmp: CovDnT(T, 2, Gdata), return (copy(ContractT(tmp, 1, 2, getdict(Gdata,["metric"])[1]))) )$ /*end proc*/ /*****>proc:DV */ /* DV(X, Y, Gdata) => D_X Y X,Y = vectors by lists => list */ DV (X, Y, Gdata) := block( [ /*scalar*/ dim, /*list*/ Z, indextype, /*dict*/ Ydict ], Ydict:[], dim: getdict(Gdata, ["dim"])[1], indextype: ['_u], if matrixp(X) and matrixtype(X)[2]=1 or matrixtype(X)[1]=1 then ( X : listme(X) ), /*end if*/ if matrixp(Y) then ( if matrixtype(Y)[2] =1 then ( indextype : ['_u] ) elseif matrixtype(Y)[1]=1 then ( indextype: ['_d] ) else ( error("Wrong data for the vector Y") ), /*end if*/ Y: listme(Y) ), /*end if*/ if not listp(X) or not listp(Y) or length(X)#dim or length(Y)#dim then ( error("Invalid dimension of the vectors") ), /*end if*/ Ydict : mkVdict( Y, indextype), Z: makelist(0, i, 1, dim), for ii : 1 thru dim do ( Z[ii]: sum(X[jj]*Covder0(Ydict, jj, [ii], Gdata), jj, 1, dim) ), /*end do*/ return (copy(Z)) )$ /*end proc*/ /*****>proc:Div */ /* Div(T, indexlist, Gdata[, k]) => D_i T^i{}_j.., . k: sum index position */ Div(T, indexlist, Gdata, [opt]) := block( [ /*list*/ indextype, indexlist1, /*matrix */ GI, G, /*misc*/ k, dim, val, Trank ], if not tensordictp(T) then ( error("Usage: Div(T:: tensordict, k::sumindex position, indexlist, Gdata") ), /*end if*/ k:1, if length(opt)>0 then ( k: opt[1] ), /*end if*/ indextype: getdict(T,["indextype"]), Trank: length(indextype), if not integerp(k) or k<1 or k>Trank then ( error("The contraction index position is invalid.") ), /*end if*/ dim: getdict(Gdata,["dim"])[1], if getdict(T, ["dim"])[1]# dim then ( error("The dimensions of the tensor T and space do not matsch.") ), /*end if*/ G: getdict(Gdata, ["metric"])[1], GI: G^^(-1), val :0, for i :1 thru dim do ( if indextype[k]='_u then ( indexlist1: append( slist(indexlist, 1, k-1), [i], slist(indexlist, k, Trank-1)), val: val+ Covder0(T, i , indexlist1, Gdata) ) else ( indexlist1: append( slist(indexlist, 1, k-1), [i], slist(indexlist, k, Trank-1)), for j:1 thru dim do ( val: val+ GI(i,j)*Covder0(T, j , indexlist1, Gdata) ) /*end do*/ ) /*end if*/ ), /*end do*/ return (val) )$ /*end proc: Div*/ /*****>proc:DivT */ /* DivT(T::tensordict, Gdata[,k]) => DivTdict:: tensordict DivTdict [,,, [a,b,, .]=sum(D_p T[p,a,b,, ],a=1, dim) */ DivT(T, Gdata, [opt]) := block( [k, tmp], if not tensordictp(T) then ( error("Usage: Div(T:: tensordict, k::sumindex position, indexlist, Gdata") ), /*end if*/ k:1, /* sum index position */ if length(opt)>0 then ( k: opt[1] ), /*end if*/ tmp: CovDT(T, Gdata), return( copy(ContractT(tmp,1, k+1,getdict(Gdata,["metric"])[1]) )) )$ /*end proc*/ /******>Connection and curvature forms */ /*****>proc:CoordBV */ /* CoordBV(i,dim)=> [0,..., 0,1,0,...,0] (i-th component) :coordinate basis vector */ CoordBV(i, dim):= block( if not ( integerp(i) and i>0 and integerp(dim) and i<=dim ) then ( error("Usage: CoordBV(i::positive integer, dim::positive integer), i<=dim) ") ), /*end if*/ return ( append( makelist(0, j, 1, (i-1)), [1], makelist(0, k, (i+1), dim) ) ) )$ /*end proc: CoordBV*/ /*****>proc: ConnectionForm */ /* ConnectionForm(FormBasis::list, Gdata[, sw]) => OmegaM=(omega[a,b])::matrix FormBasis = [FB[1]=,, .,FB[n]=] => Spin coefficients, omega[a,b] theta[i]= sum_m theta[i,m] dxx[m] /* one forms */ sw=0 => omega[a,b]=omega[a,b][mu] dx^mu sw=1 => omega[a,b]=omega[a,b][c] FB[c] */ ConnectionForm(FormBasis, Gdata,[opt]) := block( [/*scalar*/ sw, dim, /* list */ FB1, BasisList, /* matrix*/ OmegaM,ThetaM, DEM, /*dict*/ /* array*/ /*misc*/ xx, dxx, EM, nu, dxbyFB, tmp ], /* set parameters */ FB1: map(lambda([var], args(var)[1]), FormBasis), dim: getdict(Gdata, ["dim"])[1], xx: getdict(Gdata, ["varlist"]), dxx: getdict(Gdata, ["dvarlist"]), if length(FormBasis)#dim then ( print("nops(FormBasis)"=nops(FormBasis),"dim"=dim), error("Dimension of the form basis do (es not match Gdata") ), /*end if*/ sw:0, if length(opt)>0 then ( sw: opt[1] ), /*end if*/ BasisList : map( lambda([x], args(x)[2]), FormBasis), ThetaM : coefmatrix( BasisList, dxx), /* FB[a] = sum( ThetaM[a, mu] * dxx[mu], mu, 1, dim) */ EM: transpose(ThetaM^^(-1)), /*EM[a,mu]=(e_a{}^mu)*/ tmp: makelist( sum(dxx[mu]*DV(CoordBV(mu,dim), listme(row(EM,i)), Gdata), mu, 1, dim), i, 1, dim), /* DEM[i,nu]=e_j^nu omega^j{}_i = (nabla e_i)^nu */ DEM : apply(matrix, tmp), OmegaM: ThetaM. transpose(DEM), if sw>0 then ( dxbyFB: makelist(dxx[mu]=sum(FB1[i]*EM[i,mu], i, 1, dim), mu, 1, dim), OmegaM: collectterms( ev(OmegaM, dxbyFB), FB1) ) , /*end if*/ return ( copy(OmegaM)) )$ /*end proc*/ /*****>proc:RiemForm */ /* RiemForm(FormBasis,Gdata[, sw]) => Matrix( RF[a,b]) Calculating the Curvature form FormBasis,[FB[1]=theta[1],, .,FB[n]=theta[n]] theta[i]= sum_m theta[i,m] dxx[m] sw=0 => RF[a,b]=(1/2)R^a_{b mu nu} dxdx[mu,nu] sw=1 => RF[a,b]=(1/2)R^a_{b c d} FB2[c,d] */ RiemForm(FormBasis, Gdata, [opt]) := block( [ /* list */ xx, dxx,FB2list,dx2BbyFB, dx2byCB, /* matrix */ Theta, M, EM, RF, /*name*/ dx2, CB, FB2, /*dict*/ Riem13, /*misc */ sw, dim,x ], /* set data */ dim: getdict(Gdata,["dim"])[1], xx: getdict(Gdata,["varlist"]), dxx: getdict(Gdata,["dvarlist"]), declare( dx2, antisymmetric), declare(CB, antisymmetric), dx2BbyCB : apply(append, makelist(makelist( dx2[i,j]= CB[dxx[i], dxx[j]], j, i+1, dim), i, 1, dim-1)), declare( FB2, antisymmetric), if length(FormBasis)#dim then ( error("Dimension of the form basis do (es not match Gdata") ), /*end if*/ sw:0, if length(opt)>0 then ( sw: opt[1] ), /*end if*/ /* define basis matrices */ BasisList : map( lambda([x], args(x)[2]), FormBasis), ThetaM : coefmatrix( BasisList, dxx), /* FB[a] = sum( ThetaM[a, mu] * dxx[mu], mu, 1, dim) */ EM: transpose(ThetaM^^(-1)), /*EM[a,mu]=(e_a{}^mu)*/ RF: zeromatrix(dim,dim), Riem13: getdict(Gdata, ["Riem13"]), for a : 1 thru dim do ( for b : 1 thru dim do ( RF[a,b]: sum(sum(sum(sum(getdict(Riem13,[i,j,k,l])[1]*ThetaM[a,i]*EM[b,j] *dx2[k,l],i, 1, dim),j, 1, dim),l, k+1, dim),k, 1, dim-1) ) /*end do*/ ), /*end do*/ if sw=0 then ( FB2list : map(lambda([x], args(x)[1]), dx2BbyCB), RF : ev(RF, dx2BbyCB) ) else ( dx2BbyFB:[], for i : 1 thru dim-1 do ( for j : i+1 thru dim do ( dx2BbyFB: append(dx2BbyFB, [ dx2[i,j]=sum(sum(EM[a,i]*EM[b,j]*FB2[a,b], a, 1, dim), b, 1, dim)] ) ) /*end if*/ ), /*end do*/ RF: ev(RF,dx2BbyFB), FB2list: append(makelist(makelist( FB2[i,j],j, (i+1), dim), i, 1, (dim-1))) ), /*end if*/ RF: collectterms( RF, FB2list), return (copy(RF)) )$ /*end proc: RiemForm*/ /******>Basis change */ /*****>proc:CBtoFBT */ /* CBtoFBT(T, dvarlist, FormBasis[,sw]) transforming the component rep of a tensor T : the coord basis to a general basis T = tensordict with indextype info dvarlist= the list of the coordinate 1-form basis FormBasis=[FB[1]=theta[1], , .] sw=0 => no message (default) sw=1 => with "Done" message */ CBtoFBT(T, dvarlist, FormBasis, [opt]) := block( [ /* list */ indextype, index, indset,dxx, FB1,sind,Q, q, /* matrix */ ThetaM, EM, /* dict */ ThetaMdict, EMdict,Tval, newTval,Tnew, /*misc */ sw,dim,Trank, T1 ], /* set dat */ indextype: getdict(T, ["indextype"]), if not listp(indextype) then ( error("The tensor table T should increase the indextype info") ), /*end if*/ Trank: length(indextype), if Trank=0 then ( return (T) ), /*end if*/ dxx: dvarlist, dim: length(dxx), if getdict(T,["dim"])[1]#dim then ( error("The dimensions of the tensor and space do ( not match") ), /*end if*/ FB1: map(lambda([var], args(var)[1]), FormBasis), sw:0, if length(opt)>0 then ( sw: opt[1] ), /*end if*/ /* define basis matrices */ BasisList : map( lambda([x], args(x)[2]), FormBasis), ThetaM : coefmatrix( BasisList, dxx), /* FB[a] = sum( ThetaM[a, mu] * dxx[mu], mu, 1, dim) */ EM: transpose(ThetaM^^(-1)), /*EM[a,mu]=(e_a{}^mu)*/ /* trf to the general FB */ indset: mklattice(Trank,[1,dim]), Tval : getdict(T, ["value"]), newTval: newdict(), for index in indset do ( Q: [0], for x in Tval[1] do ( [sind, q] : args(x), if q[1]#0 then ( for i :1 thru Trank do ( if indextype[i]='_u then ( q: q*ThetaM[index[i], sind[i]] ) else ( q: q*EM[ index[i], sind[i]] ) /*end if*/ ), /*end do*/ Q: Q+ q ) /*end if*/ ), /*end do*/ adddict(newTval, index=Q ) ), /*end do*/ Tnew : copy(T), repldict(Tnew, ["value"]=newTval), if sw>0 then ( print("CBtoFTB: Done") ), /*end if*/ simp:true, return (copy(Tnew)) )$ /*end proc*/ /******>D=4 specific */ /*****>NP coefficients */ /****>proc:NPC */ /* NPC(Gdata,NullFormBasis[,sw]]) => NPC NullFormBasis=[FB[1]=, ., FB[4]=, ] sw=0 => no output display, sw=1 => display the simplified Psi0,, ,Psi4 NPC: dict[ ["SC"]=SCNP, ["Weyl"]=WeylNP, ["Ric"]=RicNP, ["NullTetrad"]=NullTerad, ["dim"]=[dim], ["varlist"]=varlist, ["metric"]=[G] ] */ NPC (Gdata, NullFormBasis, [opt]) := block( [ /* list */ xx, dxx,NB, FBlist, /* matrix*/ G, GI, ThetaM, EM, Ric02, FBM, /* array */ NV, NF, DNV, Weyl04, /*dict*/ Weyl04val, SCNP, WeylNP ,RicNP, NPC, /* func*/ Weylcomp, Riccomp, /*misc*/ dim, sw, RS, val, a,b,tmp ], local( NV, NF,DNV,Weyl04, G, Weylcomp, Riccomp), /* set parameters */ dim: getdict(Gdata, ["dim"])[1], dxx: getdict(Gdata, ["dvarlist"]), G : getdict(Gdata,["metric"])[1], GI: trigsimp(xthru(G^^(-1))), if dim # 4 then ( error("NP coefficients are defined only for four dimensions") ), /*end if*/ NB: map(lambda([var], args(var)[1]), NullFormBasis), sw:0, if length(opt)>0 then ( sw: opt[1] ), /*end if*/ /*null basis vectors */ FBList : map( lambda([x], args(x)[2]), NullFormBasis), eta : blockMatrices([ -matrix([0,1],[1,0]), matrix([0,1],[1,0])]), FBM : coefmatrix( FBList, dxx), ThetaM : eta . FBM, array(NV, 4), array(NF, 4), for i:1 thru 4 do ( NF[i] : row(FBM, i), /*covariant vector*/ NV[i] : transpose(trigsimp(xthru(NF[i] . GI))) /* contravariant vector */ ), /*end do*/ /** Spin coefficients, SCNP["*"] **/ SCNP : newdict(), DNV[4,1] : covect(DV( NV[4], NV[1], Gdata)), /*vector*/ DNV[4,3]: covect(DV( NV[4], NV[3], Gdata)), val: ev( 1/2*(-(NF[2] . DNV[4,1]) + (NF[4] . DNV[4,3])), simp:false)[1,1] , /* val: factor(trigsimp(xthru(val))),*/ adddict(SCNP, ["alpha"] = [ val]), DNV[3,1]: covect(DV(NV[3], NV[1], Gdata)), DNV[3,3]: covect(DV(NV[3], NV[3], Gdata)), val : ev( 1/2*(-(NF[2] . DNV[3,1]) + (NF[4] . DNV[3,3])), simp:false)[1,1], /* val: factor(trigsimp(xthru(val))), */ adddict(SCNP, ["beta"]=[ val]), val : ev( - (NF[3] . DNV[4,1]), simp:false)[1,1], /* val: factor(trigsimp(xthru(val))), */ adddict(SCNP, ["rho"]=[ val]), val : ev ( - (NF[3] . DNV[3,1]) , simp:false)[1,1], /* val: factor(trigsimp(xthru(val))),*/ adddict(SCNP, ["sigma"]=[val]), DNV[4,2] : covect(DV( NV[4], NV[2], Gdata)), val: ev( (NF[4] . DNV[4,2]), simp:false)[1,1], /* val: factor(trigsimp(xthru(val))), */ adddict(SCNP, ["lambda"]=[ val ]), DNV[3,2] : covect(DV( NV[3], NV[2], Gdata)), val : ev( (NF[4] . DNV[3,2]), simp:false)[1,1] , /* val: factor(trigsimp(xthru(val))), */ adddict(SCNP, ["mu"]=[val]), DNV[1,1] : covect(DV( NV[1], NV[1], Gdata)), val: ev( - (NF[3] . DNV[1,1]), simp:false)[1,1], /* val: factor(trigsimp(xthru(val))), */ adddict(SCNP, ["kappa"]=[ val]), DNV[1,2] : covect(DV( NV[1], NV[2], Gdata)), val: ev( (NF[4] . DNV[1,2]), simp:false)[1,1], /* val: factor(trigsimp(xthru(val))), */ adddict(SCNP, ["pi"]=[ val ]), DNV[1,3] : covect(DV( NV[1], NV[3], Gdata)), val : ev( -1/2*((NF[2] . DNV[1,1]) + (NF[4] . DNV[1,3])) , simp:false)[1,1], /* val: factor(trigsimp(xthru(val))), */ adddict(SCNP, ["epsilon"]= [val]), DNV[2,2] : covect(DV( NV[2], NV[2], Gdata)), val: ev( (NF[4] . DNV[2,2]), simp:false)[1,1], /* val: factor(trigsimp(xthru(val))), */ adddict(SCNP, ["nu"]=[ val]), DNV[2,1] : covect(DV( NV[2], NV[1], Gdata)), val : ev( - (NF[3] . DNV[2,1]) , simp:false)[1,1], /* val: factor(trigsimp(xthru(val))), */ adddict(SCNP, ["tau"]=[val]), DNV[2,3] : covect(DV( NV[2], NV[3], Gdata)), val : ev( 1/2*(-(NF[2] . DNV[2,1]) + (NF[4] . DNV[2,3])) , simp:false)[1,1], /* val: factor(trigsimp(xthru(val))), */ adddict(SCNP, ["gamma"]=[val]), if sw>0 then ( disp(showdict(SCNP)) ), /*end if*/ print("Check pt 1"), /** Weyl, Psi[j] */ Weyl04val: getdict(Gdata, ["Weyl04"]), /* dict */ array(Weyl04, 4,4,4,4), for i1:1 thru dim do ( for i2:1 thru dim do ( for i3:1 thru dim do ( for i4:1 thru dim do ( Weyl04[i1,i2,i3,i4] : getdict(Weyl04val, [i1,i2,i3,i4])[1] ) /*end do*/ ) /*end do*/ )/*end do*/ ), /*end do*/ WeylNP: newdict(), val: sum(sum(sum(sum( Weyl04[i1,i2,i3,i4]*NV[1][i1,1]*NV[3][i2,1]*NV[1][i3,1]*NV[3][i4,1], i1,1,4), i2, 1, 4), i3, 1,4), i4, 1, 4), val: factor(trigsimp(xthru(val))), adddict(WeylNP, ["Psi[0]"]=[ val ]), val: sum(sum(sum(sum( Weyl04[i1,i2,i3,i4]*NV[1][i1,1]*NV[2][i2,1]*NV[1][i3,1]*NV[3][i4,1], i1,1,4), i2, 1, 4), i3, 1,4), i4, 1, 4), val: factor(trigsimp(xthru(val))), adddict(WeylNP, ["Psi[1]"]=[val]), val: sum(sum(sum(sum( Weyl04[i1,i2,i3,i4]*NV[1][i1,1]*NV[3][i2,1]*NV[4][i3,1]*NV[2][i4,1], i1,1,4), i2, 1, 4), i3, 1,4), i4, 1, 4), val: factor(trigsimp(xthru(val))), adddict(WeylNP, ["Psi[2]"]=[val ]), val: sum(sum(sum(sum( Weyl04[i1,i2,i3,i4]*NV[2][i1,1]*NV[1][i2,1]*NV[2][i3,1]*NV[4][i4,1], i1,1,4), i2, 1, 4), i3, 1,4), i4, 1, 4), val: factor(trigsimp(xthru(val))), adddict(WeylNP, ["Psi[3]"]=[ val ]), val: sum(sum(sum(sum( Weyl04[i1,i2,i3,i4]*NV[2][i1,1]*NV[4][i2,1]*NV[2][i3,1]*NV[4][i4,1], i1,1,4), i2, 1, 4), i3, 1,4), i4, 1, 4), val: factor(trigsimp(xthru(val))), adddict(WeylNP, ["Psi[4]"]=[ val]), if sw>0 then ( disp(showdict(WeylNP)) ), /*end if*/ print("Check pt 2"), /* Ricci, Phi[i,j] */ Ric02: getdict(Gdata, ["Ric02"])[1], /* matrix*/ RS: getdict(Gdata,["RS"])[1], RicNP : newdict( ["RS"]=[RS]), for a:1 thru 4 do ( for b: a thru 4 do ( val:0, for i1:1 thru 4 do ( for i2:1 thru 4 do ( val: ev(val+ Ric02[i1,i2]*NV[a][i1,1]*NV[b][i2,1], simp:false ) ) /*end do*/ ), /*end do*/ Riccomp[a,b] : factor(trigsimp(xthru(val))) ) /*end do*/ ), /*end do*/ adddict(RicNP, ["Phi[00]"]= [(1/2*Riccomp[1,1]) ]), adddict(RicNP, ["Phi[11]"]= [(1/2*(Riccomp[1,2]+Riccomp[3,4])) ]), adddict(RicNP, ["Phi[22]"]= [(1/2*Riccomp[2,2]) ]), adddict(RicNP, ["Phi[01]"]= [(1/2*Riccomp[1,3]) ]), adddict(RicNP, ["Phi[10]"]= [(1/2*Riccomp[1,4]) ]), adddict(RicNP, ["Phi[02]"]= [(1/2*Riccomp[3,3]) ]), adddict(RicNP, ["Phi[20]"]= [(1/2*Riccomp[4,4]) ]), adddict(RicNP, ["Phi[12]"]= [(1/2*Riccomp[2,3]) ]), adddict(RicNP, ["Phi[21]"]= [(1/2*Riccomp[2,4]) ]), print("Check pt 3"), if sw>0 then ( disp(showdict(RicNP)) ), /*end if*/ remarray(NV,DNV, Weyl04,Riccomp), return ( [[ ["SCNP"]=SCNP, ["WeylNP"]=WeylNP, ["RicNP"]=RicNP, ["RS"]=[RS] ]] ) )$ /*end proc: NPC*/ /******>help */ /*****>proc:procinfo_Riemann */ procinfo([opt]) := block( [choice], if length(opt)=0 then ( choice : "help" ) else ( choice : sdowncase(opt[1]) ), /*end if*/ /****>help **/ if member(choice, ["help"]) then ( printf(true, "Usage: procinfo( arg ). ~%"), printf(true, "arg =\"help\", \"all\" or one of the following function names or key words in the string form. ~%"), printf(true, " \"Riemann\", \"print\", \"Kretschmann\", \"Newmann-Penrose\", \"ConnectionForm\", \"RiemForm\", ~%"), printf(true, " \"Frame\", \"CovDer\". \"IP\", \"Contraction\" ~%"), printf(true, " ~%") ), /*end if*/ /****> General info **/ if member(choice , [ "all", "help"]) then ( printf(true, "* General information,~%"), printf(true, " When you use this package, you have to prepare a coordinate name list, say varlist=~a, ~%", [t,r,theta,phi]), printf(true, " and a metric ds2 expressed as a quadratic equation in the differential of the coordinate names,~%"), printf(true, " say dvarlist=~a and ~%",[dt,dr,dtheta,dphi]), printf(true, " ds2 = ~a~%", -(1-2*M/r)*dt^2+dr^2/(1-2*M/r)+r^2*(dtheta^2+dphi^2*sin(theta)^2)), printf(true, "~%"), /***> Global variables */ printf(true, "* Global variables~%"), printf(true, " _u, _d, p ~%"), printf(true, "~%"), /***> Convention and notations **/ printf(true, "* Variable name convention in this help message~%"), printf(true, " dim = the dimension of the space(time) under consideration.~%"), printf(true, " dict = a special list to represent a dictionary of the following structure:~%"), printf(true, " dict=[ [ [key1]=[data1], [key2]=[data2], ... ] ] ~%"), printf(true, " The contents of a dict can be accessed by various functions. See \"dicttaglist\", \"showdict\", \"ggetdict\", \"adddict\", \"rmdict\", \"repldict\".~%"), printf(true, " T=tensor dict means that the tensor T is expressed by a dict with the following structure: ~%"), printf(true, " tensor dict = dict [ [\"dim\"]=[dim], [\"indextype\"]=indextype, [\"value\"]=Tval, ~%"), printf(true, " indextype= a list specifying the index positions, e.g. [_u, _d, _d, _d] for Riem13 ~%"), printf(true, " Tval = dict [ [ index] = [val], ...]. ~%"), printf(true, " Gdata = a special dict carrying the informaion on the geometry, such as the metric, Christoffel symbol and curvatures. See \"Riemann\" for detail. ~%"), printf(true, " ~%") ), /*end if*/ /****> The central procedure **/ if member(choice, ["all", "riemann", "gdata", "curvature", "ch", "christoffel", "weyl","ricci", "einstein"]) then ( printf(true, "* Central procedures~%"), /***> Riemann */ printf(true, " + Riemann(varlist,ds2[,sw]) => Gdata :: dict ~%"), printf(true, " procedure to calculate the connection coefficients, Riemann curvature tensors, Weyl tensor, ~%"), printf(true, " Ricci tensor, Ricci scalar and Einstein tensor for a given metric. ~%"), printf(true, " sw=0 => no message (default)~%"), printf(true, " sw=1 => diplaying progress messages~%"), printf(true, " sw=2 => displaying some results with progress messages ~%"), printf(true, " * varlist = the coordinate name list.~%"), printf(true, " * dvarlist = the list of the differential of the coordinate names~%"), printf(true, " * ds2 = the metric expressed as an equation quadratic in dvarlist~%"), printf(true, " * Gdata = a dict of the geometric data with keys: ~%"), printf(true, " [\"varlist\"], [\"dim\"], [\"metric\"], [\" Ch\"], [\"Riem13\"], [\"Riem04\"], [\"Riem22\"], ~%"), printf(true, " [\"Weyl04\"], [\"Ric02\"], [\"Ric11\"], [\"RS\"], [\"ET02\"], [\"ET11\"])~%"), printf(true, " Here, the data for each key are~%"), printf(true, " * [\"varlist\"] = varlist : the coordinate variable list, ~%"), printf(true, " * [\"dim\"] = [dim] : the spacetime dimension, ~%"), printf(true, " * [\"metric\"] = [G] : G is the matrix (G[i,j]) representing the spacetime metric,, ~%"), printf(true, " *[\"Ch\"] = Ch: Ch is the dict [ [1,1,1]=[Gamma^1_(11)],...] , ~%"), printf(true, " *[\"Riem13\"]=Riem13: Riem13 is the dict for the Riemann tensor of type (1,3) of the structure [[1,1,1,2]=[Riem^1_(112)],...], ~%"), printf(true, " *[\"Riem04\"]=Riem04 : Riem04 is the dict for the Riemann tensor of type (0,4), ~%"), printf(true, " *[ \"Riem22\"]=Riem04 : Riem04 is the dict for the Riemann tensor of type (2,2), ~%"), printf(true, " *[\"Weyl04\"]=Weyl04 : Weyl04 is the dict for Weyl tensor of type (0,4), ~%"), printf(true, " *[\"Ric02\"]=[Ric02] : Ric02 is the matrix representing the Ricci tensor of type (0,2), ~%"), printf(true, " *[, \"Ric11\"]=[Ric11] : Ric11 is the matrix representing the Ricci tensor of type (1,1), ~%"), printf(true, " *[\"RS\"]=[RS] : RS is the scalar curvature , ~%"), printf(true, " *[ \"ET02\"]=[ET02] : ET02 is the matrix representing the Einstein tensor of type (0,2), ~%"), printf(true, " *[\"ET11\"]=[ET11] : ET11 is the matrix representing the Einstein tensor of type (1,1), ~%"), printf(true, " e.g. ~%"), printf(true, " ds2 : -(1-2*M/r)*dt^2 + dr^2/(1-2*M/r) + r^2*(dtheta^2+sin(theta)^2dphi^2)$ ~%"), printf(true, " varlist : [ t, r, theta, phi]$ ~%"), printf(true, " Gdata : Riemann(varlist, ds2)$ ~%"), printf(true, " Riem13 : getdict(Gdata, [\"Riem13\"])$~%"), printf(true, " factor( getdict(Riem13, [1,2,1,2])[1] ); ~%"), printf(true, " => 2*M/(r^2*(r-2*M)) ~%"), printf(true, " ~%") ), /*end if*/ /***> Printing procedures **/ if member(choice , ["all", "print", "printch", "printriem", "printweyl", "printric", "printet", "printtensor"]) then ( printf(true, "* Printing procedures~%"), printf(true, " The following are procedures to display the data in Gdata on the monitor. ~%"), /* printCh */ printf(true, " + printCh(Gdata) =>sum(sum(Ch[a,b,c]*dx[b]*dx[c],b=1, dim),c=1, dim)~%"), /* printRiem */ printf(true, " + printRiem(Gdata, Ttype) => Riem04/Riem13/Riem22/Weyl04[*,*,*,*]#0 ~%"), printf(true, " Ttype = \"Riem13\", \"Riem04\", \"Riem22\", \"Weyl04\" ~%"), /* printRic */ printf(true, " + printRic(Gdata,Ttype) => Ric02/Ric11/ET02/ET11[*,*] #0 ~%"), printf(true, " Ttype = \"Ric02\", \"Ric11\", \"ET02\", \"ET11\" ~%"), /* printTensor*/ printf(true, " + printTensor(T=tensor dict) => T[*,*,*,...] #0 ~%"), printf(true, " This is a procedure to display all non-vanishing components of a given tensor dict T of any type. ~%") ), /*end if*/ /****> Special procedures**/ /***> Kretschmann/KreschmannS */ if member(choice, ["all", "kretschmann"]) then ( printf(true, " + Kretschmann( Gdata ) => the Kretschmann invariant K ~%"), printf(true, " K= R^(abcd)R_(abcd) ~%"), printf(true, " + KretschmannS( Gdata ) => (K1,K2) [only for a static metric]~%"), printf(true, " K1=4 R^(1i1j)R_(1i1j), K2=R^(ijkl)R_(ijkl) ~% ") ), /*end if*/ /***> Newman-Penrose */ if member(choice, ["all", "npc", "newman-penrose", "np", "spin"]) then ( printf(true, " + NPC(Gdata, NullTetrad[, sw]) => NPCdict [dim=4 only!!]~%"), printf(true, " procedure to calculate the spin coefficients and the Newman-Penrose coefficients ~%"), printf(true, " NullTetrad = list [NB[1]=k_*, NB[2]=l_*, NB[3]=m_*, NB[4]=conjugate(m_*)] ~%"), printf(true, " where k, l, m, conjugate(m) are a linearly independent null 1-forms reprensenting a null tetrad ~%"), printf(true, " normalized as g(k,k)=g(l,l)=g(m,m)=0, g(k.l)=-1, g(m,conjugate(m))=1. ~%"), printf(true, " NPCdict = dict [[\"SCNP\"]=SCNP, [\"WeylNP\"]=WeylNP, [\"RicNP\"]=RicNP, [\"RS\"]=[RS]] ~%"), printf(true, " * SCNP = dict of the NP spin coefficients [[\"alpha\"]=[alpha], ...] ~%"), printf(true, " * WeylNP = dict of the NP coefficients Psi[0], , , Psi[4] for the Weyl tensor [[\"Psi[0]\"]=[Psi[0]], ...] ~%"), printf(true, " * RicNP = dict of the NP coefficients Phi[0,0],, .,Phi[2,2] for the Ricci tensor [[\"Phi[00]\"]=[Phi[00]], ...] ~%"), printf(true, " sw=0 => no output display ~%"), printf(true, " sw=1 => display some results ~%") ), /*end if*/ /***> General frame */ /**> Frame change */ if member(choice, ["all", "frame", "cbtofbt"]) then ( printf(true, " + CBtoFBT(T, dvarlist, FormBasis[,sw]) => T:: tensor dict w.r.t. the frame FormBasis ~%"), printf(true, " transforming the tensor T from the coord basis to a general basis.~%"), printf(true, " T = a tensor dict with components in the coordinate basis~%"), printf(true, " FormBasis = list [FB[1]=theta[1], , ., FB[dim]=theta[m]] ~%"), printf(true, " sw=0 => no message ~%"), printf(true, " sw=1 => with progress messages ~%") ), /*end if*/ /**> Connection form and curvature form */ /* connection form */ if member(choice, ["all", "connection" ,"connectionform", "riemform"]) then ( printf(true, "* Connection and curvature forms~%"), printf(true, " + ConnectionForm(FormBasis, Gdata[, sw]) => CFM=(omega[a,b]) ~%"), printf(true, " CFM = the matrix whose component omega[a,b] is the connection form w.r.t. the FormBasis~%"), printf(true, " FormBasis = [FB[1]=theta[1],, ,FB[dim]=theta[dim]] ~%"), printf(true, " theta[a]= sum(theta[a,mu]*dx[mu], mu=1, dim) ~%"), printf(true, " sw=0 => omega[a,b]=sum(omega[a,b][mu]*dx[mu],mu=1, dim) ~%"), printf(true, " sw=1 => omega[a,b]=sum(omega[a,b][c]*FB[c],c=1, dim) ~%"), /* curvature form */ printf(true, " + RiemForm(FormBasis,Gdata[, sw]) => RF=matrix( RF[a,b]) ~%"), printf(true, " Calculating the curvature form RFa,b] w.r.t. the FormBasis~%"), printf(true, " sw=0 => RF[a,b]=(1/2)sum(sum(Riem13[a,b mu nu]*FB[dx[mu],dx[nu]],mu=1, dim),nu=1, dim) ~%"), printf(true, " sw=1 => RF[a,b]=(1/2)sum(sum(Riem13[a,b c d]*FB[c,d],c=1, dim),d=1, dim) ~%") ), /*end if*/ /**> Covariant derivative of tensors */ if member(choice, ["all", "covder", "covdt", "covndt", "laplacian", "dv", "div"]) then ( printf(true, "* Covariant derivative of tensors ~%"), /*> Covder */ printf(true, " + Covder(T,[a[1],, .,a[n]],[b, c, , .],Gdata)=> (D_a[1] , . D_a[n] T)[b,c,,, .] ~%"), printf(true, " calulating each component of the n-th covariant derivative of a tensor T.~%"), printf(true, " T = a tensor dict or a scalar ~%"), /*> CovDT */ printf(true, " + CovDT(T, Gdata) => DT = the covariant derivative of T~%"), printf(true, " constructing a tensor dict DT correspoinding to the covariant derivative of the given tensor T ~%"), printf(true, " T = a tensor dict or a scalar ~%"), printf(true, " DT = a tensor dict with rank=rank(T)+1 ~%"), /*> CovDnT */ printf(true, " + CovnDT(T, n, Gdata) => DnT ~%"), printf(true, " constructing a tensor dict DnT correspoinding to the n-the covariant derivative of the given tensor T ~%"), printf(true, " T, DnT = a tensor dict ~%"), /*> Laplacian */ printf(true, " + Laplacian(T, [a[1],, .,a[n]], Gdata) => tensor dict (Laplacian T)[a[1],, a[n]] ~%"), printf(true, " T = a tensor dict ~%"), /*> DV */ printf(true, " + DV(X,Y,Gdata) => D_X Y ~%"), printf(true, " calculatingthe covariant derivative of a vector Y wrt a vector X ~%"), printf(true, " X,Y = vectors expressed as lists, vectors, covectors~%"), printf(true, " Output D_X Y is given in the form of list. ~%"), /*> Div */ printf(true, " + Div(T, [a[1],... ], Gdata[,k]) => sum((D_b T)[...,a[k-1],b, a[k], ... ], b=1, dim) ~%"), printf(true, " T = a tensor dict ~%"), printf(true, " k = the sum index position (default value=1) ~%"), /*> DivT */ printf(true, " + DivT(T, Gdata[, k])=> tensor dict DivTval,~%"), printf(true, " T = tensor dict ~%"), printf(true, " k = the sum index position (default value=1) ~%"), printf(true, " DivTval[a[1], ...,]=Div(T, [a[1], ...], Gdata. k) ~%") ), /*end if*/ /***> Algebraic operations */ if member(choice, ["all", "ip", "contractt", "contraction"]) then ( printf(true, "* Algebraic procedures~%"), /**>ContractT */ printf(true, " + ContractT(T, [a,b], Gdata[,sw]) => CT ~%"), printf(true, " constructing a contracted tensor CT with two rank down.~%"), printf(true, " T = tensor dict ~%"), printf(true, " [a,b] = the index positions to contract ~%"), /**> IP */ printf(true, " + IP(V,W, Gdata) => scalar ~%"), printf(true, " calculates the inner product of two tensors/vectors/covectors of the same rank.~%"), printf(true, " V,W = tensor dicts/vectors/covectors~%") ) /*end if*/ )$ /*end proc*/ /******> Applications */ /*****> Job controle */ Job_HIU: false$ Job_SSS: false$ Job_KN: false$ Job_KN_Gdata: false$ /****> Homogeneous and isotropic universe */ if Job_HIU then ( print("Homogeneous and isotropic universe"), /***> metric & curvature */ sf: 'sf, K:'K, ds2_HIU : - dt^2 + sf(t)^2*(dr^2/(1-K*r^2)+ r^2*(dtheta^2+ sin(theta)^2* dphi^2)), varlist_S : [ t, r, theta, phi], Gdata_HIU : Riemann(varlist_S, ds2_HIU), dvarlist_HIU: getdict(Gdata_HIU,["dvarlist"]), Ric11_HIU : mkT2dict( getdict(Gdata_HIU,["Ric11"])[1], [_u,_d]), Riem13_HIU : newdict(["dim"]=[4], ["indextype"]=[_u, _d,_d, _d], ["value"]=getdict(Gdata_HIU,["Riem13"]) ), /* tensordict */ Weyl04_HIU : newdict(["dim"]=[4], ["indextype"]=[_d, _d,_d, _d], ["value"]=getdict(Gdata_HIU,["Weyl04"]) ), /* tensordict */ /***> Orthonormal basis */ ONFB_HIU: [ FB[1]= dt, FB[2]=sf(t)*dr/(1-K*r^2)^(1/2), FB[3]= sf(t)*r*dtheta, FB[4]=sf(t)* r*sin(theta)*dphi ], Riem13_HIU_ONB : CBtoFBT(Riem13_HIU, dvarlist_HIU, ONFB_HIU), /***> Null basis */ NullFB_HIU: [ NFB[2] = dt - sf(t)* dr, NFB[1]= 1/2*(dt+sf(t)*dr), NFB[4]= sf(t)*r*(dtheta - %i *sin(theta)*dphi), NFB[3]=sf(t)* r*(dtheta + %i *sin(theta)*dphi) ], Weyl04_HIU_ONB : CBtoFBT(Weyl04_HIU, dvarlist_HIU, NullFB_HIU), /***> NP coefficients */ NPC_HIU : NPC(Gdata_HIU, NullFB_HIU) )$ /*end Job_HIU*/ /****> Spherically symmetric static spacetime */ if Job_SSS then ( print("Spherically symmetric static spacetime"), /***> metric & curvature */ ff: 'ff, ds2_SSS : -ff(r)*dt^2 + dr^2/ff(r) + r^2*(dtheta^2+ sin(theta)^2* dphi^2), varlist_S : [ t, r, theta, phi], Gdata_SSS : Riemann(varlist_S, ds2_SSS), dvarlist_S: getdict(Gdata_SSS,["dvarlist"]), Ric11_SSS : mkT2dict( getdict(Gdata_SSS,["Ric11"])[1], [_u,_d]), Riem13_SSS : newdict(["dim"]=[4], ["indextype"]=[_u, _d,_d, _d], ["value"]=getdict(Gdata_SSS,["Riem13"]) ), /* tensordict */ Weyl04_SSS : newdict(["dim"]=[4], ["indextype"]=[_d, _d,_d, _d], ["value"]=getdict(Gdata_SSS,["Weyl04"]) ), /* tensordict */ /***> Orthonormal basis */ ONFB_SSS: [ FB[1]= ff(r)^(1/2)*dt, FB[2]=ff(r)^(-1/2)*dr, FB[3]= r*dtheta, FB[4]= r*sin(theta)*dphi ], Riem13_SSS_ONB : CBtoFBT(Riem13_SSS, dvarlist_S, ONFB_SSS), /***> Null basis */ NullFB_SSS: [ NFB[2] = dt - dr/ff(r), NFB[1]= 1/2*(ff(r)*dt+dr), NFB[4]= r*(dtheta - %i *sin(theta)*dphi), NFB[3]= r*(dtheta + %i *sin(theta)*dphi) ], Weyl04_SSS_ONB : CBtoFBT(Weyl04_SSS, dvarlist_S, NullFB_SSS), /***> NP coefficients */ NPC_SSS : NPC(Gdata_SSS, NullFB_SSS) )$ /*end Job_SSS*/ /****> Kerr-Newmann spacetime */ if Job_KN then ( print("Kerr-Newmann spacetime"), /***> metric & curvature */ m: 'm, a: 'a, e:0, Delta : r^2 - 2*m*r+a^2+ e^2, Sigma2: r^2+a^2*cos(theta)^2, Gamma: (r^2+a^2)^2 -a^2*Delta*sin(theta)^2, Omega : a*(2*m*r-e^2)/Gamma, ds2_KN : -Delta*Sigma2/Gamma*dt^2 + Gamma*sin(theta)^2/Sigma2 *(dphi-Omega*dt)^2 + Sigma2*(dr^2/Delta + dtheta^2), display(ds2_KN), varlist_S : [ t, r, theta, phi], if Job_KN_Gdata then ( Gdata_KN : Riemann(varlist_S, ds2_KN,1) ), /*end if*/ dvarlist_S: getdict(Gdata_KN,["dvarlist"]), /* Ric11_KN : mkT2dict( getdict(Gdata_KN,["Ric11"])[1], [_u,_d]), Riem13_KN : newdict(["dim"]=[4], ["indextype"]=[_u, _d,_d, _d], ["value"]=getdict(Gdata_KN,["Riem13"]) ), /* tensordict */ Weyl04_KN : newdict(["dim"]=[4], ["indextype"]=[_d, _d,_d, _d], ["value"]=getdict(Gdata_KN,["Weyl04"]) ), /* tensordict */ print("Calculation of curvatures: done"), */ /***> Null basis: Ker solution (e=0) */ NullFB_KN: [ NFB[1] = -dt + a*sin(theta)^2*dphi+ Sigma2/Delta*dr, NFB[2]= Delta/(2*Sigma2)*(-dt+a*sin(theta)^2*dphi)-dr/2, NFB[3]= 1/(sqrt(2)*(r+%i*a*cos(theta))) *(Sigma2*dtheta-%i*a*sin(theta)*dt + %i*(r^2+a^2)*sin(theta)*dphi), NFB[4]= 1/(sqrt(2)*(r -%i*a*cos(theta))) *(Sigma2*dtheta + %i*a*sin(theta)*dt - %i*(r^2+a^2)*sin(theta)*dphi) ] , /* FBList_KN : map( lambda([x], args(x)[2]), NullFB_KN), eta : blockMatrices([ -matrix([0,1],[1,0]), matrix([0,1],[1,0])]), FBM_KN : coefmatrix( FBList_KN, dvarlist_S), ThetaM_KN : eta . FBM_KN, GI_KN : trigsimp(xthru(getdict(Gdata_KN, ["metric"])[1]^^(-1))), array(NV_KN, 4), array(NF_KN, 4), for i:1 thru 4 do ( NF_KN[i] : row(FBM_KN, i), /*covariant vector*/ NV_KN[i] : transpose(trigsimp(xthru(NF_KN[i] . GI_KN))) /* contravariant vector */ ), /*end do*/ */ /***> NP coefficients */ NPC_KN : NPC(Gdata_KN, NullFB_KN) )$ /*end: Job_KN*/ print("Riemann.mac: load done"); /******> End */