FORM 4.0 (Jun 30 2012) 64-bits                   Run: Sat Aug  3 17:46:03 2013
    *   Title:  Elastic force density
    *    File:  elasforce.frm
    *  Author:  David N. Williams
    * License:  Creative Commons Attribution-Share Alike
    * Started:  July 23, 2012
    * Revised:  July 23-31, 2012
    *           August 1-4, 2012
    *           October 17, 2012  improved comments
    *           March 25, 2013    renamed F as U
    *           March 26, 2013    removed canonicalization experiment 
    *           July 3, 2013
    *           August 2, 2013    renamed kappa and zeta as kap and zet
    *           August 3, 2013    renamed [dL.dU/dzet] and [dR.dU/dzet] as
    *                             [d2U/dalp.dzet] and [d2U/dzet.dalp]
    *
    *   Copyright (C) 2012, 2013 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.
    *
    * Purely cosmetic changes may not be reflected in the revision
    * dates above.
    *
    *   Input:  none
    *  Output:  divE in elasforce.sav, missing factor 1/Jγ
    *           elasforce.log
    
    #if 0
    
    This program works out the four-divergence of the relativistic
    elastic energy-momentum tensor E for a continuous matter
    distribution, for a general proper potential energy density
    function
    
      U = U(ζ,α).
    
    Here α is the invariant reference three-position of an elastic
    body point, and ζ is the invariant 3x3 proper strain tensor,
    defined below.
    
    See elasemPE.frm for the use of the results to calculate the
    compensation of the elastic part of the Lorentz force density on
    a small charge.
    
    The world line four-position y is assumed to be given as a
    function of α and proper time τ.  Everything else is expressed
    in terms of y by the formulas below.  Page numbers refer to
    handwritten notes from September, 2007, which follow our
    pedagogical article, "The Elastic Energy-Momentum Tensor in
    Special Relativity," Ann. Phys. 196, 345–360 (1989), with a
    few changes in notation.
    
      y ≡ y(τ,α), u ≡ ∂y(τ,α)/∂τ, a ≡ ∂u(τ,α)/∂τ
      ηᵢ ≡ ηᵢ(τ,α) = ∂ᵢy - uu⋅∂ᵢy/c²
      ∂ᵢ ≡ ∂/∂αⁱ
    
      u⋅u = c²,  u⋅a = 0,  u⋅ηᵢ = 0,  u⋅∂ᵢu = 0
    
      ζᵢⱼ ≡ -ηᵢ⋅ηⱼ
    
      Jγ E = uu/c² U - (ηᵢηⱼ + ηⱼηᵢ)∂U/∂ζᵢⱼ (p. 29)
    
      Jγ div E = a/c² U
        + ∂U/∂ζᵢⱼ ( u/c² ∂ζᵢⱼ/∂τ + ηᵢ/c ∂κⱼ/∂τ + η�
    �/c ∂κᵢ/∂τ
                     - Dᵢ ηⱼ - Dⱼ ηᵢ )
        - ηᵢ Dⱼ ∂U/∂ζᵢⱼ - ηⱼ Dᵢ ∂U/∂ζᵢⱼ (p. 29)
    
      κᵢ ≡ u⋅∂ᵢy/c
    
      ∂κᵢ/∂τ = a⋅∂ᵢy/c = a⋅ηᵢ/c (p. 8)
    
      Dᵢ ≡ ∂ᵢ - κᵢ/c ∂/∂τ
    
      ηᵢ = ∂ᵢy - u/c κᵢ
    
      ∂ζᵢⱼ/∂τ = κᵢ∂κⱼ/∂τ + κⱼ∂κᵢ/∂τ - ∂ᵢu⋅∂�
    ��y - ∂ⱼu⋅∂ᵢy (p. 9, cf. p. 10)
    
      ∂ζᵢⱼ/∂αʳ = ∂ᵣζᵢⱼ = - ∂ᵣηᵢ⋅ηⱼ - ∂ᵣηⱼ⋅�
    �ᵢ
    
      ∂ᵢκⱼ = ∂ᵢu⋅∂ⱼy/c + u⋅∂ᵢ∂ⱼy/c = ∂ᵢu⋅ηⱼ/c + 
    u⋅∂ᵢ∂ⱼy/c
    
      ∂ᵢηⱼ = ∂ᵢ∂ⱼy - ∂ᵢu u⋅∂ⱼy/c² - u ∂ᵢu⋅∂ⱼy/c�
    � - uu⋅∂ᵢ∂ⱼy/c²
    
      ∂ηᵢ/∂τ = ∂ᵢu - au⋅∂ᵢy/c² - ua⋅∂ᵢy/c²
    
      Dᵢηⱼ = ∂ᵢ∂ⱼy - ∂ᵢκⱼ u/c - κᵢ∂ⱼu/c - κⱼ∂ᵢ u/c
     + κᵢ∂κⱼ/∂τ u/c²
             + κᵢκⱼ a/c² (p. 30)
    
      Dᵢ ∂U/∂ζᵢⱼ = ∂²U/∂ζᵢⱼ∂ζᵤᵥ ∂ᵢζᵤᵥ - κᵢ/c
     ∂²U/∂ζᵢⱼ∂ζᵤᵥ ∂ζᵤᵥ/∂τ
        + ∂²U/∂ζᵢⱼ∂αᵢ (left div, p. 12)
    
      Dⱼ ∂U/∂ζᵢⱼ = ∂²U/∂ζᵢⱼ∂ζᵤᵥ ∂ⱼζᵤᵥ - κⱼ/c
     ∂²U/∂ζᵢⱼ∂ζᵤᵥ ∂ζᵤᵥ/∂τ
        + ∂²U/∂ζᵢⱼ∂αⱼ (right div)
    
    The invariant Jγ is the Jacobian of the four-dimensional
    transformation from (τ,α) to y.
    
    The α in the proper potential energy density U(ζ,α) corresponds
    to explicit α dependence.  U has no explicit τ dependence, but
    has implicit dependence on both α and τ through ζ.
    
    The derivatives of U with respect to ζ are to be understood as
    generic matrix derivatives evaluated at the symmetric matrix ζ.
    Neither the first nor the second derivative is necessarily a
    symmetric function of ζ.  Being a symmetric function implies
    symmetry of the derivatives in the index pairs, but there are
    counterexamples to index symmetry when the function is not
    symmetric in ζ.  The index pairs do get explicitly symmetrized
    in the result for the four-divergence.
    
    The correspondence between Form names and mathematical
    quantities is the following:
    
      eta(mu,i)                  ηᵢ
      etadot(mu,i)               ∂ηᵢ/∂τ 
      zet(i,j)                   ζᵢⱼ
      zetdot(i,j)                ∂ζᵢⱼ/∂τ
      dzet(i,j,r)                ∂ᵣζᵢⱼ = ∂ζᵢⱼ/∂αᵣ
      kap(i)                     κᵢ
      kapdot(i)                  ∂κᵢ/∂τ
      dkap(i,j)                  ∂κᵢ/∂αⱼ = ∂ⱼκᵢ
      deta(mu,i,j)               ∂ⱼηᵢ = ∂ηᵢ/∂αⱼ
      Deta(mu,i,j)               Dⱼηᵢ
      [dU/dzet](i,j)             ∂U/∂ζᵢⱼ
      [d2U/dalp.dzet](i)         ∂²U/∂ζⱼᵢ∂αⱼ
      [d2U/dzet.dalp](i)         ∂²U/∂ζᵢⱼ∂αⱼ
      [DL.dU/dzet](i)            Dⱼ∂U/∂ζⱼᵢ
      [DR.dU/dzet](i)            Dⱼ∂U/∂ζᵢⱼ
      [d2U/dzet2](i,j,u,v)       ∂²U/∂ζᵢⱼ∂ζᵤᵥ
    
    #endif
    
    *** DECLARATIONS
    
    OFF statistics;
    
    S c,U,L,T,p;
    I mu,nu,la,i,j,r,s;
    V y,u,a;
    T d1y,d2y,d3y,d1u,d2u,d1a;
    T eta,etadot;
    T [dU/dzet],zetdot(symmetric),Deta,deta;
    CF [DL.dU/dzet],[DR.dU/dzet],[d2U/dalp.dzet],[d2U/dzet.dalp];
    T dkap,dzet;
    CF kap,kapdot,[d2U/dzet2];
    
    *** DEBUGGING
    
    #procedure try
    B u,a,eta,d1y,d1u,d2y;
    print +s;
    .end
    #endprocedure
    
    *** MODULES
    
    #if 0
    * Work out zet derivatives.
    L [zetdot](i,j) = zetdot(i,j); 
    L [dzet](i,j,r) = dzet(i,j,r);
    
    * ∂ζᵢⱼ/∂τ = κᵢ∂κⱼ/∂τ + κⱼ∂κᵢ/∂τ - ∂ᵢu⋅∂�
    ��y - ∂ⱼu⋅∂ᵢy
    Id zetdot(i?,j?) = kap(i)*kapdot(j) + kap(j)*kapdot(i)  
      - d1u(nu,i)*d1y(nu,j) - d1u(nu,j)*d1y(nu,i);
    Sum nu;
    
    * ∂ζᵢⱼ/∂αᵣ = ∂ᵣζᵢⱼ = - ∂ᵣηᵢ⋅ηⱼ - ∂ᵣηⱼ⋅
    ηᵢ
    Id dzet(i?,j?,r?) = - deta(nu,i,r)*eta(nu,j) - deta(nu,j,r)*eta(nu,i);
    Sum nu;
    
    * ∂ᵢηⱼ = ∂ᵢ∂ⱼy - ∂ᵢu u⋅∂ⱼy/c² - u ∂ᵢu⋅∂ⱼy/c�
    � - uu⋅∂ᵢ∂ⱼy/c²
    Id deta(nu?,j?,i?) = d2y(nu,i,j) - d1u(nu,i)*d1y(u,j)/c^2
      - u(nu)*d1u(la,i)*d1y(la,j)/c^2 - u(nu)*d2y(u,i,j)/c^2;
    Sum la;
    
    Id eta(u,i?) = 0;
    Id eta(mu?,i?) = d1y(mu,i) -u(mu)*d1y(u,i)/c^2;
    Id d1u(u,i?) = 0;
    
    #call try
    #endif
    
    #if 0
      Jγ div E = a/c² U
                 + ∂U/∂ζᵢⱼ ( u/c² ∂ζᵢⱼ/∂τ + ηᵢ/c ∂κⱼ/∂
    τ + ηⱼ/c ∂κᵢ/∂τ
                             - Dᵢηⱼ - Dⱼηᵢ )
                 - ηᵢ Dⱼ∂U/∂ζᵢⱼ - ηⱼ Dᵢ∂U/∂ζᵢⱼ
    #endif
    
    * This leaves out the overall factor 1/Jγ.
    G divE(mu) = a(mu)/c^2*U + [dU/dzet](i,j)*
      (
        u(mu)/c^2*zetdot(i,j)
        + eta(mu,i)*kapdot(j)/c + eta(mu,j)*kapdot(i)/c
        - Deta(mu,j,i) - Deta(mu,i,j)
      )
      - eta(mu,i)*[DR.dU/dzet](i) - eta(mu,i)*[DL.dU/dzet](i)
    ;
    Sum i,j;
    
    * Dᵢηⱼ = ∂ᵢ∂ⱼy - ∂ᵢκⱼ u/c - κᵢ∂ⱼu/c + κᵢ∂κⱼ/�
    ��τ u/c² + κᵢκⱼ a/c²
    * Dᵢηⱼ = ∂ᵢ∂ⱼy - ∂ᵢκⱼ u/c - κᵢ∂ⱼu/c - κⱼ∂ᵢu/c 
    + κᵢ∂κⱼ/∂τ u/c²
    *        + κᵢκⱼ a/c² (p. 30)
    Id Deta(mu?,i?,j?) = d2y(mu,i,j) -  dkap(j,i)*u(mu)/c
      - kap(i)*d1u(mu,j)/c - kap(j)*d1u(mu,i)/c
      + kap(i)*kapdot(j)*u(mu)/c^2
      + kap(i)*kap(j)*a(mu)/c^2;
    
    * ∂ζᵢⱼ/∂τ = κᵢ∂κⱼ/∂τ + κⱼ∂κᵢ/∂τ - ∂ᵢu⋅η�
    � - ∂ⱼu⋅ηᵢ
    Id zetdot(i?,j?) = kap(i)*kapdot(j) + kap(j)*kapdot(i)
      - d1u(la,i)*eta(la,j) - d1u(la,j)*eta(la,i);
    Sum la;
    
    * ∂ᵢκⱼ = ∂ᵢu⋅∂ⱼy/c + u⋅∂ᵢ∂ⱼy/c = ∂ᵢu⋅ηⱼ/c + 
    u⋅∂ᵢ∂ⱼy/c
    Id dkap(j?,i?) = d1u(la,i)*eta(la,j)/c + d2y(u,i,j)/c;
    Sum la;
    .sort
    
    * Inspection of the result so far shows that it is manifestly
    * orthogonal to u.  Nevertheless, check the conservation law,
    * dot product should vanish.
    
    L [u.divE1] = u(mu)*divE(mu);
    Id u.u = c^2;
    Id u.a = 0;
    Id eta(u,i?) = 0;
    Id d1u(u,i?) = 0;
    .sort
    
    Hide [u.divE1];
    
    * Work out D derivatives.
    
    * Dᵢ∂U/∂ζᵢⱼ = ∂²U/∂ζᵢⱼ∂ζᵤᵥ ∂ᵢζᵤᵥ - κᵢ/c 
    ∂²U/∂ζᵢⱼ∂ζᵤᵥ ∂ζᵤᵥ/∂τ
    *   + ∂²U/∂ζᵢⱼ∂αᵢ
    Id [DL.dU/dzet](j?) = [d2U/dzet2](i,j,r,s)*dzet(r,s,i)
      - kap(i)/c*[d2U/dzet2](i,j,r,s)*zetdot(r,s) + [d2U/dalp.dzet](j);
    
    * Dⱼ∂U/∂ζᵢⱼ = ∂²U/∂ζᵢⱼ∂ζᵤᵥ ∂ⱼζᵤᵥ - κⱼ/c 
    ∂²U/∂ζᵢⱼ∂ζᵤᵥ ∂ζᵤᵥ/∂τ
    *   + ∂²U/∂ζᵢⱼ∂αⱼ (right div)
    Id [DR.dU/dzet](i?) = [d2U/dzet2](i,j,r,s)*dzet(r,s,j)
      - kap(j)/c*[d2U/dzet2](i,j,r,s)*zetdot(r,s) + [d2U/dzet.dalp](i);
    
    Sum i,j;
    .sort
    
    * Work out zet derivatives.
    
    * ∂ζᵢⱼ/∂τ = κᵢ∂κⱼ/∂τ + κⱼ∂κᵢ/∂τ - ∂ᵢu⋅∂�
    ��y - ∂ⱼu⋅∂ᵢy
    Id zetdot(i?,j?) = kap(i)*kapdot(j) + kap(j)*kapdot(i)  
      - d1u(nu,i)*d1y(nu,j) - d1u(nu,j)*d1y(nu,i);
    Sum s,r,nu;
    
    * ∂ζᵢⱼ/∂αᵣ = ∂ᵣζᵢⱼ = - ∂ᵣηᵢ⋅ηⱼ - ∂ᵣηⱼ⋅
    ηᵢ
    Id dzet(i?,j?,r?) = - deta(nu,i,r)*eta(nu,j) - deta(nu,j,r)*eta(nu,i);
    Sum nu;
    
    * ∂ᵢηⱼ = ∂ᵢ∂ⱼy - ∂ᵢu u⋅∂ⱼy/c² - u ∂ᵢu⋅∂ⱼy/c�
    � - uu⋅∂ᵢ∂ⱼy/c²
    Id deta(nu?,j?,i?) = d2y(nu,i,j) - d1u(nu,i)*d1y(u,j)/c^2
      - u(nu)*d1u(la,i)*d1y(la,j)/c^2 - u(nu)*d2y(u,i,j)/c^2;
    Sum la;
    
    * Eliminate kap and kapdot.
    
    * κᵢ ≡ u⋅∂ᵢy/c
    Id kap(i?) = d1y(u,i)/c;
    * ∂κᵢ/∂τ = a⋅∂ᵢy/c
    Id kapdot(i?) = d1y(a,i)/c;
    Id u.u = c^2;
    Id u.a = 0;
    Id eta(u,i?) = 0;
    Id d1u(u,i?) = 0;
    Id d1a(u,i?) = -d1u(a,i);
    .sort
    
    * Check conservation law, should vanish.
    
    L [u.divE2] = divE(mu)*u(mu);
    Id u.u = c^2;
    Id a(u) = 0;
    Id d1u(u,i?) = 0;
    Id eta(u,i?) = 0;
    .sort
    
    Hide [u.divE2];
    Hide divE;
    
    * Check consistent units, (if no cancellations).
    
    L units = divE(mu);
    Id c^p? = (L/T)^p;
    Id a(mu) = L/T^2;
    Id u(mu) = L/T;
    Id eta(a,i?) = L/T^2;
    Id eta(mu?,i?) = 1;
    Id d1y(u,i?) = L/T;
    Id d1y(a,i?) = L/T^2;
    Id d1y(mu?,i?) = 1;
    Id d1u(mu?,i?) = 1/T;
    Id d2y(u,i?,j?) = 1/T;
    Id d2y(mu?,i?,j?) = 1/L;
    Id [d2U/dalp.dzet](i?) = U/L;
    Id [d2U/dzet.dalp](i?) = U/L;
    Id [dU/dzet](i?,j?) = U;
    .sort
    
    * Write results.
    
    Unhide [u.divE1];
    Unhide [u.divE2];
    Unhide divE;
    B u,a,eta,d1y,d1u,d2y;
    Print +s;
    .store

   divE(mu) =

       + a(mu) * (
          + c^-2*U
          )

       + d1y(u,N1_?)*d1y(u,N2_?)*a(mu) * (
          - 2*[dU/dzet](N1_?,N2_?)*c^-4
          )

       + d1y(u,N1_?)*d1y(u,N2_?)*d1y(a,N3_?)*eta(mu,N4_?) * (
          + [d2U/dzet2](N2_?,N4_?,N1_?,N3_?)*c^-4
          + [d2U/dzet2](N2_?,N4_?,N3_?,N1_?)*c^-4
          + [d2U/dzet2](N4_?,N2_?,N1_?,N3_?)*c^-4
          + [d2U/dzet2](N4_?,N2_?,N3_?,N1_?)*c^-4
          )

       + d1y(u,N1_?)*d1y(N2_?,N3_?)*d1u(N2_?,N4_?)*eta(mu,N5_?) * (
          - [d2U/dzet2](N1_?,N5_?,N3_?,N4_?)*c^-2
          - [d2U/dzet2](N1_?,N5_?,N4_?,N3_?)*c^-2
          - [d2U/dzet2](N5_?,N1_?,N3_?,N4_?)*c^-2
          - [d2U/dzet2](N5_?,N1_?,N4_?,N3_?)*c^-2
          )

       + d1y(u,N1_?)*d1u(mu,N2_?) * (
          + 2*[dU/dzet](N1_?,N2_?)*c^-2
          + 2*[dU/dzet](N2_?,N1_?)*c^-2
          )

       + d1y(u,N1_?)*d1u(N2_?,N3_?)*eta(mu,N4_?)*eta(N2_?,N5_?) * (
          - [d2U/dzet2](N3_?,N4_?,N1_?,N5_?)*c^-2
          - [d2U/dzet2](N3_?,N4_?,N5_?,N1_?)*c^-2
          - [d2U/dzet2](N4_?,N3_?,N1_?,N5_?)*c^-2
          - [d2U/dzet2](N4_?,N3_?,N5_?,N1_?)*c^-2
          )

       + d1y(a,N1_?)*eta(mu,N2_?) * (
          + [dU/dzet](N1_?,N2_?)*c^-2
          + [dU/dzet](N2_?,N1_?)*c^-2
          )

       + d2y(u,N1_?,N2_?)*u(mu) * (
          + [dU/dzet](N1_?,N2_?)*c^-2
          + [dU/dzet](N2_?,N1_?)*c^-2
          )

       + d2y(mu,N1_?,N2_?) * (
          - [dU/dzet](N1_?,N2_?)
          - [dU/dzet](N2_?,N1_?)
          )

       + d2y(N1_?,N2_?,N3_?)*eta(mu,N4_?)*eta(N1_?,N5_?) * (
          + [d2U/dzet2](N2_?,N4_?,N3_?,N5_?)
          + [d2U/dzet2](N2_?,N4_?,N5_?,N3_?)
          + [d2U/dzet2](N4_?,N2_?,N3_?,N5_?)
          + [d2U/dzet2](N4_?,N2_?,N5_?,N3_?)
          )

       + eta(mu,N1_?) * (
          - [d2U/dalp.dzet](N1_?)
          - [d2U/dzet.dalp](N1_?)
          );

   [u.divE1] = 0;

   [u.divE2] = 0;

   units =
       + 3*U*L^-1
         ;

    
    Save elasforce.sav;
    .end
  0.00 sec out of 0.00 sec