FORM 4.0 (Jun 30 2012) 64-bits                   Run: Mon Feb  4 11:47:30 2013
    *   Title:  POINT CHARGE:
    *           derivatives of retarded and advanced quantities
    *    File:  derivs.frm
    *   Tests:  derivs-test.frm
    *  Author:  David N. Williams
    * License:  Creative Commons Attribution-Share Alike
    * Started:  July 12, 2012
    * Revised:  July 12-16, 2012
    *           August 27-31, 2012
    *           September 1-3,5,8, 2012
    *
    * All parts of this program not in the public domain are:
    *
    *   Copyright (C) 1987-1990, 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.
    *
    *   Input:  none
    *  Output:  derivs.sav
    
    #if 0
    This program computes the partial derivatives of the retarded or
    advanced four-vectors r, u, a, and j described below, of their
    explicit dot products, and of the scalar function f = f(r⋅u). 
    These may be combined as tensor and scalar products in
    expressions.
    
    The output in derivs.sav is intended to be a comprehensive set
    of monomials with partial derivatives through second order and
    up to quadratic tensors in the vectors, with two exceptions:
    
    1. Only one ordering for tensor products is provided.  For
       example, dru(mu,nu,la) is provided, but dur(mu,nu,la) is not,
       because it is equal to dru(mu,la,mu).
    
    2. Partial derivatives that would produce second derivatives of
       the jerk, d²j/dτ², are omitted.
    
    Derivatives of third order are included on an as needed basis.
    
    All quantities are expressed in terms of the speed of light c
    and the following four-vectors:
    
      x       Space-time field point, x ∈ ℝ⁴.
      y(τ)    Particle four-position at proper time τ.  The
              curve in ℝ⁴ is the "world line".
      u(τ)    Particle four-velocity dy/dτ.
      a(τ)    Particle four-acceleration du/dτ. 
      j(τ)    Particle four-jerk da/dτ.
      r       Light-cone displacement r = x - y[τ(x)], where τ(x)
              is either of the two solutions of r⋅r = 0.
      ∂       Four-vector gradient with respect to x.
    
    Unless otherwise stated, it is always assumed that τ is
    evaluated at the retarded or the advanced solution of r⋅r = 0
    for τ(x).
    
    Our Lorentz metric is (+---), but is never explicit, except in
    the identification u⋅u = c².  It is implicit in the statement
    that the retarded/advanced values of sgn(r⋅u) = sgn(r⁰) are ±1,
    respectively.
    
    All derivatives computed are valid for both retarded and
    advanced evaluations of τ.
    
    NOTATION
    
      f,f1,f2  Function f(r⋅u) and its first and second derivatives
               f'(r⋅u) and f"(r⋅u).
      d        Four-vector gradient ∂.		
    
    Characteristic examples of the convention for the naming of
    results are dfru(mu,nu,la) for ∂(fru), where ∂, r, and u have
    indices mu, nu, and la respectively, and [d.fru](mu) for
    ∂⋅(fru), where the indices of ∂ and r are summed, and mu is the
    index of u.  The notation for d is an operator notation in that
    it acts on everything to the right.  Occasionally clarifying
    parentheses are provided when dot products are involved, as in
    [r.d(r.u)].
    
    USAGE
    
    Quantities that are to be differentiated by the differentiation
    procedure must be written with square bracket names from one of
    the FORM sets of cfunctions, cscal or cvec, declared below.
    Possibly a future version of this program will include a
    translator so that expressions without bracketed quantities can
    be used.
    
    As cfunctions, all dot products of bracketed vectors in the input
    expression must have explicit repeated indices, to be recognized
    for differentiation.  See the global expression for [r.d(r.u)]
    below for an example.
    
    The differentiator can be iterated on its unconverted output,
    which consists of bracketed cfunctions.  Two conversion
    procedures idvec is used to convert bracketed cfunctions in the
    final result to unbracketed FORM vectors or symbols.  There is
    also a procedure idmetric to account for the Lorentz metric. 
    See the calculations below for examples of usage.
    
    The differentiator does not work properly on powers of scalars,
    like [f]^3 or [r.u^-1]^2.  Instead, they have to be flattened,
    as in
    
      G [ddf^3(r.u)^-1)^2](mu,nu) = [f]*[f]*[f]*[r.u^-1]*[r.u^-1];
      #call partial(nu)
      #call partial(mu)
      #call idvec
      #call idmetric
    #endif
    
    *** DEBUGGING
    
    #procedure try
    print;
    .store
    save derivs.sav;
    .end
    #endprocedure
    
    *** DECLARATIONS
    
    #include derivs-ra.h
    #-
    
    off statistics;
    I mu,nu,la,ro;
    .global
    
    *** MODULES
    
    G dr(la,nu) = [r](nu);
    G du(la,nu) = [u](nu);
    G da(la,nu) = [a](nu);
    G dj(la,nu) = [j](nu);
    
    #call partial(la)
    #call idvec;
    print;
    .store

   dr(la,nu) =
       - r(la)*u(nu)*r.u^-1 + g(nu,la);

   du(la,nu) =
      r(la)*a(nu)*r.u^-1;

   da(la,nu) =
      r(la)*j(nu)*r.u^-1;

   dj(la,nu) =
      r(la)*jdot(nu)*r.u^-1;

    
    G [d(r.u)](la) = [r](mu)*[u](mu);
    G [d(r.a)](la) = [r](mu)*[a](mu);
    
    #call partial(la)
    #call idvec
    #call idmetric
    print;
    .store

   [d(r.u)](la) =
      r(la)*r.u^-1*r.a - r(la)*r.u^-1*c^2 + u(la);

   [d(r.a)](la) =
      r(la)*r.u^-1*r.j + a(la);

    
    G [d.r] = [r](mu);
    G [d.u] = [u](mu);
    G [d.a] = [a](mu);
    G [d.j] = [j](mu);
    
    #call partial(mu)
    #call idvec
    #call idmetric
    print;
    .store

   [d.r] =
      3;

   [d.u] =
      r.u^-1*r.a;

   [d.a] =
      r.u^-1*r.j;

   [d.j] =
      r.u^-1*r.jdot;

    
    G [r.dr](mu) = r(nu)*[r](mu); 
    G [r.du](mu) = r(nu)*[u](mu); 
    G [r.da](mu) = r(nu)*[a](mu); 
    G [r.dj](mu) = r(nu)*[j](mu); 
    
    G [u.dr](mu) = u(nu)*[r](mu); 
    G [u.du](mu) = u(nu)*[u](mu); 
    G [u.da](mu) = u(nu)*[a](mu); 
    
    G [a.dr](mu) = a(nu)*[r](mu); 
    G [a.du](mu) = a(nu)*[u](mu); 
    G [a.da](mu) = a(nu)*[a](mu); 
    
    #call partial(nu)
    #call idvec
    #call idmetric
    print;
    .store

   [r.dr](mu) =
      r(mu);

   [r.du](mu) = 0;

   [r.da](mu) = 0;

   [r.dj](mu) = 0;

   [u.dr](mu) = 0;

   [u.du](mu) =
      a(mu);

   [u.da](mu) =
      j(mu);

   [a.dr](mu) =
       - u(mu)*r.u^-1*r.a + a(mu);

   [a.du](mu) =
      a(mu)*r.u^-1*r.a;

   [a.da](mu) =
      j(mu)*r.u^-1*r.a;

    
    G df(la)  = [f];
    G df1(la) = [f1];
    #call partial(la)
    #call idvec
    #call idmetric
    print;
    .store

   df(la) =
      r(la)*r.u^-1*r.a*f1 - r(la)*r.u^-1*c^2*f1 + u(la)*f1;

   df1(la) =
      r(la)*r.u^-1*r.a*f2 - r(la)*r.u^-1*c^2*f2 + u(la)*f2;

    
    G [r.d(r.u)] = r(mu)*[r](nu)*[u](nu);
    G [r.d(r.a)] = r(mu)*[r](nu)*[a](nu);
    G [u.d(r.u)] = u(mu)*[r](nu)*[u](nu);
    G [u.d(r.a)] = u(mu)*[r](nu)*[a](nu);
    G [a.d(r.u)] = a(mu)*[r](nu)*[u](nu);
    G [a.d(r.a)] = a(mu)*[r](nu)*[a](nu);
    #call partial(la)
    multiply g(la,mu);
    #call idvec
    #call idmetric
    print;
    .store

   [r.d(r.u)] =
      r.u;

   [r.d(r.a)] =
      r.a;

   [u.d(r.u)] =
      r.a;

   [u.d(r.a)] =
      r.j;

   [a.d(r.u)] =
       - r.u^-1*r.a*c^2 + r.u^-1*r.a^2;

   [a.d(r.a)] =
      r.u^-1*r.a*r.j + a.a;

    
    G [r.df] = r(mu)*[f];
    G [u.df] = u(mu)*[f];
    #call partial(mu)
    #call idvec
    #call idmetric
    print;
    .store

   [r.df] =
      r.u*f1;

   [u.df] =
      r.a*f1;

    
    G dfr(la,nu)  = [f]*[r](nu);
    G dfu(la,nu)  = [f]*[u](nu);
    G dfa(la,nu)  = [f]*[a](nu);
    G dfj(la,nu)  = [f]*[j](nu);
    G df1r(la,nu) = [f1]*[r](nu);
    G df1u(la,nu) = [f1]*[u](nu);
    G df1a(la,nu) = [f1]*[a](nu);
    G df1j(la,nu) = [f1]*[j](nu);
    #call partial(la)
    #call idvec
    #call idmetric
    print;
    .store

   dfr(la,nu) =
      r(nu)*r(la)*r.u^-1*r.a*f1 - r(nu)*r(la)*r.u^-1*c^2*f1 + r(nu)*u(la)*f1
       - r(la)*u(nu)*r.u^-1*f + g(nu,la)*f;

   dfu(la,nu) =
      r(la)*u(nu)*r.u^-1*r.a*f1 - r(la)*u(nu)*r.u^-1*c^2*f1 + r(la)*a(nu)*
      r.u^-1*f + u(nu)*u(la)*f1;

   dfa(la,nu) =
      r(la)*a(nu)*r.u^-1*r.a*f1 - r(la)*a(nu)*r.u^-1*c^2*f1 + r(la)*j(nu)*
      r.u^-1*f + u(la)*a(nu)*f1;

   dfj(la,nu) =
      r(la)*j(nu)*r.u^-1*r.a*f1 - r(la)*j(nu)*r.u^-1*c^2*f1 + r(la)*jdot(nu)*
      r.u^-1*f + u(la)*j(nu)*f1;

   df1r(la,nu) =
      r(nu)*r(la)*r.u^-1*r.a*f2 - r(nu)*r(la)*r.u^-1*c^2*f2 + r(nu)*u(la)*f2
       - r(la)*u(nu)*r.u^-1*f1 + g(nu,la)*f1;

   df1u(la,nu) =
      r(la)*u(nu)*r.u^-1*r.a*f2 - r(la)*u(nu)*r.u^-1*c^2*f2 + r(la)*a(nu)*
      r.u^-1*f1 + u(nu)*u(la)*f2;

   df1a(la,nu) =
      r(la)*a(nu)*r.u^-1*r.a*f2 - r(la)*a(nu)*r.u^-1*c^2*f2 + r(la)*j(nu)*
      r.u^-1*f1 + u(la)*a(nu)*f2;

   df1j(la,nu) =
      r(la)*j(nu)*r.u^-1*r.a*f2 - r(la)*j(nu)*r.u^-1*c^2*f2 + r(la)*jdot(nu)*
      r.u^-1*f1 + u(la)*j(nu)*f2;

    
    G drr(la,mu,nu) = [r](mu)*[r](nu);
    G dru(la,mu,nu) = [r](mu)*[u](nu);
    G dra(la,mu,nu) = [r](mu)*[a](nu);
    G drj(la,mu,nu) = [r](mu)*[j](nu);
    G duu(la,mu,nu) = [u](mu)*[u](nu);
    G dua(la,mu,nu) = [u](mu)*[a](nu);
    G duj(la,mu,nu) = [u](mu)*[j](nu);
    G daa(la,mu,nu) = [a](mu)*[a](nu);
    G daj(la,mu,nu) = [a](mu)*[j](nu);
    
    G dfrr(la,mu,nu) = [f]*[r](mu)*[r](nu);
    G dfru(la,mu,nu) = [f]*[r](mu)*[u](nu);
    G dfra(la,mu,nu) = [f]*[r](mu)*[a](nu);
    G dfrj(la,mu,nu) = [f]*[r](mu)*[j](nu);
    G dfuu(la,mu,nu) = [f]*[u](mu)*[u](nu);
    G dfua(la,mu,nu) = [f]*[u](mu)*[a](nu);
    G dfuj(la,mu,nu) = [f]*[u](mu)*[j](nu);
    G dfaa(la,mu,nu) = [f]*[a](mu)*[a](nu);
    G dfaj(la,mu,nu) = [f]*[a](mu)*[j](nu);
    
    #call partial(la)
    #call idvec
    #call idmetric
    print;
    .store

   drr(la,mu,nu) =
       - r(mu)*r(la)*u(nu)*r.u^-1 - r(nu)*r(la)*u(mu)*r.u^-1 + g(mu,la)*r(nu)
       + g(nu,la)*r(mu);

   dru(la,mu,nu) =
      r(mu)*r(la)*a(nu)*r.u^-1 - r(la)*u(mu)*u(nu)*r.u^-1 + g(mu,la)*u(nu);

   dra(la,mu,nu) =
      r(mu)*r(la)*j(nu)*r.u^-1 - r(la)*u(mu)*a(nu)*r.u^-1 + g(mu,la)*a(nu);

   drj(la,mu,nu) =
      r(mu)*r(la)*jdot(nu)*r.u^-1 - r(la)*u(mu)*j(nu)*r.u^-1 + g(mu,la)*j(nu);

   duu(la,mu,nu) =
      r(la)*u(mu)*a(nu)*r.u^-1 + r(la)*u(nu)*a(mu)*r.u^-1;

   dua(la,mu,nu) =
      r(la)*u(mu)*j(nu)*r.u^-1 + r(la)*a(mu)*a(nu)*r.u^-1;

   duj(la,mu,nu) =
      r(la)*u(mu)*jdot(nu)*r.u^-1 + r(la)*a(mu)*j(nu)*r.u^-1;

   daa(la,mu,nu) =
      r(la)*a(mu)*j(nu)*r.u^-1 + r(la)*a(nu)*j(mu)*r.u^-1;

   daj(la,mu,nu) =
      r(la)*a(mu)*jdot(nu)*r.u^-1 + r(la)*j(mu)*j(nu)*r.u^-1;

   dfrr(la,mu,nu) =
      r(mu)*r(nu)*r(la)*r.u^-1*r.a*f1 - r(mu)*r(nu)*r(la)*r.u^-1*c^2*f1 + 
      r(mu)*r(nu)*u(la)*f1 - r(mu)*r(la)*u(nu)*r.u^-1*f - r(nu)*r(la)*u(mu)*
      r.u^-1*f + g(mu,la)*r(nu)*f + g(nu,la)*r(mu)*f;

   dfru(la,mu,nu) =
      r(mu)*r(la)*u(nu)*r.u^-1*r.a*f1 - r(mu)*r(la)*u(nu)*r.u^-1*c^2*f1 + 
      r(mu)*r(la)*a(nu)*r.u^-1*f + r(mu)*u(nu)*u(la)*f1 - r(la)*u(mu)*u(nu)*
      r.u^-1*f + g(mu,la)*u(nu)*f;

   dfra(la,mu,nu) =
      r(mu)*r(la)*a(nu)*r.u^-1*r.a*f1 - r(mu)*r(la)*a(nu)*r.u^-1*c^2*f1 + 
      r(mu)*r(la)*j(nu)*r.u^-1*f + r(mu)*u(la)*a(nu)*f1 - r(la)*u(mu)*a(nu)*
      r.u^-1*f + g(mu,la)*a(nu)*f;

   dfrj(la,mu,nu) =
      r(mu)*r(la)*j(nu)*r.u^-1*r.a*f1 - r(mu)*r(la)*j(nu)*r.u^-1*c^2*f1 + 
      r(mu)*r(la)*jdot(nu)*r.u^-1*f + r(mu)*u(la)*j(nu)*f1 - r(la)*u(mu)*j(nu)
      *r.u^-1*f + g(mu,la)*j(nu)*f;

   dfuu(la,mu,nu) =
      r(la)*u(mu)*u(nu)*r.u^-1*r.a*f1 - r(la)*u(mu)*u(nu)*r.u^-1*c^2*f1 + 
      r(la)*u(mu)*a(nu)*r.u^-1*f + r(la)*u(nu)*a(mu)*r.u^-1*f + u(mu)*u(nu)*
      u(la)*f1;

   dfua(la,mu,nu) =
      r(la)*u(mu)*a(nu)*r.u^-1*r.a*f1 - r(la)*u(mu)*a(nu)*r.u^-1*c^2*f1 + 
      r(la)*u(mu)*j(nu)*r.u^-1*f + r(la)*a(mu)*a(nu)*r.u^-1*f + u(mu)*u(la)*
      a(nu)*f1;

   dfuj(la,mu,nu) =
      r(la)*u(mu)*j(nu)*r.u^-1*r.a*f1 - r(la)*u(mu)*j(nu)*r.u^-1*c^2*f1 + 
      r(la)*u(mu)*jdot(nu)*r.u^-1*f + r(la)*a(mu)*j(nu)*r.u^-1*f + u(mu)*u(la)
      *j(nu)*f1;

   dfaa(la,mu,nu) =
      r(la)*a(mu)*a(nu)*r.u^-1*r.a*f1 - r(la)*a(mu)*a(nu)*r.u^-1*c^2*f1 + 
      r(la)*a(mu)*j(nu)*r.u^-1*f + r(la)*a(nu)*j(mu)*r.u^-1*f + u(la)*a(mu)*
      a(nu)*f1;

   dfaj(la,mu,nu) =
      r(la)*a(mu)*j(nu)*r.u^-1*r.a*f1 - r(la)*a(mu)*j(nu)*r.u^-1*c^2*f1 + 
      r(la)*a(mu)*jdot(nu)*r.u^-1*f + r(la)*j(mu)*j(nu)*r.u^-1*f + u(la)*a(mu)
      *j(nu)*f1;

    
    G [d.fr] = dfr(mu,mu);
    G [d.fu] = dfu(mu,mu);
    G [d.fa] = dfa(mu,mu);
    G [d.fj] = dfj(mu,mu);
    print;
    .store

   [d.fr] =
       - f + r.r*r.u^-1*r.a*f1 - r.r*r.u^-1*c^2*f1 + r.u*f1 + g(mu,mu)*f;

   [d.fu] =
       - c^2*f1 + r.u^-1*r.a*f + r.a*f1 + u.u*f1;

   [d.fa] =
       - r.u^-1*r.a*c^2*f1 + r.u^-1*r.a^2*f1 + r.u^-1*r.j*f + u.a*f1;

   [d.fj] =
      r.u^-1*r.a*r.j*f1 - r.u^-1*r.j*c^2*f1 + r.u^-1*r.jdot*f + u.j*f1;

    
    G [d.frr](mu) = dfrr(nu,nu,mu); 
    G [d.fru](mu) = dfru(nu,nu,mu); 
    G [d.fra](mu) = dfra(nu,nu,mu); 
    G [d.fur](mu) = dfru(nu,mu,nu); 
    G [d.fuu](mu) = dfuu(nu,nu,mu); 
    G [d.fua](mu) = dfua(nu,nu,mu); 
    G [d.far](mu) = dfra(nu,mu,nu); 
    G [d.fau](mu) = dfua(nu,mu,nu); 
    G [d.faa](mu) = dfaa(nu,nu,mu); 
    #call idmetric
    print;
    .store

   [d.frr](mu) =
      4*r(mu)*f + r(mu)*r.u*f1;

   [d.fru](mu) =
      3*u(mu)*f + u(mu)*r.u*f1;

   [d.fra](mu) =
      3*a(mu)*f + a(mu)*r.u*f1;

   [d.fur](mu) =
      r(mu)*r.u^-1*r.a*f + r(mu)*r.a*f1;

   [d.fuu](mu) =
      u(mu)*r.u^-1*r.a*f + u(mu)*r.a*f1 + a(mu)*f;

   [d.fua](mu) =
      a(mu)*r.u^-1*r.a*f + a(mu)*r.a*f1 + j(mu)*f;

   [d.far](mu) =
       - r(mu)*r.u^-1*r.a*c^2*f1 + r(mu)*r.u^-1*r.a^2*f1 + r(mu)*r.u^-1*r.j*f
       - u(mu)*r.u^-1*r.a*f + a(mu)*f;

   [d.fau](mu) =
       - u(mu)*r.u^-1*r.a*c^2*f1 + u(mu)*r.u^-1*r.a^2*f1 + u(mu)*r.u^-1*r.j*f
       + a(mu)*r.u^-1*r.a*f;

   [d.faa](mu) =
       - a(mu)*r.u^-1*r.a*c^2*f1 + a(mu)*r.u^-1*r.a^2*f1 + a(mu)*r.u^-1*r.j*f
       + j(mu)*r.u^-1*r.a*f;

    
    G ddf(mu,nu) = [f];
    
    G ddr(mu,nu,la) = [r](la);
    G ddu(mu,nu,la) = [u](la);
    G dda(mu,nu,la) = [a](la);
    
    G ddfr(mu,nu,la) = [f]*[r](la);
    G ddfu(mu,nu,la) = [f]*[u](la);
    G ddfa(mu,nu,la) = [f]*[a](la);
    
    G ddrr(mu,nu,la,ro) = [r](la)*[r](ro);
    G ddru(mu,nu,la,ro) = [r](la)*[u](ro);
    G ddra(mu,nu,la,ro) = [r](la)*[a](ro);
    G dduu(mu,nu,la,ro) = [u](la)*[u](ro);
    G ddua(mu,nu,la,ro) = [u](la)*[a](ro);
    G ddaa(mu,nu,la,ro) = [a](la)*[a](ro);
    
    G ddfrr(mu,nu,la,ro) = [f]*[r](la)*[r](ro);
    G ddfru(mu,nu,la,ro) = [f]*[r](la)*[u](ro);
    G ddfra(mu,nu,la,ro) = [f]*[r](la)*[a](ro);
    G ddfuu(mu,nu,la,ro) = [f]*[u](la)*[u](ro);
    G ddfua(mu,nu,la,ro) = [f]*[u](la)*[a](ro);
    G ddfaa(mu,nu,la,ro) = [f]*[a](la)*[a](ro);
    
    #call partial(nu)
    #call partial(mu)
    #call idvec
    #call idmetric
    print;
    .store

   ddf(mu,nu) =
      2*r(mu)*r(nu)*r.u^-3*r.a*c^2*f1 - r(mu)*r(nu)*r.u^-3*r.a^2*f1 - r(mu)*
      r(nu)*r.u^-3*c^4*f1 - 2*r(mu)*r(nu)*r.u^-2*r.a*c^2*f2 + r(mu)*r(nu)*
      r.u^-2*r.a^2*f2 + r(mu)*r(nu)*r.u^-2*r.j*f1 + r(mu)*r(nu)*r.u^-2*c^4*f2
       - r(mu)*u(nu)*r.u^-2*r.a*f1 + r(mu)*u(nu)*r.u^-2*c^2*f1 + r(mu)*u(nu)*
      r.u^-1*r.a*f2 - r(mu)*u(nu)*r.u^-1*c^2*f2 + r(mu)*a(nu)*r.u^-1*f1 - 
      r(nu)*u(mu)*r.u^-2*r.a*f1 + r(nu)*u(mu)*r.u^-2*c^2*f1 + r(nu)*u(mu)*
      r.u^-1*r.a*f2 - r(nu)*u(mu)*r.u^-1*c^2*f2 + r(nu)*a(mu)*r.u^-1*f1 + 
      u(mu)*u(nu)*f2 + g(mu,nu)*r.u^-1*r.a*f1 - g(mu,nu)*r.u^-1*c^2*f1;

   ddr(mu,nu,la) =
      r(mu)*r(nu)*u(la)*r.u^-3*r.a - r(mu)*r(nu)*u(la)*r.u^-3*c^2 - r(mu)*
      r(nu)*a(la)*r.u^-2 + r(mu)*u(nu)*u(la)*r.u^-2 + r(nu)*u(mu)*u(la)*r.u^-2
       - g(mu,nu)*u(la)*r.u^-1;

   ddu(mu,nu,la) =
       - r(mu)*r(nu)*a(la)*r.u^-3*r.a + r(mu)*r(nu)*a(la)*r.u^-3*c^2 + r(mu)*
      r(nu)*j(la)*r.u^-2 - r(mu)*u(nu)*a(la)*r.u^-2 - r(nu)*u(mu)*a(la)*r.u^-2
       + g(mu,nu)*a(la)*r.u^-1;

   dda(mu,nu,la) =
       - r(mu)*r(nu)*j(la)*r.u^-3*r.a + r(mu)*r(nu)*j(la)*r.u^-3*c^2 + r(mu)*
      r(nu)*jdot(la)*r.u^-2 - r(mu)*u(nu)*j(la)*r.u^-2 - r(nu)*u(mu)*j(la)*
      r.u^-2 + g(mu,nu)*j(la)*r.u^-1;

   ddfr(mu,nu,la) =
      2*r(mu)*r(nu)*r(la)*r.u^-3*r.a*c^2*f1 - r(mu)*r(nu)*r(la)*r.u^-3*r.a^2*
      f1 - r(mu)*r(nu)*r(la)*r.u^-3*c^4*f1 - 2*r(mu)*r(nu)*r(la)*r.u^-2*r.a*
      c^2*f2 + r(mu)*r(nu)*r(la)*r.u^-2*r.a^2*f2 + r(mu)*r(nu)*r(la)*r.u^-2*
      r.j*f1 + r(mu)*r(nu)*r(la)*r.u^-2*c^4*f2 + r(mu)*r(nu)*u(la)*r.u^-3*r.a*
      f - r(mu)*r(nu)*u(la)*r.u^-3*c^2*f - 2*r(mu)*r(nu)*u(la)*r.u^-2*r.a*f1
       + 2*r(mu)*r(nu)*u(la)*r.u^-2*c^2*f1 - r(mu)*r(nu)*a(la)*r.u^-2*f - 
      r(mu)*r(la)*u(nu)*r.u^-2*r.a*f1 + r(mu)*r(la)*u(nu)*r.u^-2*c^2*f1 + 
      r(mu)*r(la)*u(nu)*r.u^-1*r.a*f2 - r(mu)*r(la)*u(nu)*r.u^-1*c^2*f2 + 
      r(mu)*r(la)*a(nu)*r.u^-1*f1 + r(mu)*u(nu)*u(la)*r.u^-2*f - r(mu)*u(nu)*
      u(la)*r.u^-1*f1 - r(nu)*r(la)*u(mu)*r.u^-2*r.a*f1 + r(nu)*r(la)*u(mu)*
      r.u^-2*c^2*f1 + r(nu)*r(la)*u(mu)*r.u^-1*r.a*f2 - r(nu)*r(la)*u(mu)*
      r.u^-1*c^2*f2 + r(nu)*r(la)*a(mu)*r.u^-1*f1 + r(nu)*u(mu)*u(la)*r.u^-2*f
       - r(nu)*u(mu)*u(la)*r.u^-1*f1 + r(la)*u(mu)*u(nu)*f2 + g(mu,nu)*r(la)*
      r.u^-1*r.a*f1 - g(mu,nu)*r(la)*r.u^-1*c^2*f1 - g(mu,nu)*u(la)*r.u^-1*f
       + g(mu,la)*r(nu)*r.u^-1*r.a*f1 - g(mu,la)*r(nu)*r.u^-1*c^2*f1 + g(mu,la
      )*u(nu)*f1 + g(nu,la)*r(mu)*r.u^-1*r.a*f1 - g(nu,la)*r(mu)*r.u^-1*c^2*f1
       + g(nu,la)*u(mu)*f1;

   ddfu(mu,nu,la) =
      2*r(mu)*r(nu)*u(la)*r.u^-3*r.a*c^2*f1 - r(mu)*r(nu)*u(la)*r.u^-3*r.a^2*
      f1 - r(mu)*r(nu)*u(la)*r.u^-3*c^4*f1 - 2*r(mu)*r(nu)*u(la)*r.u^-2*r.a*
      c^2*f2 + r(mu)*r(nu)*u(la)*r.u^-2*r.a^2*f2 + r(mu)*r(nu)*u(la)*r.u^-2*
      r.j*f1 + r(mu)*r(nu)*u(la)*r.u^-2*c^4*f2 - r(mu)*r(nu)*a(la)*r.u^-3*r.a*
      f + r(mu)*r(nu)*a(la)*r.u^-3*c^2*f + 2*r(mu)*r(nu)*a(la)*r.u^-2*r.a*f1
       - 2*r(mu)*r(nu)*a(la)*r.u^-2*c^2*f1 + r(mu)*r(nu)*j(la)*r.u^-2*f - 
      r(mu)*u(nu)*u(la)*r.u^-2*r.a*f1 + r(mu)*u(nu)*u(la)*r.u^-2*c^2*f1 + 
      r(mu)*u(nu)*u(la)*r.u^-1*r.a*f2 - r(mu)*u(nu)*u(la)*r.u^-1*c^2*f2 - 
      r(mu)*u(nu)*a(la)*r.u^-2*f + r(mu)*u(nu)*a(la)*r.u^-1*f1 + r(mu)*u(la)*
      a(nu)*r.u^-1*f1 - r(nu)*u(mu)*u(la)*r.u^-2*r.a*f1 + r(nu)*u(mu)*u(la)*
      r.u^-2*c^2*f1 + r(nu)*u(mu)*u(la)*r.u^-1*r.a*f2 - r(nu)*u(mu)*u(la)*
      r.u^-1*c^2*f2 - r(nu)*u(mu)*a(la)*r.u^-2*f + r(nu)*u(mu)*a(la)*r.u^-1*f1
       + r(nu)*u(la)*a(mu)*r.u^-1*f1 + u(mu)*u(nu)*u(la)*f2 + g(mu,nu)*u(la)*
      r.u^-1*r.a*f1 - g(mu,nu)*u(la)*r.u^-1*c^2*f1 + g(mu,nu)*a(la)*r.u^-1*f;

   ddfa(mu,nu,la) =
      2*r(mu)*r(nu)*a(la)*r.u^-3*r.a*c^2*f1 - r(mu)*r(nu)*a(la)*r.u^-3*r.a^2*
      f1 - r(mu)*r(nu)*a(la)*r.u^-3*c^4*f1 - 2*r(mu)*r(nu)*a(la)*r.u^-2*r.a*
      c^2*f2 + r(mu)*r(nu)*a(la)*r.u^-2*r.a^2*f2 + r(mu)*r(nu)*a(la)*r.u^-2*
      r.j*f1 + r(mu)*r(nu)*a(la)*r.u^-2*c^4*f2 - r(mu)*r(nu)*j(la)*r.u^-3*r.a*
      f + r(mu)*r(nu)*j(la)*r.u^-3*c^2*f + 2*r(mu)*r(nu)*j(la)*r.u^-2*r.a*f1
       - 2*r(mu)*r(nu)*j(la)*r.u^-2*c^2*f1 + r(mu)*r(nu)*jdot(la)*r.u^-2*f - 
      r(mu)*u(nu)*a(la)*r.u^-2*r.a*f1 + r(mu)*u(nu)*a(la)*r.u^-2*c^2*f1 + 
      r(mu)*u(nu)*a(la)*r.u^-1*r.a*f2 - r(mu)*u(nu)*a(la)*r.u^-1*c^2*f2 - 
      r(mu)*u(nu)*j(la)*r.u^-2*f + r(mu)*u(nu)*j(la)*r.u^-1*f1 + r(mu)*a(nu)*
      a(la)*r.u^-1*f1 - r(nu)*u(mu)*a(la)*r.u^-2*r.a*f1 + r(nu)*u(mu)*a(la)*
      r.u^-2*c^2*f1 + r(nu)*u(mu)*a(la)*r.u^-1*r.a*f2 - r(nu)*u(mu)*a(la)*
      r.u^-1*c^2*f2 - r(nu)*u(mu)*j(la)*r.u^-2*f + r(nu)*u(mu)*j(la)*r.u^-1*f1
       + r(nu)*a(mu)*a(la)*r.u^-1*f1 + u(mu)*u(nu)*a(la)*f2 + g(mu,nu)*a(la)*
      r.u^-1*r.a*f1 - g(mu,nu)*a(la)*r.u^-1*c^2*f1 + g(mu,nu)*j(la)*r.u^-1*f;

   ddrr(mu,nu,la,ro) =
      r(mu)*r(nu)*r(la)*u(ro)*r.u^-3*r.a - r(mu)*r(nu)*r(la)*u(ro)*r.u^-3*c^2
       - r(mu)*r(nu)*r(la)*a(ro)*r.u^-2 + r(mu)*r(nu)*r(ro)*u(la)*r.u^-3*r.a
       - r(mu)*r(nu)*r(ro)*u(la)*r.u^-3*c^2 - r(mu)*r(nu)*r(ro)*a(la)*r.u^-2
       + 2*r(mu)*r(nu)*u(la)*u(ro)*r.u^-2 + r(mu)*r(la)*u(nu)*u(ro)*r.u^-2 + 
      r(mu)*r(ro)*u(nu)*u(la)*r.u^-2 + r(nu)*r(la)*u(mu)*u(ro)*r.u^-2 + r(nu)*
      r(ro)*u(mu)*u(la)*r.u^-2 - g(mu,nu)*r(la)*u(ro)*r.u^-1 - g(mu,nu)*r(ro)*
      u(la)*r.u^-1 - g(mu,la)*r(nu)*u(ro)*r.u^-1 + g(mu,la)*g(nu,ro) - g(mu,ro
      )*r(nu)*u(la)*r.u^-1 + g(mu,ro)*g(nu,la) - g(nu,la)*r(mu)*u(ro)*r.u^-1
       - g(nu,ro)*r(mu)*u(la)*r.u^-1;

   ddru(mu,nu,la,ro) =
       - r(mu)*r(nu)*r(la)*a(ro)*r.u^-3*r.a + r(mu)*r(nu)*r(la)*a(ro)*r.u^-3*
      c^2 + r(mu)*r(nu)*r(la)*j(ro)*r.u^-2 + r(mu)*r(nu)*u(la)*u(ro)*r.u^-3*
      r.a - r(mu)*r(nu)*u(la)*u(ro)*r.u^-3*c^2 - 2*r(mu)*r(nu)*u(la)*a(ro)*
      r.u^-2 - r(mu)*r(nu)*u(ro)*a(la)*r.u^-2 - r(mu)*r(la)*u(nu)*a(ro)*r.u^-2
       + r(mu)*u(nu)*u(la)*u(ro)*r.u^-2 - r(nu)*r(la)*u(mu)*a(ro)*r.u^-2 + 
      r(nu)*u(mu)*u(la)*u(ro)*r.u^-2 + g(mu,nu)*r(la)*a(ro)*r.u^-1 - g(mu,nu)*
      u(la)*u(ro)*r.u^-1 + g(mu,la)*r(nu)*a(ro)*r.u^-1 + g(nu,la)*r(mu)*a(ro)*
      r.u^-1;

   ddra(mu,nu,la,ro) =
       - r(mu)*r(nu)*r(la)*j(ro)*r.u^-3*r.a + r(mu)*r(nu)*r(la)*j(ro)*r.u^-3*
      c^2 + r(mu)*r(nu)*r(la)*jdot(ro)*r.u^-2 + r(mu)*r(nu)*u(la)*a(ro)*r.u^-3
      *r.a - r(mu)*r(nu)*u(la)*a(ro)*r.u^-3*c^2 - 2*r(mu)*r(nu)*u(la)*j(ro)*
      r.u^-2 - r(mu)*r(nu)*a(la)*a(ro)*r.u^-2 - r(mu)*r(la)*u(nu)*j(ro)*r.u^-2
       + r(mu)*u(nu)*u(la)*a(ro)*r.u^-2 - r(nu)*r(la)*u(mu)*j(ro)*r.u^-2 + 
      r(nu)*u(mu)*u(la)*a(ro)*r.u^-2 + g(mu,nu)*r(la)*j(ro)*r.u^-1 - g(mu,nu)*
      u(la)*a(ro)*r.u^-1 + g(mu,la)*r(nu)*j(ro)*r.u^-1 + g(nu,la)*r(mu)*j(ro)*
      r.u^-1;

   dduu(mu,nu,la,ro) =
       - r(mu)*r(nu)*u(la)*a(ro)*r.u^-3*r.a + r(mu)*r(nu)*u(la)*a(ro)*r.u^-3*
      c^2 + r(mu)*r(nu)*u(la)*j(ro)*r.u^-2 - r(mu)*r(nu)*u(ro)*a(la)*r.u^-3*
      r.a + r(mu)*r(nu)*u(ro)*a(la)*r.u^-3*c^2 + r(mu)*r(nu)*u(ro)*j(la)*
      r.u^-2 + 2*r(mu)*r(nu)*a(la)*a(ro)*r.u^-2 - r(mu)*u(nu)*u(la)*a(ro)*
      r.u^-2 - r(mu)*u(nu)*u(ro)*a(la)*r.u^-2 - r(nu)*u(mu)*u(la)*a(ro)*r.u^-2
       - r(nu)*u(mu)*u(ro)*a(la)*r.u^-2 + g(mu,nu)*u(la)*a(ro)*r.u^-1 + g(mu,
      nu)*u(ro)*a(la)*r.u^-1;

   ddua(mu,nu,la,ro) =
       - r(mu)*r(nu)*u(la)*j(ro)*r.u^-3*r.a + r(mu)*r(nu)*u(la)*j(ro)*r.u^-3*
      c^2 + r(mu)*r(nu)*u(la)*jdot(ro)*r.u^-2 - r(mu)*r(nu)*a(la)*a(ro)*r.u^-3
      *r.a + r(mu)*r(nu)*a(la)*a(ro)*r.u^-3*c^2 + 2*r(mu)*r(nu)*a(la)*j(ro)*
      r.u^-2 + r(mu)*r(nu)*a(ro)*j(la)*r.u^-2 - r(mu)*u(nu)*u(la)*j(ro)*r.u^-2
       - r(mu)*u(nu)*a(la)*a(ro)*r.u^-2 - r(nu)*u(mu)*u(la)*j(ro)*r.u^-2 - 
      r(nu)*u(mu)*a(la)*a(ro)*r.u^-2 + g(mu,nu)*u(la)*j(ro)*r.u^-1 + g(mu,nu)*
      a(la)*a(ro)*r.u^-1;

   ddaa(mu,nu,la,ro) =
       - r(mu)*r(nu)*a(la)*j(ro)*r.u^-3*r.a + r(mu)*r(nu)*a(la)*j(ro)*r.u^-3*
      c^2 + r(mu)*r(nu)*a(la)*jdot(ro)*r.u^-2 - r(mu)*r(nu)*a(ro)*j(la)*r.u^-3
      *r.a + r(mu)*r(nu)*a(ro)*j(la)*r.u^-3*c^2 + r(mu)*r(nu)*a(ro)*jdot(la)*
      r.u^-2 + 2*r(mu)*r(nu)*j(la)*j(ro)*r.u^-2 - r(mu)*u(nu)*a(la)*j(ro)*
      r.u^-2 - r(mu)*u(nu)*a(ro)*j(la)*r.u^-2 - r(nu)*u(mu)*a(la)*j(ro)*r.u^-2
       - r(nu)*u(mu)*a(ro)*j(la)*r.u^-2 + g(mu,nu)*a(la)*j(ro)*r.u^-1 + g(mu,
      nu)*a(ro)*j(la)*r.u^-1;

   ddfrr(mu,nu,la,ro) =
      2*r(mu)*r(nu)*r(la)*r(ro)*r.u^-3*r.a*c^2*f1 - r(mu)*r(nu)*r(la)*r(ro)*
      r.u^-3*r.a^2*f1 - r(mu)*r(nu)*r(la)*r(ro)*r.u^-3*c^4*f1 - 2*r(mu)*r(nu)*
      r(la)*r(ro)*r.u^-2*r.a*c^2*f2 + r(mu)*r(nu)*r(la)*r(ro)*r.u^-2*r.a^2*f2
       + r(mu)*r(nu)*r(la)*r(ro)*r.u^-2*r.j*f1 + r(mu)*r(nu)*r(la)*r(ro)*
      r.u^-2*c^4*f2 + r(mu)*r(nu)*r(la)*u(ro)*r.u^-3*r.a*f - r(mu)*r(nu)*r(la)
      *u(ro)*r.u^-3*c^2*f - 2*r(mu)*r(nu)*r(la)*u(ro)*r.u^-2*r.a*f1 + 2*r(mu)*
      r(nu)*r(la)*u(ro)*r.u^-2*c^2*f1 - r(mu)*r(nu)*r(la)*a(ro)*r.u^-2*f + 
      r(mu)*r(nu)*r(ro)*u(la)*r.u^-3*r.a*f - r(mu)*r(nu)*r(ro)*u(la)*r.u^-3*
      c^2*f - 2*r(mu)*r(nu)*r(ro)*u(la)*r.u^-2*r.a*f1 + 2*r(mu)*r(nu)*r(ro)*
      u(la)*r.u^-2*c^2*f1 - r(mu)*r(nu)*r(ro)*a(la)*r.u^-2*f + 2*r(mu)*r(nu)*
      u(la)*u(ro)*r.u^-2*f - r(mu)*r(la)*r(ro)*u(nu)*r.u^-2*r.a*f1 + r(mu)*
      r(la)*r(ro)*u(nu)*r.u^-2*c^2*f1 + r(mu)*r(la)*r(ro)*u(nu)*r.u^-1*r.a*f2
       - r(mu)*r(la)*r(ro)*u(nu)*r.u^-1*c^2*f2 + r(mu)*r(la)*r(ro)*a(nu)*
      r.u^-1*f1 + r(mu)*r(la)*u(nu)*u(ro)*r.u^-2*f - r(mu)*r(la)*u(nu)*u(ro)*
      r.u^-1*f1 + r(mu)*r(ro)*u(nu)*u(la)*r.u^-2*f - r(mu)*r(ro)*u(nu)*u(la)*
      r.u^-1*f1 - r(nu)*r(la)*r(ro)*u(mu)*r.u^-2*r.a*f1 + r(nu)*r(la)*r(ro)*
      u(mu)*r.u^-2*c^2*f1 + r(nu)*r(la)*r(ro)*u(mu)*r.u^-1*r.a*f2 - r(nu)*
      r(la)*r(ro)*u(mu)*r.u^-1*c^2*f2 + r(nu)*r(la)*r(ro)*a(mu)*r.u^-1*f1 + 
      r(nu)*r(la)*u(mu)*u(ro)*r.u^-2*f - r(nu)*r(la)*u(mu)*u(ro)*r.u^-1*f1 + 
      r(nu)*r(ro)*u(mu)*u(la)*r.u^-2*f - r(nu)*r(ro)*u(mu)*u(la)*r.u^-1*f1 + 
      r(la)*r(ro)*u(mu)*u(nu)*f2 + g(mu,nu)*r(la)*r(ro)*r.u^-1*r.a*f1 - g(mu,
      nu)*r(la)*r(ro)*r.u^-1*c^2*f1 - g(mu,nu)*r(la)*u(ro)*r.u^-1*f - g(mu,nu)
      *r(ro)*u(la)*r.u^-1*f + g(mu,la)*r(nu)*r(ro)*r.u^-1*r.a*f1 - g(mu,la)*
      r(nu)*r(ro)*r.u^-1*c^2*f1 - g(mu,la)*r(nu)*u(ro)*r.u^-1*f + g(mu,la)*
      r(ro)*u(nu)*f1 + g(mu,la)*g(nu,ro)*f + g(mu,ro)*r(nu)*r(la)*r.u^-1*r.a*
      f1 - g(mu,ro)*r(nu)*r(la)*r.u^-1*c^2*f1 - g(mu,ro)*r(nu)*u(la)*r.u^-1*f
       + g(mu,ro)*r(la)*u(nu)*f1 + g(mu,ro)*g(nu,la)*f + g(nu,la)*r(mu)*r(ro)*
      r.u^-1*r.a*f1 - g(nu,la)*r(mu)*r(ro)*r.u^-1*c^2*f1 - g(nu,la)*r(mu)*
      u(ro)*r.u^-1*f + g(nu,la)*r(ro)*u(mu)*f1 + g(nu,ro)*r(mu)*r(la)*r.u^-1*
      r.a*f1 - g(nu,ro)*r(mu)*r(la)*r.u^-1*c^2*f1 - g(nu,ro)*r(mu)*u(la)*
      r.u^-1*f + g(nu,ro)*r(la)*u(mu)*f1;

   ddfru(mu,nu,la,ro) =
      2*r(mu)*r(nu)*r(la)*u(ro)*r.u^-3*r.a*c^2*f1 - r(mu)*r(nu)*r(la)*u(ro)*
      r.u^-3*r.a^2*f1 - r(mu)*r(nu)*r(la)*u(ro)*r.u^-3*c^4*f1 - 2*r(mu)*r(nu)*
      r(la)*u(ro)*r.u^-2*r.a*c^2*f2 + r(mu)*r(nu)*r(la)*u(ro)*r.u^-2*r.a^2*f2
       + r(mu)*r(nu)*r(la)*u(ro)*r.u^-2*r.j*f1 + r(mu)*r(nu)*r(la)*u(ro)*
      r.u^-2*c^4*f2 - r(mu)*r(nu)*r(la)*a(ro)*r.u^-3*r.a*f + r(mu)*r(nu)*r(la)
      *a(ro)*r.u^-3*c^2*f + 2*r(mu)*r(nu)*r(la)*a(ro)*r.u^-2*r.a*f1 - 2*r(mu)*
      r(nu)*r(la)*a(ro)*r.u^-2*c^2*f1 + r(mu)*r(nu)*r(la)*j(ro)*r.u^-2*f + 
      r(mu)*r(nu)*u(la)*u(ro)*r.u^-3*r.a*f - r(mu)*r(nu)*u(la)*u(ro)*r.u^-3*
      c^2*f - 2*r(mu)*r(nu)*u(la)*u(ro)*r.u^-2*r.a*f1 + 2*r(mu)*r(nu)*u(la)*
      u(ro)*r.u^-2*c^2*f1 - 2*r(mu)*r(nu)*u(la)*a(ro)*r.u^-2*f - r(mu)*r(nu)*
      u(ro)*a(la)*r.u^-2*f - r(mu)*r(la)*u(nu)*u(ro)*r.u^-2*r.a*f1 + r(mu)*
      r(la)*u(nu)*u(ro)*r.u^-2*c^2*f1 + r(mu)*r(la)*u(nu)*u(ro)*r.u^-1*r.a*f2
       - r(mu)*r(la)*u(nu)*u(ro)*r.u^-1*c^2*f2 - r(mu)*r(la)*u(nu)*a(ro)*
      r.u^-2*f + r(mu)*r(la)*u(nu)*a(ro)*r.u^-1*f1 + r(mu)*r(la)*u(ro)*a(nu)*
      r.u^-1*f1 + r(mu)*u(nu)*u(la)*u(ro)*r.u^-2*f - r(mu)*u(nu)*u(la)*u(ro)*
      r.u^-1*f1 - r(nu)*r(la)*u(mu)*u(ro)*r.u^-2*r.a*f1 + r(nu)*r(la)*u(mu)*
      u(ro)*r.u^-2*c^2*f1 + r(nu)*r(la)*u(mu)*u(ro)*r.u^-1*r.a*f2 - r(nu)*
      r(la)*u(mu)*u(ro)*r.u^-1*c^2*f2 - r(nu)*r(la)*u(mu)*a(ro)*r.u^-2*f + 
      r(nu)*r(la)*u(mu)*a(ro)*r.u^-1*f1 + r(nu)*r(la)*u(ro)*a(mu)*r.u^-1*f1 + 
      r(nu)*u(mu)*u(la)*u(ro)*r.u^-2*f - r(nu)*u(mu)*u(la)*u(ro)*r.u^-1*f1 + 
      r(la)*u(mu)*u(nu)*u(ro)*f2 + g(mu,nu)*r(la)*u(ro)*r.u^-1*r.a*f1 - g(mu,
      nu)*r(la)*u(ro)*r.u^-1*c^2*f1 + g(mu,nu)*r(la)*a(ro)*r.u^-1*f - g(mu,nu)
      *u(la)*u(ro)*r.u^-1*f + g(mu,la)*r(nu)*u(ro)*r.u^-1*r.a*f1 - g(mu,la)*
      r(nu)*u(ro)*r.u^-1*c^2*f1 + g(mu,la)*r(nu)*a(ro)*r.u^-1*f + g(mu,la)*
      u(nu)*u(ro)*f1 + g(nu,la)*r(mu)*u(ro)*r.u^-1*r.a*f1 - g(nu,la)*r(mu)*
      u(ro)*r.u^-1*c^2*f1 + g(nu,la)*r(mu)*a(ro)*r.u^-1*f + g(nu,la)*u(mu)*
      u(ro)*f1;

   ddfra(mu,nu,la,ro) =
      2*r(mu)*r(nu)*r(la)*a(ro)*r.u^-3*r.a*c^2*f1 - r(mu)*r(nu)*r(la)*a(ro)*
      r.u^-3*r.a^2*f1 - r(mu)*r(nu)*r(la)*a(ro)*r.u^-3*c^4*f1 - 2*r(mu)*r(nu)*
      r(la)*a(ro)*r.u^-2*r.a*c^2*f2 + r(mu)*r(nu)*r(la)*a(ro)*r.u^-2*r.a^2*f2
       + r(mu)*r(nu)*r(la)*a(ro)*r.u^-2*r.j*f1 + r(mu)*r(nu)*r(la)*a(ro)*
      r.u^-2*c^4*f2 - r(mu)*r(nu)*r(la)*j(ro)*r.u^-3*r.a*f + r(mu)*r(nu)*r(la)
      *j(ro)*r.u^-3*c^2*f + 2*r(mu)*r(nu)*r(la)*j(ro)*r.u^-2*r.a*f1 - 2*r(mu)*
      r(nu)*r(la)*j(ro)*r.u^-2*c^2*f1 + r(mu)*r(nu)*r(la)*jdot(ro)*r.u^-2*f + 
      r(mu)*r(nu)*u(la)*a(ro)*r.u^-3*r.a*f - r(mu)*r(nu)*u(la)*a(ro)*r.u^-3*
      c^2*f - 2*r(mu)*r(nu)*u(la)*a(ro)*r.u^-2*r.a*f1 + 2*r(mu)*r(nu)*u(la)*
      a(ro)*r.u^-2*c^2*f1 - 2*r(mu)*r(nu)*u(la)*j(ro)*r.u^-2*f - r(mu)*r(nu)*
      a(la)*a(ro)*r.u^-2*f - r(mu)*r(la)*u(nu)*a(ro)*r.u^-2*r.a*f1 + r(mu)*
      r(la)*u(nu)*a(ro)*r.u^-2*c^2*f1 + r(mu)*r(la)*u(nu)*a(ro)*r.u^-1*r.a*f2
       - r(mu)*r(la)*u(nu)*a(ro)*r.u^-1*c^2*f2 - r(mu)*r(la)*u(nu)*j(ro)*
      r.u^-2*f + r(mu)*r(la)*u(nu)*j(ro)*r.u^-1*f1 + r(mu)*r(la)*a(nu)*a(ro)*
      r.u^-1*f1 + r(mu)*u(nu)*u(la)*a(ro)*r.u^-2*f - r(mu)*u(nu)*u(la)*a(ro)*
      r.u^-1*f1 - r(nu)*r(la)*u(mu)*a(ro)*r.u^-2*r.a*f1 + r(nu)*r(la)*u(mu)*
      a(ro)*r.u^-2*c^2*f1 + r(nu)*r(la)*u(mu)*a(ro)*r.u^-1*r.a*f2 - r(nu)*
      r(la)*u(mu)*a(ro)*r.u^-1*c^2*f2 - r(nu)*r(la)*u(mu)*j(ro)*r.u^-2*f + 
      r(nu)*r(la)*u(mu)*j(ro)*r.u^-1*f1 + r(nu)*r(la)*a(mu)*a(ro)*r.u^-1*f1 + 
      r(nu)*u(mu)*u(la)*a(ro)*r.u^-2*f - r(nu)*u(mu)*u(la)*a(ro)*r.u^-1*f1 + 
      r(la)*u(mu)*u(nu)*a(ro)*f2 + g(mu,nu)*r(la)*a(ro)*r.u^-1*r.a*f1 - g(mu,
      nu)*r(la)*a(ro)*r.u^-1*c^2*f1 + g(mu,nu)*r(la)*j(ro)*r.u^-1*f - g(mu,nu)
      *u(la)*a(ro)*r.u^-1*f + g(mu,la)*r(nu)*a(ro)*r.u^-1*r.a*f1 - g(mu,la)*
      r(nu)*a(ro)*r.u^-1*c^2*f1 + g(mu,la)*r(nu)*j(ro)*r.u^-1*f + g(mu,la)*
      u(nu)*a(ro)*f1 + g(nu,la)*r(mu)*a(ro)*r.u^-1*r.a*f1 - g(nu,la)*r(mu)*
      a(ro)*r.u^-1*c^2*f1 + g(nu,la)*r(mu)*j(ro)*r.u^-1*f + g(nu,la)*u(mu)*
      a(ro)*f1;

   ddfuu(mu,nu,la,ro) =
      2*r(mu)*r(nu)*u(la)*u(ro)*r.u^-3*r.a*c^2*f1 - r(mu)*r(nu)*u(la)*u(ro)*
      r.u^-3*r.a^2*f1 - r(mu)*r(nu)*u(la)*u(ro)*r.u^-3*c^4*f1 - 2*r(mu)*r(nu)*
      u(la)*u(ro)*r.u^-2*r.a*c^2*f2 + r(mu)*r(nu)*u(la)*u(ro)*r.u^-2*r.a^2*f2
       + r(mu)*r(nu)*u(la)*u(ro)*r.u^-2*r.j*f1 + r(mu)*r(nu)*u(la)*u(ro)*
      r.u^-2*c^4*f2 - r(mu)*r(nu)*u(la)*a(ro)*r.u^-3*r.a*f + r(mu)*r(nu)*u(la)
      *a(ro)*r.u^-3*c^2*f + 2*r(mu)*r(nu)*u(la)*a(ro)*r.u^-2*r.a*f1 - 2*r(mu)*
      r(nu)*u(la)*a(ro)*r.u^-2*c^2*f1 + r(mu)*r(nu)*u(la)*j(ro)*r.u^-2*f - 
      r(mu)*r(nu)*u(ro)*a(la)*r.u^-3*r.a*f + r(mu)*r(nu)*u(ro)*a(la)*r.u^-3*
      c^2*f + 2*r(mu)*r(nu)*u(ro)*a(la)*r.u^-2*r.a*f1 - 2*r(mu)*r(nu)*u(ro)*
      a(la)*r.u^-2*c^2*f1 + r(mu)*r(nu)*u(ro)*j(la)*r.u^-2*f + 2*r(mu)*r(nu)*
      a(la)*a(ro)*r.u^-2*f - r(mu)*u(nu)*u(la)*u(ro)*r.u^-2*r.a*f1 + r(mu)*
      u(nu)*u(la)*u(ro)*r.u^-2*c^2*f1 + r(mu)*u(nu)*u(la)*u(ro)*r.u^-1*r.a*f2
       - r(mu)*u(nu)*u(la)*u(ro)*r.u^-1*c^2*f2 - r(mu)*u(nu)*u(la)*a(ro)*
      r.u^-2*f + r(mu)*u(nu)*u(la)*a(ro)*r.u^-1*f1 - r(mu)*u(nu)*u(ro)*a(la)*
      r.u^-2*f + r(mu)*u(nu)*u(ro)*a(la)*r.u^-1*f1 + r(mu)*u(la)*u(ro)*a(nu)*
      r.u^-1*f1 - r(nu)*u(mu)*u(la)*u(ro)*r.u^-2*r.a*f1 + r(nu)*u(mu)*u(la)*
      u(ro)*r.u^-2*c^2*f1 + r(nu)*u(mu)*u(la)*u(ro)*r.u^-1*r.a*f2 - r(nu)*
      u(mu)*u(la)*u(ro)*r.u^-1*c^2*f2 - r(nu)*u(mu)*u(la)*a(ro)*r.u^-2*f + 
      r(nu)*u(mu)*u(la)*a(ro)*r.u^-1*f1 - r(nu)*u(mu)*u(ro)*a(la)*r.u^-2*f + 
      r(nu)*u(mu)*u(ro)*a(la)*r.u^-1*f1 + r(nu)*u(la)*u(ro)*a(mu)*r.u^-1*f1 + 
      u(mu)*u(nu)*u(la)*u(ro)*f2 + g(mu,nu)*u(la)*u(ro)*r.u^-1*r.a*f1 - g(mu,
      nu)*u(la)*u(ro)*r.u^-1*c^2*f1 + g(mu,nu)*u(la)*a(ro)*r.u^-1*f + g(mu,nu)
      *u(ro)*a(la)*r.u^-1*f;

   ddfua(mu,nu,la,ro) =
      2*r(mu)*r(nu)*u(la)*a(ro)*r.u^-3*r.a*c^2*f1 - r(mu)*r(nu)*u(la)*a(ro)*
      r.u^-3*r.a^2*f1 - r(mu)*r(nu)*u(la)*a(ro)*r.u^-3*c^4*f1 - 2*r(mu)*r(nu)*
      u(la)*a(ro)*r.u^-2*r.a*c^2*f2 + r(mu)*r(nu)*u(la)*a(ro)*r.u^-2*r.a^2*f2
       + r(mu)*r(nu)*u(la)*a(ro)*r.u^-2*r.j*f1 + r(mu)*r(nu)*u(la)*a(ro)*
      r.u^-2*c^4*f2 - r(mu)*r(nu)*u(la)*j(ro)*r.u^-3*r.a*f + r(mu)*r(nu)*u(la)
      *j(ro)*r.u^-3*c^2*f + 2*r(mu)*r(nu)*u(la)*j(ro)*r.u^-2*r.a*f1 - 2*r(mu)*
      r(nu)*u(la)*j(ro)*r.u^-2*c^2*f1 + r(mu)*r(nu)*u(la)*jdot(ro)*r.u^-2*f - 
      r(mu)*r(nu)*a(la)*a(ro)*r.u^-3*r.a*f + r(mu)*r(nu)*a(la)*a(ro)*r.u^-3*
      c^2*f + 2*r(mu)*r(nu)*a(la)*a(ro)*r.u^-2*r.a*f1 - 2*r(mu)*r(nu)*a(la)*
      a(ro)*r.u^-2*c^2*f1 + 2*r(mu)*r(nu)*a(la)*j(ro)*r.u^-2*f + r(mu)*r(nu)*
      a(ro)*j(la)*r.u^-2*f - r(mu)*u(nu)*u(la)*a(ro)*r.u^-2*r.a*f1 + r(mu)*
      u(nu)*u(la)*a(ro)*r.u^-2*c^2*f1 + r(mu)*u(nu)*u(la)*a(ro)*r.u^-1*r.a*f2
       - r(mu)*u(nu)*u(la)*a(ro)*r.u^-1*c^2*f2 - r(mu)*u(nu)*u(la)*j(ro)*
      r.u^-2*f + r(mu)*u(nu)*u(la)*j(ro)*r.u^-1*f1 - r(mu)*u(nu)*a(la)*a(ro)*
      r.u^-2*f + r(mu)*u(nu)*a(la)*a(ro)*r.u^-1*f1 + r(mu)*u(la)*a(nu)*a(ro)*
      r.u^-1*f1 - r(nu)*u(mu)*u(la)*a(ro)*r.u^-2*r.a*f1 + r(nu)*u(mu)*u(la)*
      a(ro)*r.u^-2*c^2*f1 + r(nu)*u(mu)*u(la)*a(ro)*r.u^-1*r.a*f2 - r(nu)*
      u(mu)*u(la)*a(ro)*r.u^-1*c^2*f2 - r(nu)*u(mu)*u(la)*j(ro)*r.u^-2*f + 
      r(nu)*u(mu)*u(la)*j(ro)*r.u^-1*f1 - r(nu)*u(mu)*a(la)*a(ro)*r.u^-2*f + 
      r(nu)*u(mu)*a(la)*a(ro)*r.u^-1*f1 + r(nu)*u(la)*a(mu)*a(ro)*r.u^-1*f1 + 
      u(mu)*u(nu)*u(la)*a(ro)*f2 + g(mu,nu)*u(la)*a(ro)*r.u^-1*r.a*f1 - g(mu,
      nu)*u(la)*a(ro)*r.u^-1*c^2*f1 + g(mu,nu)*u(la)*j(ro)*r.u^-1*f + g(mu,nu)
      *a(la)*a(ro)*r.u^-1*f;

   ddfaa(mu,nu,la,ro) =
      2*r(mu)*r(nu)*a(la)*a(ro)*r.u^-3*r.a*c^2*f1 - r(mu)*r(nu)*a(la)*a(ro)*
      r.u^-3*r.a^2*f1 - r(mu)*r(nu)*a(la)*a(ro)*r.u^-3*c^4*f1 - 2*r(mu)*r(nu)*
      a(la)*a(ro)*r.u^-2*r.a*c^2*f2 + r(mu)*r(nu)*a(la)*a(ro)*r.u^-2*r.a^2*f2
       + r(mu)*r(nu)*a(la)*a(ro)*r.u^-2*r.j*f1 + r(mu)*r(nu)*a(la)*a(ro)*
      r.u^-2*c^4*f2 - r(mu)*r(nu)*a(la)*j(ro)*r.u^-3*r.a*f + r(mu)*r(nu)*a(la)
      *j(ro)*r.u^-3*c^2*f + 2*r(mu)*r(nu)*a(la)*j(ro)*r.u^-2*r.a*f1 - 2*r(mu)*
      r(nu)*a(la)*j(ro)*r.u^-2*c^2*f1 + r(mu)*r(nu)*a(la)*jdot(ro)*r.u^-2*f - 
      r(mu)*r(nu)*a(ro)*j(la)*r.u^-3*r.a*f + r(mu)*r(nu)*a(ro)*j(la)*r.u^-3*
      c^2*f + 2*r(mu)*r(nu)*a(ro)*j(la)*r.u^-2*r.a*f1 - 2*r(mu)*r(nu)*a(ro)*
      j(la)*r.u^-2*c^2*f1 + r(mu)*r(nu)*a(ro)*jdot(la)*r.u^-2*f + 2*r(mu)*
      r(nu)*j(la)*j(ro)*r.u^-2*f - r(mu)*u(nu)*a(la)*a(ro)*r.u^-2*r.a*f1 + 
      r(mu)*u(nu)*a(la)*a(ro)*r.u^-2*c^2*f1 + r(mu)*u(nu)*a(la)*a(ro)*r.u^-1*
      r.a*f2 - r(mu)*u(nu)*a(la)*a(ro)*r.u^-1*c^2*f2 - r(mu)*u(nu)*a(la)*j(ro)
      *r.u^-2*f + r(mu)*u(nu)*a(la)*j(ro)*r.u^-1*f1 - r(mu)*u(nu)*a(ro)*j(la)*
      r.u^-2*f + r(mu)*u(nu)*a(ro)*j(la)*r.u^-1*f1 + r(mu)*a(nu)*a(la)*a(ro)*
      r.u^-1*f1 - r(nu)*u(mu)*a(la)*a(ro)*r.u^-2*r.a*f1 + r(nu)*u(mu)*a(la)*
      a(ro)*r.u^-2*c^2*f1 + r(nu)*u(mu)*a(la)*a(ro)*r.u^-1*r.a*f2 - r(nu)*
      u(mu)*a(la)*a(ro)*r.u^-1*c^2*f2 - r(nu)*u(mu)*a(la)*j(ro)*r.u^-2*f + 
      r(nu)*u(mu)*a(la)*j(ro)*r.u^-1*f1 - r(nu)*u(mu)*a(ro)*j(la)*r.u^-2*f + 
      r(nu)*u(mu)*a(ro)*j(la)*r.u^-1*f1 + r(nu)*a(mu)*a(la)*a(ro)*r.u^-1*f1 + 
      u(mu)*u(nu)*a(la)*a(ro)*f2 + g(mu,nu)*a(la)*a(ro)*r.u^-1*r.a*f1 - g(mu,
      nu)*a(la)*a(ro)*r.u^-1*c^2*f1 + g(mu,nu)*a(la)*j(ro)*r.u^-1*f + g(mu,nu)
      *a(ro)*j(la)*r.u^-1*f;

    
    *** AD HOC UNTESTED MODULES
    
    G dddfr(mu,nu,la,ro) = [f]*[r](ro);
    #call partial(la)
    #call partial(nu)
    #call partial(mu)
    #call idvec
    #call idmetric
    print +s;
    .store

   dddfr(mu,nu,la,ro) =
       + 9*r(mu)*r(nu)*r(la)*r(ro)*r.u^-5*r.a*c^4*f1
       - 9*r(mu)*r(nu)*r(la)*r(ro)*r.u^-5*r.a^2*c^2*f1
       + 3*r(mu)*r(nu)*r(la)*r(ro)*r.u^-5*r.a^3*f1
       - 3*r(mu)*r(nu)*r(la)*r(ro)*r.u^-5*c^6*f1
       - 4*r(mu)*r(nu)*r(la)*r(ro)*r.u^-4*r.a*r.j*f1
       - 9*r(mu)*r(nu)*r(la)*r(ro)*r.u^-4*r.a*c^4*f2
       + 9*r(mu)*r(nu)*r(la)*r(ro)*r.u^-4*r.a^2*c^2*f2
       - 3*r(mu)*r(nu)*r(la)*r(ro)*r.u^-4*r.a^3*f2
       + 4*r(mu)*r(nu)*r(la)*r(ro)*r.u^-4*r.j*c^2*f1
       + 3*r(mu)*r(nu)*r(la)*r(ro)*r.u^-4*c^6*f2
       + 3*r(mu)*r(nu)*r(la)*r(ro)*r.u^-3*r.a*r.j*f2
       + 3*r(mu)*r(nu)*r(la)*r(ro)*r.u^-3*r.a*c^4*f3
       - 3*r(mu)*r(nu)*r(la)*r(ro)*r.u^-3*r.a^2*c^2*f3
       + r(mu)*r(nu)*r(la)*r(ro)*r.u^-3*r.a^3*f3
       - 3*r(mu)*r(nu)*r(la)*r(ro)*r.u^-3*r.j*c^2*f2
       + r(mu)*r(nu)*r(la)*r(ro)*r.u^-3*r.jdot*f1
       + r(mu)*r(nu)*r(la)*r(ro)*r.u^-3*a.a*f1
       - r(mu)*r(nu)*r(la)*r(ro)*r.u^-3*c^6*f3
       + 6*r(mu)*r(nu)*r(la)*u(ro)*r.u^-5*r.a*c^2*f
       - 3*r(mu)*r(nu)*r(la)*u(ro)*r.u^-5*r.a^2*f
       - 3*r(mu)*r(nu)*r(la)*u(ro)*r.u^-5*c^4*f
       - 12*r(mu)*r(nu)*r(la)*u(ro)*r.u^-4*r.a*c^2*f1
       + 6*r(mu)*r(nu)*r(la)*u(ro)*r.u^-4*r.a^2*f1
       + r(mu)*r(nu)*r(la)*u(ro)*r.u^-4*r.j*f
       + 6*r(mu)*r(nu)*r(la)*u(ro)*r.u^-4*c^4*f1
       + 6*r(mu)*r(nu)*r(la)*u(ro)*r.u^-3*r.a*c^2*f2
       - 3*r(mu)*r(nu)*r(la)*u(ro)*r.u^-3*r.a^2*f2
       - 3*r(mu)*r(nu)*r(la)*u(ro)*r.u^-3*r.j*f1
       - 3*r(mu)*r(nu)*r(la)*u(ro)*r.u^-3*c^4*f2
       + 3*r(mu)*r(nu)*r(la)*a(ro)*r.u^-4*r.a*f
       - 3*r(mu)*r(nu)*r(la)*a(ro)*r.u^-4*c^2*f
       - 3*r(mu)*r(nu)*r(la)*a(ro)*r.u^-3*r.a*f1
       + 3*r(mu)*r(nu)*r(la)*a(ro)*r.u^-3*c^2*f1
       - r(mu)*r(nu)*r(la)*j(ro)*r.u^-3*f
       - 6*r(mu)*r(nu)*r(ro)*u(la)*r.u^-4*r.a*c^2*f1
       + 3*r(mu)*r(nu)*r(ro)*u(la)*r.u^-4*r.a^2*f1
       + 3*r(mu)*r(nu)*r(ro)*u(la)*r.u^-4*c^4*f1
       + 6*r(mu)*r(nu)*r(ro)*u(la)*r.u^-3*r.a*c^2*f2
       - 3*r(mu)*r(nu)*r(ro)*u(la)*r.u^-3*r.a^2*f2
       - 2*r(mu)*r(nu)*r(ro)*u(la)*r.u^-3*r.j*f1
       - 3*r(mu)*r(nu)*r(ro)*u(la)*r.u^-3*c^4*f2
       - 2*r(mu)*r(nu)*r(ro)*u(la)*r.u^-2*r.a*c^2*f3
       + r(mu)*r(nu)*r(ro)*u(la)*r.u^-2*r.a^2*f3
       + r(mu)*r(nu)*r(ro)*u(la)*r.u^-2*r.j*f2
       + r(mu)*r(nu)*r(ro)*u(la)*r.u^-2*c^4*f3
       - 2*r(mu)*r(nu)*r(ro)*a(la)*r.u^-3*r.a*f1
       + 2*r(mu)*r(nu)*r(ro)*a(la)*r.u^-3*c^2*f1
       + 2*r(mu)*r(nu)*r(ro)*a(la)*r.u^-2*r.a*f2
       - 2*r(mu)*r(nu)*r(ro)*a(la)*r.u^-2*c^2*f2
       + r(mu)*r(nu)*r(ro)*j(la)*r.u^-2*f1
       - 3*r(mu)*r(nu)*u(la)*u(ro)*r.u^-4*r.a*f
       + 3*r(mu)*r(nu)*u(la)*u(ro)*r.u^-4*c^2*f
       + 5*r(mu)*r(nu)*u(la)*u(ro)*r.u^-3*r.a*f1
       - 5*r(mu)*r(nu)*u(la)*u(ro)*r.u^-3*c^2*f1
       - 2*r(mu)*r(nu)*u(la)*u(ro)*r.u^-2*r.a*f2
       + 2*r(mu)*r(nu)*u(la)*u(ro)*r.u^-2*c^2*f2
       + 2*r(mu)*r(nu)*u(la)*a(ro)*r.u^-3*f
       - r(mu)*r(nu)*u(la)*a(ro)*r.u^-2*f1
       + r(mu)*r(nu)*u(ro)*a(la)*r.u^-3*f
       - 2*r(mu)*r(nu)*u(ro)*a(la)*r.u^-2*f1
       - 6*r(mu)*r(la)*r(ro)*u(nu)*r.u^-4*r.a*c^2*f1
       + 3*r(mu)*r(la)*r(ro)*u(nu)*r.u^-4*r.a^2*f1
       + 3*r(mu)*r(la)*r(ro)*u(nu)*r.u^-4*c^4*f1
       + 6*r(mu)*r(la)*r(ro)*u(nu)*r.u^-3*r.a*c^2*f2
       - 3*r(mu)*r(la)*r(ro)*u(nu)*r.u^-3*r.a^2*f2
       - 2*r(mu)*r(la)*r(ro)*u(nu)*r.u^-3*r.j*f1
       - 3*r(mu)*r(la)*r(ro)*u(nu)*r.u^-3*c^4*f2
       - 2*r(mu)*r(la)*r(ro)*u(nu)*r.u^-2*r.a*c^2*f3
       + r(mu)*r(la)*r(ro)*u(nu)*r.u^-2*r.a^2*f3
       + r(mu)*r(la)*r(ro)*u(nu)*r.u^-2*r.j*f2
       + r(mu)*r(la)*r(ro)*u(nu)*r.u^-2*c^4*f3
       - 2*r(mu)*r(la)*r(ro)*a(nu)*r.u^-3*r.a*f1
       + 2*r(mu)*r(la)*r(ro)*a(nu)*r.u^-3*c^2*f1
       + 2*r(mu)*r(la)*r(ro)*a(nu)*r.u^-2*r.a*f2
       - 2*r(mu)*r(la)*r(ro)*a(nu)*r.u^-2*c^2*f2
       + r(mu)*r(la)*r(ro)*j(nu)*r.u^-2*f1
       - 3*r(mu)*r(la)*u(nu)*u(ro)*r.u^-4*r.a*f
       + 3*r(mu)*r(la)*u(nu)*u(ro)*r.u^-4*c^2*f
       + 5*r(mu)*r(la)*u(nu)*u(ro)*r.u^-3*r.a*f1
       - 5*r(mu)*r(la)*u(nu)*u(ro)*r.u^-3*c^2*f1
       - 2*r(mu)*r(la)*u(nu)*u(ro)*r.u^-2*r.a*f2
       + 2*r(mu)*r(la)*u(nu)*u(ro)*r.u^-2*c^2*f2
       + 2*r(mu)*r(la)*u(nu)*a(ro)*r.u^-3*f
       - r(mu)*r(la)*u(nu)*a(ro)*r.u^-2*f1
       + r(mu)*r(la)*u(ro)*a(nu)*r.u^-3*f
       - 2*r(mu)*r(la)*u(ro)*a(nu)*r.u^-2*f1
       + 2*r(mu)*r(ro)*u(nu)*u(la)*r.u^-3*r.a*f1
       - 2*r(mu)*r(ro)*u(nu)*u(la)*r.u^-3*c^2*f1
       - 2*r(mu)*r(ro)*u(nu)*u(la)*r.u^-2*r.a*f2
       + 2*r(mu)*r(ro)*u(nu)*u(la)*r.u^-2*c^2*f2
       + r(mu)*r(ro)*u(nu)*u(la)*r.u^-1*r.a*f3
       - r(mu)*r(ro)*u(nu)*u(la)*r.u^-1*c^2*f3
       - r(mu)*r(ro)*u(nu)*a(la)*r.u^-2*f1
       + r(mu)*r(ro)*u(nu)*a(la)*r.u^-1*f2
       - r(mu)*r(ro)*u(la)*a(nu)*r.u^-2*f1
       + r(mu)*r(ro)*u(la)*a(nu)*r.u^-1*f2
       - 2*r(mu)*u(nu)*u(la)*u(ro)*r.u^-3*f
       + 2*r(mu)*u(nu)*u(la)*u(ro)*r.u^-2*f1
       - r(mu)*u(nu)*u(la)*u(ro)*r.u^-1*f2
       - 6*r(nu)*r(la)*r(ro)*u(mu)*r.u^-4*r.a*c^2*f1
       + 3*r(nu)*r(la)*r(ro)*u(mu)*r.u^-4*r.a^2*f1
       + 3*r(nu)*r(la)*r(ro)*u(mu)*r.u^-4*c^4*f1
       + 6*r(nu)*r(la)*r(ro)*u(mu)*r.u^-3*r.a*c^2*f2
       - 3*r(nu)*r(la)*r(ro)*u(mu)*r.u^-3*r.a^2*f2
       - 2*r(nu)*r(la)*r(ro)*u(mu)*r.u^-3*r.j*f1
       - 3*r(nu)*r(la)*r(ro)*u(mu)*r.u^-3*c^4*f2
       - 2*r(nu)*r(la)*r(ro)*u(mu)*r.u^-2*r.a*c^2*f3
       + r(nu)*r(la)*r(ro)*u(mu)*r.u^-2*r.a^2*f3
       + r(nu)*r(la)*r(ro)*u(mu)*r.u^-2*r.j*f2
       + r(nu)*r(la)*r(ro)*u(mu)*r.u^-2*c^4*f3
       - 2*r(nu)*r(la)*r(ro)*a(mu)*r.u^-3*r.a*f1
       + 2*r(nu)*r(la)*r(ro)*a(mu)*r.u^-3*c^2*f1
       + 2*r(nu)*r(la)*r(ro)*a(mu)*r.u^-2*r.a*f2
       - 2*r(nu)*r(la)*r(ro)*a(mu)*r.u^-2*c^2*f2
       + r(nu)*r(la)*r(ro)*j(mu)*r.u^-2*f1
       - 3*r(nu)*r(la)*u(mu)*u(ro)*r.u^-4*r.a*f
       + 3*r(nu)*r(la)*u(mu)*u(ro)*r.u^-4*c^2*f
       + 5*r(nu)*r(la)*u(mu)*u(ro)*r.u^-3*r.a*f1
       - 5*r(nu)*r(la)*u(mu)*u(ro)*r.u^-3*c^2*f1
       - 2*r(nu)*r(la)*u(mu)*u(ro)*r.u^-2*r.a*f2
       + 2*r(nu)*r(la)*u(mu)*u(ro)*r.u^-2*c^2*f2
       + 2*r(nu)*r(la)*u(mu)*a(ro)*r.u^-3*f
       - r(nu)*r(la)*u(mu)*a(ro)*r.u^-2*f1
       + r(nu)*r(la)*u(ro)*a(mu)*r.u^-3*f
       - 2*r(nu)*r(la)*u(ro)*a(mu)*r.u^-2*f1
       + 2*r(nu)*r(ro)*u(mu)*u(la)*r.u^-3*r.a*f1
       - 2*r(nu)*r(ro)*u(mu)*u(la)*r.u^-3*c^2*f1
       - 2*r(nu)*r(ro)*u(mu)*u(la)*r.u^-2*r.a*f2
       + 2*r(nu)*r(ro)*u(mu)*u(la)*r.u^-2*c^2*f2
       + r(nu)*r(ro)*u(mu)*u(la)*r.u^-1*r.a*f3
       - r(nu)*r(ro)*u(mu)*u(la)*r.u^-1*c^2*f3
       - r(nu)*r(ro)*u(mu)*a(la)*r.u^-2*f1
       + r(nu)*r(ro)*u(mu)*a(la)*r.u^-1*f2
       - r(nu)*r(ro)*u(la)*a(mu)*r.u^-2*f1
       + r(nu)*r(ro)*u(la)*a(mu)*r.u^-1*f2
       - 2*r(nu)*u(mu)*u(la)*u(ro)*r.u^-3*f
       + 2*r(nu)*u(mu)*u(la)*u(ro)*r.u^-2*f1
       - r(nu)*u(mu)*u(la)*u(ro)*r.u^-1*f2
       + 2*r(la)*r(ro)*u(mu)*u(nu)*r.u^-3*r.a*f1
       - 2*r(la)*r(ro)*u(mu)*u(nu)*r.u^-3*c^2*f1
       - 2*r(la)*r(ro)*u(mu)*u(nu)*r.u^-2*r.a*f2
       + 2*r(la)*r(ro)*u(mu)*u(nu)*r.u^-2*c^2*f2
       + r(la)*r(ro)*u(mu)*u(nu)*r.u^-1*r.a*f3
       - r(la)*r(ro)*u(mu)*u(nu)*r.u^-1*c^2*f3
       - r(la)*r(ro)*u(mu)*a(nu)*r.u^-2*f1
       + r(la)*r(ro)*u(mu)*a(nu)*r.u^-1*f2
       - r(la)*r(ro)*u(nu)*a(mu)*r.u^-2*f1
       + r(la)*r(ro)*u(nu)*a(mu)*r.u^-1*f2
       - 2*r(la)*u(mu)*u(nu)*u(ro)*r.u^-3*f
       + 2*r(la)*u(mu)*u(nu)*u(ro)*r.u^-2*f1
       - r(la)*u(mu)*u(nu)*u(ro)*r.u^-1*f2
       + r(ro)*u(mu)*u(nu)*u(la)*f3
       + 2*g(mu,nu)*r(la)*r(ro)*r.u^-3*r.a*c^2*f1
       - g(mu,nu)*r(la)*r(ro)*r.u^-3*r.a^2*f1
       - g(mu,nu)*r(la)*r(ro)*r.u^-3*c^4*f1
       - 2*g(mu,nu)*r(la)*r(ro)*r.u^-2*r.a*c^2*f2
       + g(mu,nu)*r(la)*r(ro)*r.u^-2*r.a^2*f2
       + g(mu,nu)*r(la)*r(ro)*r.u^-2*r.j*f1
       + g(mu,nu)*r(la)*r(ro)*r.u^-2*c^4*f2
       + g(mu,nu)*r(la)*u(ro)*r.u^-3*r.a*f
       - g(mu,nu)*r(la)*u(ro)*r.u^-3*c^2*f
       - 2*g(mu,nu)*r(la)*u(ro)*r.u^-2*r.a*f1
       + 2*g(mu,nu)*r(la)*u(ro)*r.u^-2*c^2*f1
       - g(mu,nu)*r(la)*a(ro)*r.u^-2*f
       - g(mu,nu)*r(ro)*u(la)*r.u^-2*r.a*f1
       + g(mu,nu)*r(ro)*u(la)*r.u^-2*c^2*f1
       + g(mu,nu)*r(ro)*u(la)*r.u^-1*r.a*f2
       - g(mu,nu)*r(ro)*u(la)*r.u^-1*c^2*f2
       + g(mu,nu)*r(ro)*a(la)*r.u^-1*f1
       + g(mu,nu)*u(la)*u(ro)*r.u^-2*f
       - g(mu,nu)*u(la)*u(ro)*r.u^-1*f1
       + g(mu,nu)*g(la,ro)*r.u^-1*r.a*f1
       - g(mu,nu)*g(la,ro)*r.u^-1*c^2*f1
       + 2*g(mu,la)*r(nu)*r(ro)*r.u^-3*r.a*c^2*f1
       - g(mu,la)*r(nu)*r(ro)*r.u^-3*r.a^2*f1
       - g(mu,la)*r(nu)*r(ro)*r.u^-3*c^4*f1
       - 2*g(mu,la)*r(nu)*r(ro)*r.u^-2*r.a*c^2*f2
       + g(mu,la)*r(nu)*r(ro)*r.u^-2*r.a^2*f2
       + g(mu,la)*r(nu)*r(ro)*r.u^-2*r.j*f1
       + g(mu,la)*r(nu)*r(ro)*r.u^-2*c^4*f2
       + g(mu,la)*r(nu)*u(ro)*r.u^-3*r.a*f
       - g(mu,la)*r(nu)*u(ro)*r.u^-3*c^2*f
       - 2*g(mu,la)*r(nu)*u(ro)*r.u^-2*r.a*f1
       + 2*g(mu,la)*r(nu)*u(ro)*r.u^-2*c^2*f1
       - g(mu,la)*r(nu)*a(ro)*r.u^-2*f
       - g(mu,la)*r(ro)*u(nu)*r.u^-2*r.a*f1
       + g(mu,la)*r(ro)*u(nu)*r.u^-2*c^2*f1
       + g(mu,la)*r(ro)*u(nu)*r.u^-1*r.a*f2
       - g(mu,la)*r(ro)*u(nu)*r.u^-1*c^2*f2
       + g(mu,la)*r(ro)*a(nu)*r.u^-1*f1
       + g(mu,la)*u(nu)*u(ro)*r.u^-2*f
       - g(mu,la)*u(nu)*u(ro)*r.u^-1*f1
       + g(mu,la)*g(nu,ro)*r.u^-1*r.a*f1
       - g(mu,la)*g(nu,ro)*r.u^-1*c^2*f1
       + 2*g(mu,ro)*r(nu)*r(la)*r.u^-3*r.a*c^2*f1
       - g(mu,ro)*r(nu)*r(la)*r.u^-3*r.a^2*f1
       - g(mu,ro)*r(nu)*r(la)*r.u^-3*c^4*f1
       - 2*g(mu,ro)*r(nu)*r(la)*r.u^-2*r.a*c^2*f2
       + g(mu,ro)*r(nu)*r(la)*r.u^-2*r.a^2*f2
       + g(mu,ro)*r(nu)*r(la)*r.u^-2*r.j*f1
       + g(mu,ro)*r(nu)*r(la)*r.u^-2*c^4*f2
       - g(mu,ro)*r(nu)*u(la)*r.u^-2*r.a*f1
       + g(mu,ro)*r(nu)*u(la)*r.u^-2*c^2*f1
       + g(mu,ro)*r(nu)*u(la)*r.u^-1*r.a*f2
       - g(mu,ro)*r(nu)*u(la)*r.u^-1*c^2*f2
       + g(mu,ro)*r(nu)*a(la)*r.u^-1*f1
       - g(mu,ro)*r(la)*u(nu)*r.u^-2*r.a*f1
       + g(mu,ro)*r(la)*u(nu)*r.u^-2*c^2*f1
       + g(mu,ro)*r(la)*u(nu)*r.u^-1*r.a*f2
       - g(mu,ro)*r(la)*u(nu)*r.u^-1*c^2*f2
       + g(mu,ro)*r(la)*a(nu)*r.u^-1*f1
       + g(mu,ro)*u(nu)*u(la)*f2
       + g(mu,ro)*g(nu,la)*r.u^-1*r.a*f1
       - g(mu,ro)*g(nu,la)*r.u^-1*c^2*f1
       + 2*g(nu,la)*r(mu)*r(ro)*r.u^-3*r.a*c^2*f1
       - g(nu,la)*r(mu)*r(ro)*r.u^-3*r.a^2*f1
       - g(nu,la)*r(mu)*r(ro)*r.u^-3*c^4*f1
       - 2*g(nu,la)*r(mu)*r(ro)*r.u^-2*r.a*c^2*f2
       + g(nu,la)*r(mu)*r(ro)*r.u^-2*r.a^2*f2
       + g(nu,la)*r(mu)*r(ro)*r.u^-2*r.j*f1
       + g(nu,la)*r(mu)*r(ro)*r.u^-2*c^4*f2
       + g(nu,la)*r(mu)*u(ro)*r.u^-3*r.a*f
       - g(nu,la)*r(mu)*u(ro)*r.u^-3*c^2*f
       - 2*g(nu,la)*r(mu)*u(ro)*r.u^-2*r.a*f1
       + 2*g(nu,la)*r(mu)*u(ro)*r.u^-2*c^2*f1
       - g(nu,la)*r(mu)*a(ro)*r.u^-2*f
       - g(nu,la)*r(ro)*u(mu)*r.u^-2*r.a*f1
       + g(nu,la)*r(ro)*u(mu)*r.u^-2*c^2*f1
       + g(nu,la)*r(ro)*u(mu)*r.u^-1*r.a*f2
       - g(nu,la)*r(ro)*u(mu)*r.u^-1*c^2*f2
       + g(nu,la)*r(ro)*a(mu)*r.u^-1*f1
       + g(nu,la)*u(mu)*u(ro)*r.u^-2*f
       - g(nu,la)*u(mu)*u(ro)*r.u^-1*f1
       + 2*g(nu,ro)*r(mu)*r(la)*r.u^-3*r.a*c^2*f1
       - g(nu,ro)*r(mu)*r(la)*r.u^-3*r.a^2*f1
       - g(nu,ro)*r(mu)*r(la)*r.u^-3*c^4*f1
       - 2*g(nu,ro)*r(mu)*r(la)*r.u^-2*r.a*c^2*f2
       + g(nu,ro)*r(mu)*r(la)*r.u^-2*r.a^2*f2
       + g(nu,ro)*r(mu)*r(la)*r.u^-2*r.j*f1
       + g(nu,ro)*r(mu)*r(la)*r.u^-2*c^4*f2
       - g(nu,ro)*r(mu)*u(la)*r.u^-2*r.a*f1
       + g(nu,ro)*r(mu)*u(la)*r.u^-2*c^2*f1
       + g(nu,ro)*r(mu)*u(la)*r.u^-1*r.a*f2
       - g(nu,ro)*r(mu)*u(la)*r.u^-1*c^2*f2
       + g(nu,ro)*r(mu)*a(la)*r.u^-1*f1
       - g(nu,ro)*r(la)*u(mu)*r.u^-2*r.a*f1
       + g(nu,ro)*r(la)*u(mu)*r.u^-2*c^2*f1
       + g(nu,ro)*r(la)*u(mu)*r.u^-1*r.a*f2
       - g(nu,ro)*r(la)*u(mu)*r.u^-1*c^2*f2
       + g(nu,ro)*r(la)*a(mu)*r.u^-1*f1
       + g(nu,ro)*u(mu)*u(la)*f2
       + 2*g(la,ro)*r(mu)*r(nu)*r.u^-3*r.a*c^2*f1
       - g(la,ro)*r(mu)*r(nu)*r.u^-3*r.a^2*f1
       - g(la,ro)*r(mu)*r(nu)*r.u^-3*c^4*f1
       - 2*g(la,ro)*r(mu)*r(nu)*r.u^-2*r.a*c^2*f2
       + g(la,ro)*r(mu)*r(nu)*r.u^-2*r.a^2*f2
       + g(la,ro)*r(mu)*r(nu)*r.u^-2*r.j*f1
       + g(la,ro)*r(mu)*r(nu)*r.u^-2*c^4*f2
       - g(la,ro)*r(mu)*u(nu)*r.u^-2*r.a*f1
       + g(la,ro)*r(mu)*u(nu)*r.u^-2*c^2*f1
       + g(la,ro)*r(mu)*u(nu)*r.u^-1*r.a*f2
       - g(la,ro)*r(mu)*u(nu)*r.u^-1*c^2*f2
       + g(la,ro)*r(mu)*a(nu)*r.u^-1*f1
       - g(la,ro)*r(nu)*u(mu)*r.u^-2*r.a*f1
       + g(la,ro)*r(nu)*u(mu)*r.u^-2*c^2*f1
       + g(la,ro)*r(nu)*u(mu)*r.u^-1*r.a*f2
       - g(la,ro)*r(nu)*u(mu)*r.u^-1*c^2*f2
       + g(la,ro)*r(nu)*a(mu)*r.u^-1*f1
       + g(la,ro)*u(mu)*u(nu)*f2
      ;

    
    save derivs.sav;
    .end
  0.04 sec out of 0.04 sec