#-
*   Title:  Retarded/advanced derivative declarations
*           and procedures
*    File:  derivs-ra.h
*  Author:  David N. Williams
* License:  Creative Commons Attribution-Share Alike
* Started:  August 27, 2012
* Revised:  August 27-31, 2012
*           September 1-3,5,8,10, 2012
*
* All parts of this program not in the public domain are:
*
*   Copyright (C) 2012 David N. Williams
*
* This work is licensed under the Creative Commons Attribution-
* Share Alike 2.5 License.  To view a copy of this license, visit
* http://creativecommons.org/licenses/by-sa/2.5/ or send a letter
* to Creative Commons, 543 Howard Street, 5th Floor, San
* Francisco, California, 94105, USA.

S c,f,f1,f2,f3,[1/r.u];
V r,u,a,j,jdot;
CF g(symmetric);

* for differentiation procedures (up to three derivs)
CF [f],[f1],[f2],[f3],[r.u^-1];
CF [df],[df1],[df2],[df3],[dr.u^-1];
CF [r],[u],[a],[j];
CF [dr],[du],[da],[dj];
F  [Nf],[Nf1],[Nf2],[Nf3],[Nr.u^-1];
F  [Ndf],[Ndf1],[Ndf2],[Ndf3],[Ndr.u^-1];
F  [Nr],[Nu],[Na],[Nj],dtau;
F  [Ndr],[Ndu],[Nda],[Ndj];

* dummies for exclusive use by procedures
S  dumi;
I  dum1,dum2,dum3;
CF [s],[ds],[v],[dv];

Set scal:    f,    f1,    f2,    f3,    [1/r.u];
Set cscal:   [f],  [f1],  [f2],  [f3],  [r.u^-1];
Set ncscal:  [Nf], [Nf1], [Nf2], [Nf3], [Nr.u^-1];
Set cdscal:  [df], [df1], [df2], [df3], [dr.u^-1];
Set ncdscal: [Ndf],[Ndf1],[Ndf2],[Ndf3],[Ndr.u^-1];

Set vec:    r,    u,    a,    j;
Set cvec:   [r],  [u],  [a],  [j];
Set ncvec:  [Nr], [Nu], [Na], [Nj];
Set cdvec:  [dr], [du], [da], [dj];
Set ncdvec: [Ndr],[Ndu],[Nda],[Ndj];

* The partial differentiation procedure is based on André Heck's
* diff procedure example from "FORM for Pedestrians", revised
* for tensor products of vectors.  We have removed his
* differential place holder dx from the argument in diff, and
* replaced diff(x,dx) by a procedure partial(mu), where the
* differential place holder is hard-wired as dtau.

#procedure partial(mu)
* Flatten powers of [r.u^-1] left by previous diff's.
  Id [r.u^-1]^2 = [r.u^-1]*[r.u^-1];

* Replace commuting functions with noncommuting functions.
  Id [s]?cscal?ncscal = [s];
  Id [v]?cvec?ncvec(mu?) = [v](mu);

* Initialize differentiation.
  Multiply left dtau;

* Commute dtau and leave formal derivs of scalars.
  Repeat;
    Id dtau*[s]?ncscal[dumi] = ncdscal[dumi](`mu') + [s]*dtau;
  Endrepeat;

* Commute dtau and leave formal derivs of vectors.
  Repeat;
    Id dtau*[v]?ncvec[dumi](dum1?) = ncdvec[dumi](`mu',dum1) + [v](dum1)*dtau;
  Endrepeat;

* Remove leftover terms with dtau on the right.
  Id dtau = 0;

* Replace noncommuting functions with commuting functions.
  Id [s]?ncscal?cscal = [s];
  Id [ds]?ncdscal?cdscal(dum1?) = [ds](dum1);
  Id [v]?ncvec?cvec(dum1?) = [v](dum1);
  Id [dv]?ncdvec?cdvec(dum1?,dum2?) = [dv](dum1,dum2);

* Expand primitive formal derivs.
  Id [dr](dum1?,dum2?) =  g(dum1,dum2) - [r](dum1)*[u](dum2)*[r.u^-1];
  Id [du](dum1?,dum2?) = [r](dum1)*[a](dum2)*[r.u^-1];
  Id [da](dum1?,dum2?) = [r](dum1)*[j](dum2)*[r.u^-1];
  Id [dj](dum1?,dum2?) = [r](dum1)*jdot(dum2)*[r.u^-1];
  Id [df](dum1?)  = [f1]*( [u](dum1)
    + [r](dum1)*( [r](dum2)*[a](dum2) - c^2 )*[r.u^-1] );
  Sum dum2;
  Id [df1](dum1?) = [f2]*( [u](dum1)
    + [r](dum1)*( [r](dum2)*[a](dum2) - c^2 )*[r.u^-1] );
  Sum dum2;
  Id [df2](dum1?) = [f3]*( [u](dum1)
    + [r](dum1)*( [r](dum2)*[a](dum2) - c^2 )*[r.u^-1] );
  Sum dum2;
  Id [dr.u^-1](dum1?) = -[r.u^-1]^2*( [u](dum1)
    + [r](dum1)*( [r](dum2)*[a](dum2) - c^2 )*[r.u^-1] );
  Sum dum2;
#endprocedure

#procedure idvec
* Replace functions with vectors and symbols.
  Repeat;
    Id [v]?cvec[dumi](dum1?) = vec[dumi](dum1);
  Endrepeat;
  Repeat;
    Id [s]?cscal[dumi] = scal[dumi];
  Endrepeat;
  Id [1/r.u] = r.u^-1;
#endprocedure

#procedure idmetric
  Id g(dum1?,dum2?)*g(dum1?,dum3?) = g(dum2,dum3);
  Id g(dum1?,dum1?) = 4;
  Id g(dum1?,u?) = u(dum1);
  Id r.r = 0;
  Id u.u = c^2;
  Id u.a = 0;
  Id u.j = -a.a;
#endprocedure2

#procedure idf
  Id f  =    r.u^-1;
  Id f1 =   -r.u^-2;
  Id f2 =  2*r.u^-3;
  Id f3 = -6*r.u^-4;
#endprocedure

#procedure idfp2
  Id f  =     r.u^-2;
  Id f1 =  -2*r.u^-3;
  Id f2 =   6*r.u^-4;
  Id f3 = -24*r.u^-5;
#endprocedure

#procedure idfp3
  Id f  =     r.u^-3;
  Id f1 =  -3*r.u^-4;
  Id f2 =  12*r.u^-5;
  Id f3 = -60*r.u^-6;
#endprocedure

#procedure idfp4
  Id f  =      r.u^-4;
  Id f1 =   -4*r.u^-5;
  Id f2 =   20*r.u^-6;
  Id f3 = -120*r.u^-7;
#endprocedure

#procedure idfp5
  Id f  =      r.u^-5;
  Id f1 =   -5*r.u^-6;
  Id f2 =   30*r.u^-7;
  Id f3 = -210*r.u^-8;
#endprocedure

#procedure idfp6
  Id f  =      r.u^-6;
  Id f1 =   -6*r.u^-7;
  Id f2 =   42*r.u^-8;
  Id f3 = -336*r.u^-9;
#endprocedure
#+