/* -*- Mode: MACSYMA -*- */ /* (c) Copyright 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985 Richard Pavelle */ EVAL_WHEN(TRANSLATE,TRANSBIND:TRUE,PACKAGEFILE:TRUE)$ EVAL_WHEN([TRANSLATE,BATCH,DEMO], IF GET('PACKG,'VERSION)=FALSE THEN LOADFILE(PACKG,FASL,DSK,SHAREM))$ HERALD_PACKAGE(CTENSR)$ EVAL_WHEN(TRANSLATE, MODEDECLARE(FUNCTION(SYMMETRICP,DIAGMATRIXP),BOOLEAN))$ TENSORKILL:TRUE$ DEFINE_VARIABLE(DIM,4,FIXNUM,"The dimension of the problem.")$ DEFINE_VARIABLE(DIAGMETRIC,FALSE,BOOLEAN, "True if the metric entered is a diagonal matrix.")$ READVALUE(MESSAGE,PRED,[BADBOY]):= BLOCK([VALUE], LOOP, VALUE:READ(MESSAGE), IF APPLY(PRED,[VALUE])=TRUE THEN RETURN(VALUE), IF BADBOY#[] THEN APPLY('PRINT,BADBOY), GO(LOOP))$ MODEDECLARE(FUNCTION(YESP),BOOLEAN)$ YESP([MESSAGES]):= BLOCK([REPLY], LOOP, REPLY:APPLY('READ,MESSAGES), IF MEMBER(REPLY,['YES,'YE,'Y]) THEN RETURN(TRUE), IF MEMBER(REPLY,['NO,'N]) THEN RETURN(FALSE), PRINT("Please reply YES or NO."), GO(LOOP))$ SWAPOUT(FILE,[FUNCTIONS]):=APPLY('KILL,FUNCTIONS)$ /* Temporary until SWAPOUT is written by GJC */ RESIMP(EXP):=APPLY('EV,[EXP,'NOEVAL,'SIMP])$ /* 62 lines between ^L's can be used with GACHA10 ,68 with GACHA9 */ /* The convention for the sign of the Riemann tensor is the same as in Landau-Lifshits, Misner-Thorne-Wheeler */ /* Kronecker delta function */ KDELT(I,J):=(MODEDECLARE([I,J],FIXNUM),IF I = J THEN 1 ELSE 0)$ /* TTRANSFORM transforms second rank covariant tensors given the components and the transformation law. The latter is in the form YI=YI(X1,...,XDIM) */ TTRANSFORM(QXYZ):=BLOCK([],LOCAL(TTEMP,OMTEMP,VARTEMP,NEWTEMP), FOR I THRU DIM DO (FOR J THRU DIM DO (TTEMP[I,J]:QXYZ[I,J], FOR K THRU DIM DO TTEMP[I,J]:SUBST(OMTEMP[K],OMEGA[K],TTEMP[I,J]))), FOR I THRU DIM DO VARTEMP[I]:READ("Transform #",I), FOR I THRU DIM DO (FOR J THRU DIM DO (FOR K THRU DIM DO TTEMP[I,J]:SUBST(VARTEMP[K],OMTEMP[K],TTEMP[I,J]))), FOR A THRU DIM DO (FOR B THRU DIM DO NEWTEMP[A,B]:TXYZSUM(VARTEMP,OMEGA,A,B,TTEMP)), GENMATRIX(NEWTEMP,DIM,DIM))$ TXYZSUM(VARTEMP,OMEGA,A,B,TTEMP):=BLOCK([TEMP:0],FOR I FROM 1 THRU DIM DO FOR J FROM 1 THRU DIM DO TEMP:TEMP+DIFF(VARTEMP[I],OMEGA[A])* DIFF(VARTEMP[J],OMEGA[B])* TTEMP[I,J],TEMP)$ /* Setup function */ DEFINE_VARIABLE(DERIVABBREV,TRUE,BOOLEAN)$ DEFINE_VARIABLE(TETRADCALEQ,FALSE,BOOLEAN)$ DEFINE_VARIABLE(TAYSWITCH,FALSE,BOOLEAN)$ DEFINE_VARIABLE(RATCHRISTOF,TRUE,BOOLEAN)$ DEFINE_VARIABLE(RATEINSTEIN,TRUE,BOOLEAN)$ DEFINE_VARIABLE(RATRIEMAN,TRUE,BOOLEAN)$ DEFINE_VARIABLE(RATWEYL,TRUE,BOOLEAN)$ SETFLAGS():=(DERIVABBREV:TRUE,TETRADCALEQ:FALSE,TAYSWITCH:FALSE, RATCHRISTOF:TRUE,RATEINSTEIN:TRUE,RATRIEMAN:TRUE,RATWEYL:TRUE)$ TSETUP():=BLOCK([],IF TENSORKILL # TRUE THEN ERROR("Type KILL(ALL)\; and then TENSORKILL:TRUE\; before you enter a new metric."), TENSORKILL:FALSE, SETFLAGS(), DIM:MODE_IDENTITY(FIXNUM, READVALUE("Enter the dimension of the coordinate system:", LAMBDA([V], IF INTEGERP(V) THEN BLOCK([U:V],MODEDECLARE(U,FIXNUM), IF U>0 THEN TRUE)), "Must be a positive integer!")), OMEGA:IF DIM = 2 THEN [X,Y] ELSE (IF DIM = 3 THEN [X,Y,Z] ELSE (IF DIM = 4 THEN [X,Y,Z,T] ELSE READ("Enter a list containing the names of the coordinates in order"))), IF YESP("Do you wish to change the coordinate names?") THEN OMEGA:READ("Enter a list containing the names of the coordinates in order"), IF LENGTH(OMEGA) # DIM THEN ERROR("Length of the coordinate list is not equal to the dimension"), READVALUE("Do you want to 1. Enter a new metric? 2. Enter a metric from a file? 3. Approximate a metric with a Taylor series?", LAMBDA([OPT], IF OPT = 1 THEN(NEWMET(),TRUE) ELSE IF OPT = 2 THEN(FILEMET(),TRUE) ELSE IF OPT = 3 THEN(SERMET(),TRUE) ELSE FALSE), "Invalid option, please enter 1, 2, or 3."), DONE)$ /* Enter a new metric */ NEWMET():=BLOCK([],LG:ENTERMATRIX(DIM,DIM), READ("Enter functional dependencies with the DEPENDS function or 'N' if none"), IF YESP("Do you wish to see the metric?") THEN DISP(LG),METRIC())$ /* Enter a metric from a file */ FILEMET():=BLOCK([FILE,FPOS], FILE:READ("Specify the file as [filename1,filename2,directory]?"), FPOS:(PRINT("What is the ordinal position of the metric in this file?"), READ("(Note, the name LG must be assigned to your metric in the file)")), APPLY(BATCH,[FILE,OFF,FPOS,FPOS]),METRIC())$ /* Approximate a metric with a Taylor series */ SERMET():=BLOCK([],TAYSWITCH:TRUE, PARAM:READ("Enter expansion parameter"), MINP:READ("Enter minimum power to drop"), TAYPT:READ("Enter the point to expand the series around"), IF YESP("Is the metric in a file?") THEN FILEMET() ELSE NEWMET())$ /* Checks diagonality of a matrix or 2D array */ DIAGMATRIXP(MAT,NLIM):=(MODEDECLARE(NLIM,FIXNUM), BLOCK([DIAGFLAG:TRUE], FOR I THRU NLIM DO (FOR J THRU NLIM DO (IF I # J AND MAT[I,J] # 0 THEN RETURN(DIAGFLAG:FALSE))),DIAGFLAG))$ /* Used for Taylor series approximation */ DLGTAYLOR(X):=IF TAYSWITCH THEN RATDISREP(TAYLOR(X,PARAM,TAYPT,MINP-1)) ELSE X$ /* Routine for computing contravariant metric from the covariant but if the metric is diagonal then no out of core files are loaded. UG is defined so that for display purposes it appears with DETOUT and is then evaluated again */ METRIC():=BLOCK([], IF LENGTH(LG) # DIM OR LENGTH(TRANSPOSE(LG)) # DIM THEN ERROR("The rank of the metric is not equal to the dimension of the space") ELSE (IF NOT SYMMETRICP(LG,DIM) THEN ERROR( "You must be working in a new gravity theory not supported by this program")), DIAGMETRIC:DIAGMATRIXP(LG,DIM),GDET:FACTOR(DETERMINANT(LG)), UG:IDENT(LENGTH(FIRST(LG))), IF DIAGMETRIC THEN FOR II:1 THRU LENGTH(FIRST(LG)) DO (UG[II,II]:1/LG[II,II]) ELSE UG:BLOCK([DETOUT:TRUE],LG^^(-1)), IF YESP("Do you wish to see the metric inverse?") THEN LDISP(UG),UG:RESIMP(UG),DONE)$ /* Following returns TRUE if M is an NxN symmetric matrix or array */ SYMMETRICP(M,N):=(MODEDECLARE(N,FIXNUM), BLOCK([SYMFLAG:TRUE], IF N # 1 THEN FOR I THRU N-1 DO (FOR J FROM I+1 THRU N DO (IF M[J,I] # M[I,J] THEN SYMFLAG:FALSE)), SYMFLAG))$ /* Computes geodesic equations of motion */ MOTION(DIS):=BLOCK([S],DEPENDS(OMEGA,S), FOR I THRU DIM DO EM[I]:IF DIAGMETRIC THEN RATSIMP(1/2*SUM( DIFF(LG[A,A],OMEGA[I]) *DIFF(OMEGA[A],S)^2,A,1,DIM)) ELSE 1/2*SUM(SUM(DIFF(LG[A,B],OMEGA[I]) *DIFF(OMEGA[A],S)*DIFF(OMEGA[B],S),A, 1,DIM),B,1,DIM), IF DIS#FALSE THEN FOR I THRU DIM DO LDISPLAY(EM[I]),DONE)$ /* Computes d'Alembertian of a 2nd rank covariant tensor AND DOES NOT WORK FEB. 10, 1980 */ /* Computes covariant and contravariant gradients where the :: allows the user to define the array name*/ COGRAD(F,XYZ):=BLOCK([], FOR MM THRU DIM DO (ARRAYSETAPPLY(XYZ,[MM],RATSIMP(DIFF(F,OMEGA[MM])))))$ CONTRAGRAD(F,XYZ):=BLOCK([], IF DIAGMETRIC THEN FOR MM THRU DIM DO (ARRAYSETAPPLY(XYZ,[MM],RATSIMP(RATSIMP(UG[MM,MM]*DIFF(F,OMEGA[MM]))))) ELSE FOR MM THRU DIM DO (ARRAYSETAPPLY(XYZ,[MM],SUM(UG[M,N]*DIFF(F,OMEGA[NN]),NN,1,DIM))))$ /* DALEM(P,I,J):=IF DIAGMETRIC THEN RATSIMP(SUM(CONTRAGRAD(COGRAD(P[I,J],M),M),M,1,DIM) +1/SQRT(-GDET)*SUM(COGRAD(SQRT(-GDET)*UG[M,M],M)*COGRAD(P[I,J],M),M,1,DIM)) ELSE SUM(CONTRAGRAD(COGRAD(P[I,J],M),M),M,1,DIM) +1/SQRT(-GDET)*SUM(SUM(COGRAD(SQRT(-GDET)*UG[M,N],N)*COGRAD(P[I,J],M),N,1, DIM),M,1,DIM)$ */ /* Routine for computing the Christoffel symbols COMMENT: Change routine so GDET is not evaluated until the end */ CHRISTOF(DIS):=BLOCK([], FOR I THRU DIM DO (FOR J FROM I THRU DIM DO (FOR K THRU DIM DO LCS[J,I,K]:LCS[I,J,K] :(DIFF(LG[J,K],OMEGA[I]) +DIFF(LG[I,K],OMEGA[J]) -DIFF(LG[I,J],OMEGA[K]))/2, FOR K THRU DIM DO (MCS[J,I,K]:MCS[I,J,K] :DLGTAYLOR( IF DIAGMETRIC THEN RATSIMP(UG[K,K]*LCS[I,J,K]) ELSE SUM(UG[K,L]*LCS[I,J,L],L,1,DIM)), IF RATCHRISTOF THEN MCS[J,I,K] :MCS[I,J,K]:RATSIMP(MCS[I,J,K])))), IF DIS = ALL OR DIS = LCS THEN FOR I THRU DIM DO (FOR J:I THRU DIM DO (FOR K THRU DIM DO (IF LCS[I,J,K] # 0 THEN LDISPLAY(LCS[I,J,K])))), REMARRAY(LCS), IF DIS = MCS OR DIS = ALL THEN FOR I THRU DIM DO (FOR J:I THRU DIM DO (FOR K THRU DIM DO (IF MCS[I,J,K] # 0 THEN LDISPLAY(MCS[I,J,K])))),DONE)$ /* Covariant components of the Ricci tensor */ LRICCICOM(DIS):=BLOCK([SUMA,SUMB,FLAT:TRUE], MODEDECLARE(FLAT,BOOLEAN), IF MEMBER('MCS,ARRAYS) THEN (TRUE) ELSE (CHRISTOF(FALSE)), FOR I THRU DIM DO (FOR J FROM I THRU DIM DO (SUMA:0,SUMB:0, FOR K THRU DIM DO (IF K # I THEN SUMA:SUMA +(DIFF(MCS[I,J,K],OMEGA[K]) -DIFF(MCS[J,K,K],OMEGA[I])), SUMB:SUMB+SUM( MCS[K,L,K]*MCS[I,J,L] -MCS[I,K,L]*MCS[J,L,K],L,1,DIM)), LR[I,J]:DLGTAYLOR(SUMA+SUMB), IF RATFAC THEN LR[I,J]:FACTOR(DLGTAYLOR(RATSIMP(LR[I,J]))), LR[J,I]:LR[I,J])), IF DIS#FALSE THEN (FOR I THRU DIM DO (FOR J FROM I THRU DIM DO (IF LR[I,J] # 0 THEN (FLAT:FALSE,LDISPLAY(LR[I,J])))), IF FLAT THEN PRINT("This spacetime is empty and/or flat")), TRACER:IF DIAGMETRIC THEN DLGTAYLOR(RATSIMP(SUM(LR[I,I]*UG[I,I],I,1,DIM))) ELSE DLGTAYLOR(SUM(SUM(LR[I,J]*UG[I,J],I,1,DIM),J,1,DIM)),DONE)$ /* Forms a scalar from the inner product of the two arguments assumed to be second rank matrices */ CONTRACT(BLU,BLA):= BLOCK(MODEDECLARE([I,J],FIXNUM), RATSIMP(DLGTAYLOR(SUM(SUM(BLU[I,J]*BLA[I,J],I,1,DIM),J,1,DIM)),DONE))$ /* Computes mixed Ricci tensor where the first index is covariant */ RICCICOM(DIS):=BLOCK([FLAT:TRUE],MODEDECLARE(FLAT,BOOLEAN), IF MEMBER('LR,ARRAYS) THEN (TRUE) ELSE (LRICCICOM(FALSE)), FOR I THRU DIM DO FOR J THRU DIM DO (RICCI[I,J]:0), FOR I THRU DIM DO (FOR J THRU DIM DO (RICCI[I,J]:IF DIAGMETRIC THEN RATSIMP(DLGTAYLOR(LR[I,J]*UG[J,J])) ELSE DLGTAYLOR(SUM(LR[I,K]*UG[K,J],K,1,DIM)), IF RATFAC THEN RICCI[I,J] :FACTOR(DLGTAYLOR(RATSIMP(RICCI[I,J]))))), IF DIS#FALSE THEN (FOR I THRU DIM DO (FOR J THRU DIM DO (IF RICCI[I,J] # 0 THEN (FLAT:FALSE,LDISPLAY(RICCI[I,J])))), IF FLAT THEN PRINT("This spacetime is empty and/or flat")), DONE)$ /* Computes scalar curvature */ SCURVATURE():=IF RATFAC THEN FACTOR(TRACER) ELSE TRACER$ /* Computes mixed Einstein tensor where the first index is covariant */ EINSTEIN(DIS):=BLOCK([FLAT:TRUE],MODEDECLARE(FLAT,BOOLEAN), IF MEMBER('RICCI,ARRAYS) THEN (TRUE) ELSE (RICCICOM(FALSE)), FOR I THRU DIM DO (FOR J THRU DIM DO (IF I = J THEN G[I,J]:DLGTAYLOR(RICCI[I,J]-TRACER/2) ELSE G[I,J]:DLGTAYLOR(RICCI[I,J]), IF RATFAC THEN G[I,J]:FACTOR(RATSIMP(G[I,J])) ELSE (IF RATEINSTEIN THEN G[I,J]:RATSIMP(G[I,J])), IF DIS#FALSE AND G[I,J] # 0 THEN (FLAT:FALSE,LDISPLAY(G[I,J])))), IF DIS#FALSE AND FLAT THEN PRINT("This spacetime is empty and/or flat"),DONE)$ /* Computes covariant Riemann curvature tensor */ RIEMANN(DIS):=BLOCK([FLAT:TRUE], MODEDECLARE(FLAT,BOOLEAN), IF MEMBER('MCS,ARRAYS) THEN (TRUE) ELSE (CHRISTOF(FALSE)), FOR I THRU DIM DO (FOR J THRU DIM DO (FOR K THRU DIM DO (FOR L THRU DIM DO R[I,J,K,L]:0))), FOR I THRU DIM DO (FOR J FROM I+1 THRU DIM DO (FOR K FROM I THRU DIM DO (FOR L FROM K+1 THRU (IF K = I THEN J ELSE DIM) DO (R[I,J,K,L]:DLGTAYLOR( 1/2 *(DIFF(LG[I,L],OMEGA[J],1,OMEGA[K],1) +DIFF(LG[J,K],OMEGA[I],1,OMEGA[L],1) -DIFF(LG[I,K],OMEGA[J],1,OMEGA[L],1) -DIFF(LG[J,L],OMEGA[I],1,OMEGA[K],1)) +(IF DIAGMETRIC THEN RATSIMP( SUM( LG[U,U] *(MCS[J,K,U]*MCS[I,L,U] -MCS[J,L,U]*MCS[I,K,U]),U,1, DIM)) ELSE SUM( SUM( LG[U,V] *(MCS[J,K,U]*MCS[I,L,V] -MCS[J,L,U]*MCS[I,K,V]),V,1, DIM),U,1,DIM))), IF RATFAC THEN R[I,J,K,L] :FACTOR(RATSIMP(R[I,J,K,L])) ELSE (IF RATRIEMAN THEN R[I,J,K,L] :RATSIMP(R[I,J,K,L])), R[J,I,L,K]:R[I,J,K,L], R[I,J,L,K]:R[J,I,K,L]:-R[I,J,K,L], IF I # K OR J > L THEN (R[K,L,I,J]:R[L,K,J,I]:R[I,J,K,L], R[L,K,I,J]:R[K,L,J,I]:-R[I,J,K,L]), IF DIS#FALSE AND R[I,J,K,L] # 0 THEN (FLAT:FALSE,LDISPLAY(R[I,J,K,L])))))), IF DIS#FALSE AND FLAT THEN PRINT("This spacetime is flat"),DONE)$ /* Computes contravariant Riemann tensor from covariant form */ RAISERIEMANN(DIS):=BLOCK([], IF MEMBER('R,ARRAYS) THEN (TRUE) ELSE (RIEMANN(FALSE)), FOR I THRU DIM DO (FOR J THRU DIM DO (FOR K THRU DIM DO (FOR L THRU DIM DO UR[I,J,K,L]:0))), FOR I THRU DIM DO (FOR J FROM I+1 THRU DIM DO (FOR K FROM I THRU DIM DO (FOR L FROM K+1 THRU (IF K = I THEN J ELSE DIM) DO (UR[I,J,K,L] :IF DIAGMETRIC THEN RATSIMP( R[I,J,K,L]*UG[I,I]*UG[J,J]*UG[K,K] *UG[L,L]) ELSE SUM( SUM( SUM( SUM( R[A,B,C,D]*UG[I,A]*UG[J,B]*UG[K,C] *UG[L,D],A,1,DIM),B,1, DIM),C,1,DIM),D,1,DIM), IF RATRIEMAN THEN UR[I,J,K,L]:RATSIMP(UR[I,J,K,L]), UR[J,I,L,K]:UR[I,J,K,L], UR[I,J,L,K]:UR[J,I,K,L]:-UR[I,J,K,L], IF I # K OR J > L THEN ( UR[K,L,I,J]:UR[L,K,J,I]:UR[I,J,K,L], UR[L,K,I,J]:UR[K,L,J,I]:-UR[I,J,K,L]), IF DIS#FALSE AND UR[I,J,K,L] # 0 THEN LDISPLAY(UR[I,J,K,L]))))))$ /* Computes the Kretchmann invariant R[i,j,k,l]*UR[i,j,k,l] */ /* old definition RINVARIANT():=KINVARIANT:SUM( SUM(SUM(SUM(R[I,J,K,L]*UR[I,J,K,L],I,1,DIM),J, 1,DIM),K,1,DIM),L,1,DIM)$ */ RINVARIANT():=IF MEMBER('UR,ARRAYS) THEN KINVARIANT:SUM( SUM(SUM(SUM(R[I,J,K,L]*UR[I,J,K,L],I,1,DIM),J, 1,DIM),K,1,DIM),L,1,DIM) ELSE (RAISERIEMANN(FALSE),KINVARIANT:SUM( SUM(SUM(SUM(R[I,J,K,L]*UR[I,J,K,L],I,1,DIM),J, 1,DIM),K,1,DIM),L,1,DIM))$ /* Computes covariant Weyl tensor */ WEYL(DIS):=BLOCK([FLAT:TRUE],MODEDECLARE(FLAT,BOOLEAN), IF DIM = 2 THEN (PRINT("All 2 dimensional spacetimes are conformally flat"), RETURN(DONE)), IF (MEMBER('LR,ARRAYS),MEMBER('R,ARRAYS)) THEN TRUE ELSE (LRICCICOM(FALSE),RIEMANN(FALSE)), FOR I THRU DIM DO (FOR J THRU DIM DO (FOR K THRU DIM DO (FOR L THRU DIM DO W[I,J,K,L]:0))), FOR I THRU DIM DO (FOR J FROM I+1 THRU DIM DO (FOR K FROM I THRU DIM DO (FOR L FROM K+1 THRU (IF K = I THEN J ELSE DIM) DO (W[I,J,K,L]:DLGTAYLOR( R[I,J,K,L] +1/(DIM-2) *(LG[I,L]*LR[J,K]-LG[I,K]*LR[L,J] +LG[J,K]*LR[L,I] -LG[J,L]*LR[K,I]) -TRACER/((DIM-1)*(DIM-2)) *(LG[I,L]*LG[K,J]-LG[I,K]*LG[L,J])), IF RATFAC THEN W[I,J,K,L] :FACTOR(RATSIMP(W[I,J,K,L])) ELSE (IF RATWEYL THEN W[I,J,K,L]:RATSIMP(W[I,J,K,L])), W[J,I,L,K]:W[I,J,K,L], W[I,J,L,K]:W[J,I,K,L]:-W[I,J,K,L], IF I # K OR J > L THEN (W[K,L,I,J]:W[L,K,J,I]:W[I,J,K,L], W[L,K,I,J]:W[K,L,J,I]:-W[I,J,K,L]), IF DIS#FALSE AND W[I,J,K,L] # 0 THEN (FLAT:FALSE,LDISPLAY(W[I,J,K,L])))))), IF DIS#FALSE AND FLAT THEN PRINT("This spacetime is conformally flat"),DONE)$ /* Number of terms per component of the array F */ NTERMST(F):=FOR I THRU DIM DO (FOR J THRU DIM DO PRINT([[I,J],NTERMS(F[I,J])]))$ /* Computes d'Alembertian of the scalar PHI */ DSCALAR(PHI):=IF DIAGMETRIC THEN RATSIMP(1/SQRT(-GDET)*SUM( DIFF( UG[I,I]*SQRT(-GDET)*DIFF(PHI,OMEGA[I]), OMEGA[I]),I,1,DIM)) ELSE RATSIMP(1/SQRT(-GDET)*SUM( SUM( DIFF( UG[I,J]*SQRT(-GDET)*DIFF(PHI,OMEGA[J]), OMEGA[I]),I,1,DIM),J,1,DIM))$ /* Computes the covariant divergence of the object GXYZ which must be a mixed 2nd rank tensor whose first index is covariant- a check should be put in to see if GXYZ has the correct dimensions Apr.9,80 */ CHECKDIV(GXYZ):=BLOCK(MODEDECLARE([I,J],FIXNUM), IF DIAGMATRIXP(GXYZ,DIM) THEN FOR I THRU DIM DO PRINT(DIV[I]:RATSIMP(DLGTAYLOR(1/SQRT(-GDET) *DIFF(SQRT(-GDET)*GXYZ[I,I],OMEGA[I]) -SUM(MCS[I,J,J]*GXYZ[J,J],J,1,DIM)))) ELSE FOR I THRU DIM DO PRINT(DIV[I]:RATSIMP(DLGTAYLOR(1/SQRT(-GDET) *SUM(DIFF(SQRT(-GDET)*GXYZ[I,J],OMEGA[J]),J,1,DIM) -SUM(SUM(MCS[I,J,A]*GXYZ[A,J],A,1,DIM),J,1,DIM)))))$ /* FINDDE returns a list of the unique differential equations in the list or 2 or 3 dimensional array A. DEINDEX is a global list containing the indices of A corresponding to these unique differential equations. */ FINDDE1(LIST):=BLOCK([INFLAG:TRUE,DERIV:NOUNIFY('DERIVATIVE),L:[],L1,Q], DEINDEX:[], FOR I WHILE LIST # [] DO (L1:FACTOR(NUM(FIRST(LIST))),LIST:REST(LIST), IF NOT FREEOF(DERIV,L1) THEN (DEINDEX:CONS(I,DEINDEX), IF INPART(L1,0) # "*" THEN L:CONS(L1,L) ELSE (Q:1, FOR J THRU LENGTH(L1) DO (IF NOT FREEOF(DERIV,INPART(L1,J)) THEN Q:Q*INPART(L1,J)), L:CONS(Q,L)))), CLEANUP(L))$ FINDDE2(A):=BLOCK([INFLAG:TRUE,DERIV:NOUNIFY('DERIVATIVE),L:[],T,Q], DEINDEX:[], FOR I THRU DIM DO (FOR J THRU DIM DO (T:FACTOR(NUM(A[I,J])), IF NOT FREEOF(DERIV,T) THEN (DEINDEX:CONS([I,J],DEINDEX), IF INPART(T,0) # "*" THEN L:CONS(T,L) ELSE (Q:1, FOR N THRU LENGTH(T) DO (IF NOT FREEOF(DERIV,INPART(T,N)) THEN Q:Q*INPART(T,N)), L:CONS(Q,L))))), CLEANUP(L))$ FINDDE3(A):=BLOCK([INFLAG:TRUE,DERIV:NOUNIFY('DERIVATIVE),L:[],T,Q], DEINDEX:[], FOR I THRU DIM DO (FOR J THRU DIM DO (FOR K THRU DIM DO (T:FACTOR(NUM(A[I,J,K])), IF NOT FREEOF(DERIV,T) THEN (DEINDEX:CONS([I,J,K],DEINDEX), IF INPART(T,0) # "*" THEN L:CONS(T,L) ELSE (Q:1, FOR N THRU LENGTH(T) DO (IF NOT FREEOF(DERIV,INPART(T,N)) THEN Q:Q*INPART(T,N)), L:CONS(Q,L)))))), CLEANUP(L))$ CLEANUP(LL):=BLOCK([A,L:[],INDEX:[]], WHILE LL # [] DO (A:FIRST(LL),LL:REST(LL), IF NOT MEMBER(A,LL) THEN (L:CONS(A,L),INDEX:CONS(FIRST(DEINDEX),INDEX)), DEINDEX:REST(DEINDEX)),DEINDEX:INDEX,L)$ FINDDE(A,N):=(MODEDECLARE(N,FIXNUM), IF N = 1 THEN FINDDE1(A) ELSE (IF N = 2 THEN FINDDE2(A) ELSE (IF N = 3 THEN FINDDE3(A) ELSE ERROR("Invalid dimension:",N))))$ DELETEN(L,N):=(MODEDECLARE(N,FIXNUM), BLOCK([LEN],MODEDECLARE(LEN,FIXNUM), IF L = [] THEN L ELSE (LEN:LENGTH(L), IF N > LEN OR N < 1 THEN ERROR("Second argument out of range") ELSE (IF N = 1 THEN REST(L) ELSE (IF N = LEN THEN REST(L,-1) ELSE APPEND(REST(L,N-LEN-1), REST(L,N)))))))$ GETRID():= (MAPLIST(LAMBDA([U],U::FALSE), '[TAYSWITCH,RATCHRISTOF,EXPANDCHRISTOF,RATEINSTEIN,RATRIEMAN, HALFC,CHRSUB,MOTION,DLGTAYLOR,TSETUP,CHRISTOF,RICCICOM,TESTINDEX,EINSTEIN, TTRANSFORM,RIEMANN,DIAGMATRIXP,RAISERIEMANN,RSCALAR,LRICCICOM,WEYL,METRIC]), SWAPOUT([], TETRADCALEQ,RATWEYL,NICEPRINT,KDELT,SETFLAGS,BDVAC,INVARIANT1,INVARIANT2, TSETUP,NEWMET,FILEMET,SERMET,SYMMETRICP,DL,DU,DALEM,SCURVATURE,RINVARIANT, NTERMST,DSCALAR,CHECKDIV,SETUPTETRAD,CONTRACT4,PSI,PETROV,FINDDE1,FINDDE2, FINDDE3,CLEANUP,FINDDE,DELETEN,GETRID))$ /* The 4 dimensional Brans- Dicke vacuum field equations are represented by the array BD and the scalar equation is generated by setting the d'Alembertian of the scalar field to zero. That is one calls the function DSCALAR on the scalar field. */ BDVAC():=BLOCK([], IF DIM # 4 THEN ERROR(" THIS PROGRAM IS RESTRICTED TO 4 DIMENSIONS"), ZZ:READ("give a name to the scalar field and declare its functional dependencies"), BOXQ:0, FOR I:1 THRU 4 DO FOR J:1 THRU 4 DO (ADDD[I,J]: W/ZZ^2*( DIFF(ZZ,OMEGA[I])*DIFF(ZZ,OMEGA[J])- LG[I,J]*SUM(DIFF(ZZ,OMEGA[KK])*DIFF(ZZ,OMEGA[KK])*UG[KK,KK],KK,1,4)/2)+ (DIFF(DIFF(ZZ,OMEGA[I]),OMEGA[J])-SUM(MCS[I,J,KK]*DIFF(ZZ,OMEGA[KK]),KK,1,4)- LG[I,J]*BOXQ)/ZZ), FOR I:1 THRU 4 DO FOR J:1 THRU 4 DO (BD[I,J]:RATSIMP( LR[I,J]-R*LG[I,J]/2-0*T[I,J]-ADDD[I,J])), REMARRAY(ADDD))$ /* This gives the Euler- Lagrange equations for the density of the invariant R^2. The form is (where D is the Kronecker delta) b 2 b b bc (1/2)*D R - 2 R R + 2 D []R -2 g R a a a ;ac The equations are sometimes less complex if TRACER is given a parametric name with dependencies corresponding to those of the Scalar Curvature */ INVARIANT1():=BLOCK([], FOR AA THRU DIM DO FOR BB THRU DIM DO (INV1[AA,BB]:0), IF DIAGMETRIC THEN FOR AA THRU DIM DO FOR BB THRU DIM DO (INV1[AA,BB]:RATSIMP( KDELT(AA,BB)*TRACER^2/2- 2*TRACER*RICCI[AA,BB]- 2*KDELT(AA,BB)*DSCALAR(TRACER)+ 2*UG[AA,AA]*( DIFF(DIFF(TRACER,OMEGA[BB]),OMEGA[AA]) -SUM(MCS[BB,AA,KK]*DIFF(TRACER,OMEGA[KK]),KK,1,DIM)))) ELSE FOR AA THRU DIM DO FOR BB THRU DIM DO (INV1[AA,BB]:RATSIMP( KDELT(AA,BB)*TRACER^2/2- 2*TRACER*RICCI[AA,BB]- 2*KDELT(AA,BB)*DSCALAR(TRACER)+ 2*SUM(UG[BB,CC]*( DIFF(DIFF(TRACER,OMEGA[AA]),OMEGA[CC]) -SUM(MCS[AA,CC,KK]*DIFF(TRACER,OMEGA[KK]),KK,1,DIM)),CC,1,DIM))))$ INVARIANT2():="NOT YET IMPLEMENTED"; BIMETRIC():="NOT YET IMPLEMENTED";