FORM 4.0 (Jun 30 2012) 64-bits                   Run: Sat Aug  3 17:45:52 2013
    *   Title:  SMALL CHARGED SPHERE:  Lorentz force
    *           F(mu,nu)*U(nu) singular orders
    *    File:  lfsing.frm
    *  Author:  David N. Williams
    * License:  Creative Commons Attribution-Share Alike
    * Started:  November 7, 1987 (Schoonschip)
    *           July 9, 2012     (Form)
    * Revised:  July 9,10,27-29, 2012
    * Revised:  August 2, 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:  TAU1, TAU2, TAU3 from dtau.sav
    *           RdotA from rdota.sav
    *           [1/RdotU^2], [1/RdotU^3] from denoms.sav
    *  Output:  lfsing in lfsing.sav
    
    #if 0
    
    The electromagnetic units are rationalized Gaussian, the same as
    Jackson.
    
    The quantity F(mu,nu) below is the proper electromagnetic field
    density 𝔽(x,β).  It remains to be integrated over the proper
    three-position, so the field itself is
    
      F(x) = ∫d³β ρ₀(β) 𝔽(x,β),   x = y(τ,α),   ∫d³β ρ₀(β) = Q
    .
    
    The Lorentz force density, which we call the Lorentz force
    field, is
    
      f(x) = F(x)⋅j(x)/c,   j(x) = ρ₀(α)/Jγ u(τ,α),
    
    where γ = u⁰(τ,α), and J is the Jacobian of the transfomation
    from the proper coordinates (τ,α) to x.
    
    The quantity fsingular is the singular part of the proper
    density 𝕗(x,β) of the Lorentz force field, without the integral
    over β in F and without the factor ρ₀(α)/Jγ in j.  The singular
    part of f(x) is
    
      ∫d³β ρ₀(β) 𝕗(x,β)
    
        =  ∫d³β ρ₀(β) 𝔽[y(τ,α),β]⋅u(τ,α) ρ₀(α)/Jγc
    
        = ∫d³β ρ₀(β) ρ₀(α)/Jγ fsingular
    
    #endif
    
    *** DECLARATIONS
    
    OFF statistics;
     
    S eps,c,snorm,p,L,T;
    V u,a,j,n;
    T d1y,d2y,d3y,d1u,d2u,d1a;
    T eta;
    
    Load dtau.sav;
 TAU1 loaded
 TAU2 loaded
 TAU3 loaded
    Load rdota.sav;
 RdotA loaded
    Load denoms.sav;
 [1/RdotU^2] loaded
 [1/RdotU^3] loaded
    
    S dtau;
    S [1/d^2],[1/d^3];
    I mu,nu,la,l,m;
    CF R,U,A;
    
    *** DEBUGGING
    
    #procedure try
    .sort
    B eps,snorm;
    print +s;
    .end
    #endprocedure
    
    *** MODULES
    
    * leading nonsingular order in this section:  eps^3
    
    L F(mu,nu) = ( R(mu)*U(nu) - R(nu)*U(mu) ) * [1/d^3] * ( c^2 - RdotA )
    
                 + ( R(mu)*A(nu) - R(nu)*A(mu) ) * [1/d^2];
    
    
    Id R(mu?) = - d1y(mu,n)*eps - d2y(mu,n,n)*eps^2/2
    
                    - d3y(mu,n,n,n)*eps^3/6
    
                - u(mu)*dtau - d1u(mu,n)*eps*dtau
    
                    - d2u(mu,n,n)*eps^2*dtau/2
    
                - a(mu)*dtau^2/2 - d1a(mu,n)*eps*dtau^2/2
    
                - j(mu)*dtau^3/6;
    
    Id U(mu?) =   u(mu) + d1u(mu,n)*eps + d2u(mu,n,n)*eps^2/2
    
                + a(mu)*dtau + d1a(mu,n)*eps*dtau
    
                + j(mu)*dtau^2/2;
    
    Id A(mu?) =   a(mu) + d1a(mu,n)*eps
    
                + j(mu)*dtau;
    
    Id eps^4 = 0;
    
    .sort
    
    * Lorentz force density
    G lfsing(mu) = F(mu,nu)*u(nu)/c;
    Id eps^4 = 0;
    Id eps^3 * dtau = 0;
    Id eps^2 * dtau^2 = 0;
    Id eps * dtau^3 = 0;
    Id dtau^4 = 0;
    
    Id u.u = c^2;
    Id u.a = 0;
    Id u.j = -a.a;
    
    Id d1u(u,l?) = 0;
    Id d1a(u,l?) = -d1u(a,l);
    Id, d2u(u,l?,m?) = -d1u(la,l)*d1u(la,m);
    Sum la;
    
    Id dtau = TAU1*eps + TAU2*eps^2 + TAU3*eps^3;
    Id eps^4 = 0;
    .sort
    
    Id [1/d^2] = [1/RdotU^2];
    Id [1/d^3] = [1/RdotU^3];
    .sort
    
    * kill superfluous order from RdotA
    Id eps = 0;
    
    * ready to kill nonsingular order
    multiply eps;
    
    * kill nonsingular order
    Id eps = 0;
    
    * restore eps powers
    multiply eps^-1;
    
    Id d1y(mu,n) = eta(mu,n) + d1y(u,n)*u(mu)*c^-2;
    .sort
    
    * check units
    
    Hide lfsing;
    L lfsingunits = lfsing(mu);
    Id c^p? = (L/T)^p;
    Id a(mu) = L/T^2;
    Id u(mu) = L/T;
    Id eta(mu,l?) = 1;
    Id d1y(u,l?) = L/T;
    Id d1u(mu?,l?) = 1/T;
    Id d1y(a,l?) = L/T^2;
    Id d1y(mu?,l?) = 1;
    Id d2y(u,l?,m?) = 1/T;
    Id d2y(mu?,l?,m?) = 1/L;
    Id eps^p? = L^p;
    Id snorm^p? = 1;
    .sort
    
    Unhide lfsing;
    B eps,snorm;
    print +s;
    .store

   F(mu,nu) =
       - 5/2*L^-2
          + 3/2*u(nu)*L^-3*T
          + 2*a(nu)*L^-3*T^2
         ;

   lfsing(mu) =

       + eps^-2*snorm^-3 * (
          - eta(mu,n)
          )

       + eps^-1*snorm^-5 * (
          - 3/2*d1y(u,n)*d1y(u,n)*d1y(a,n)*eta(mu,n)*c^-4
          + 3*d1y(u,n)*d1y(N1_?,n)*d1u(N1_?,n)*eta(mu,n)*c^-2
          + 3/2*d1y(u,n)*d2y(u,n,n)*eta(mu,n)*c^-2
          - 3/2*d1y(N1_?,n)*d2y(N1_?,n,n)*eta(mu,n)
          )

       + eps^-1*snorm^-3 * (
          - 1/2*d1y(u,n)*d1y(u,n)*a(mu)*c^-4
          + d1y(u,n)*d1u(mu,n)*c^-2
          + 1/2*d1y(a,n)*eta(mu,n)*c^-2
          + 1/2*d2y(u,n,n)*u(mu)*c^-2
          - 1/2*d2y(mu,n,n)
          )

       + eps^-1*snorm^-1 * (
          - 1/2*a(mu)*c^-2
          );

   lfsingunits =
       + L^-2
         ;

    
    Save lfsing.sav lfsing;
    .end
  5.18 sec out of 5.18 sec