(* There are four rules defining the action of d on forms.  I copy these
straight out of Flanders.  They are a bit abstract, but that makes them ideal
for Mathematica.  The first rule is that d is distributive over addition: *)


        d[w1_ + w2_] := d[w1] + d[w2]

(* The second rule is a wedge product rule: d[w1^w2] = d[w1]^w2 +- w1^d[w2],
where the sign is (-1) to the length of the forms w1.  Formally,
functions are zero forms, so this rule becomes three *)

        d[Wedge[w1_,w2__]] := d[w1] ^ w2 - w1 ^ d[w2]
        d[f_ d[w_]] := d[f] ^ d[w]
        d[f_ Wedge[w1__] ] := d[f] ^ Wedge[w1] + f d[Wedge[w1]]

(* The third rule is that d squared is zero. This is Poincare's lemma. *)

        d[d[w_]] := 0

(* The fourth rule defines d acting on a function to be the gradient.  This
depends on the coordinates we use.  In three dimensions, ... *)

        Clear[x]; Clear[y]; Clear[z]; Clear[dx]; Clear[dy]; Clear[dz]
        d[f_] := D[f,x] d[x] + D[f,y] d[y] + D[f,z] d[z] /; !MemberQ[{x,y,z},f]

(* Finally, we can recognize forms as d[something], so we can pull out
functions from wedge products. *)

(*      Wedge[f_ d[x_],w2_] := f Wedge[d[x],w2]
        Wedge[w1_,f_ d[x_]] := f Wedge[w1,d[x]] 
*)

        Wedge[f_ Wedge[x_,y_],w2_] := f Wedge[Wedge[x,y],w2]
        Wedge[w1_,f_ Wedge[x_,y_]] := f Wedge[w1,Wedge[x,y]]
(* These two lines are modified from the now-commented lines above.
They ensure that things like dx ^ (f dy ^ dz) will be simplified into 
f dx ^ dy ^ dz. The old version doesn't simplify this: dy ^ dz is not a
one d-form and the rule doesn't apply. Now the rule will be applied
recursively until all forms have been grouped together. *)

(* ************************************************************** *)
(* IMPORTANT NOTE *)
(* Note that there's a bug in this package: we can never use Exp[],
because we have stolen the symbol ^ to mean wedge instead of the usual
power, and Mathematica interprets Exp[x] as E^x. *)
(* To be fixed.  Oct 22, 95. *)