/*
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of
 * the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be
 * useful, but WITHOUT ANY WARRANTY; without even the implied
 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 * PURPOSE.  See the GNU General Public License for more details.
 *
 * Maxwell's equations
 */

if get('itensor,'version)=false then load(itensor);
idim(4);
imetric(g);
/*dim:4;
("We need to declare some properties of the Levi-Civita symbol")$
declare(nounify(levi_civita),constant);
decsym(nounify(levi_civita),4,0,[anti(all)],[]);
decsym(nounify(levi_civita),0,4,[],[anti(all)]);*/
decsym(g,2,0,[sym(all)],[]);
decsym(g,0,2,[],[sym(all)]);
("The EM field tensor is the exterior derivative of a vector field")$
remcomps(F);
components(F([i,j],[]),extdiff(A([i],[]),j));
components(F([],[i,j]),F([k,l],[])*g([],[i,k])*g([],[j,l]));
showcomps(F([i,j],[]));
("Two of Maxwell's equations reduce to simple geometric identities")$
ishow('extdiff('F([i,j],[]),k)=extdiff(F([i,j],[]),k))$
("The other two Maxwell equations define the 4-current")$
remcomps(J);
components(J([],[i]),idiff(F([],[i,j]),j));
("We limit ourselves to the constant metric of special relativity...")$
declare(g,constant);
("Here, the 4-current satisfies a conservation equation")$
covdiff(J([],[i]),i)$
%,ichr2$
undiff(%)$
%,idiff$
canform(%)$
rename(%);
("The EM field tensor is invariant under a gauge transformation.")$
("We can add the gradient of an arbitrary scalar function to A:")$
extdiff(A([i],[])+f([],[],i),j)-F([i,j],[]);
("The dual of the EM field tensor shall be denoted xF here:")$
remcomps(xF);
components(xF([a,b],[]),'levi_civita([a,b,c,d],[])*F([],[c,d]));
("We also define mixed-index components")$
components(F([a],[b]),F([a,c],[])*g([],[b,c]));
components(xF([a],[b]),xF([a,c],[])*g([],[b,c]));
("This is what we need to define the energy-momentum tensor")$
remcomps(T);
components(T([a,b],[]),(F([a,c],[])*F([b],[c])+xF([a,c],[])*xF([b],[c]))/8/%pi);
components(T([a],[b]),T([a,c],[])*g([],[b,c]));
components(T([],[a,b]),T([c],[b])*g([],[a,c]));
("T satisfies the equation T_[a,b;c] = 0")$
covdiff(canform(T([a,b],[])),c)-
covdiff(canform(T([b,a],[])),c)+
covdiff(canform(T([b,c],[])),a)-
covdiff(canform(T([c,b],[])),a)+
covdiff(canform(T([c,a],[])),b)-
covdiff(canform(T([a,c],[])),b)$
%,ichr2$
undiff(%)$
%,idiff$
canform(%);
("The energy-momentum tensor is symmetric in its indices")$
remove(g,constant);
T([a,b],[])-T([b,a],[])$
canform(lc2kdt(%))$
%,kdelta$
canform(rename(contract(expand(%))));
("Let us examine the 4-current again")$
declare(g,constant);
rename(J([],[a]))$
undiff(%)$
%,idiff$
ishow(rename(%))$
("In empty space, the 4-current is identically zero.")$
("This can only happen if A([a],[],b,c)-A([b],[],a,c)=0.")$
("In this case, the energy-momentum tensor should be conserved")$
covdiff(T([a],[b]),b)$
%,ichr2$
undiff(%)$
%,idiff$
canform(lc2kdt(%))$
%,kdelta$
canform(%)$
ishow(canform(rename(contract(expand(%)))))$
("Let us apply what we learned about the 4-current:")$
subst(A([%2],[],a,%1),A([a],[],%1,%2),%th(2))$
canform(%)$
ishow(contract(%))$
("and apply it again, resolving subexpressions 'by hand'")$
ratsubst(A([],[%3],%2,a)*g([],[%1,%2]),A([],[%1],%2,a)*g([],[%2,%3]),%th(2))$
ishow(canform(contract(expand(%))))$
ratsubst(A([],[%3],%2,%3)*g([],[%1,%2]),A([],[%1],%2,%3)*g([],[%2,%3]),%)$
ishow(canform(contract(expand(%))))$
subst(A([%2],[],%1,a),A([%1],[],%2,a),%)$
subst(A([%2],[],%1,%3),A([%1],[],%2,%3),%)$
("and we find that energy is, indeed, conserved for a free EM field:")$
canform(contract(%th(2)));

/* End of demo -- comment line needed by MAXIMA to resume demo menu */