Source-code of the DFEM-algorithm:


Program Oszillator_im_DFEM_mit_OVER_UNITY;
{$APPTYPE CONSOLE}
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs;

Var epo,muo     : Double;  {Constants of nature}
    c           : Double;  {speed of propagation of the waves and fields}
    D           : Double;  {Hooke’s spring force constant}
    m1,m2       : Double;  {Masses of both bodides}
    Q1,Q2       : Double;  {electrical charges of both bodides }
    RLL,FL      : Double;  {relaxed length of the unloaded helical spring}
    r           : Double;  {distance with regard to the finite speed of propagation of the fields}
    diff,ds,ds1 : Double;  {some variables}
    FK1,FK2     : Double;  {spring forces acting on body No.1 and 2}
    FEL1,FEL2   : Double;  {electrical forces acting on body No.1 and 2}
    delt        : Double;  {time-steps for the motion of the bodies and fields}
    x1,x2,v1,v2 : Array [0..200000] of Real48; {time, position, velocity of the bodies}
    t           : Double;  {variable from the propagation-time of the fields}
    a1,a2       : Double;  {acceleration of the bodies}
    i           : Integer; {counter-variable}
    tj,ts,tr    : Extended;{variable for the determination of the field-propagation-duration in part 3}
    ianf,iend   : Integer; {begin and end of the time under analysis}
    Abstd       : Integer; {distance of the data-points being plotted}
    Ukp,UkpAlt  : Double;  {for part 3}
    unten,neu   : Boolean; {for part 3}
    AmplAnf,AmplEnd : Double; {for the determination of the enhancement of amplitude}
    Reib      : Double; {force of friction}
    P         : Double; {machine power}
    Pn        : Double; {for the determination of the average value of the machine power}
    
Procedure Wait;
Var Ki : Char;
begin                        
  Write('<W>'); Read(Ki); Write(Ki);
  If Ki='e' then Halt;       
end;

Procedure Excel_Datenausgabe(Name:String);
Var fout  : Text;    {file to write a results for excel}
    Zahl  : String;
    i,j   : Integer; { counter-variables}
  begin  {data-output for excel:}
  Assign(fout,Name); Rewrite(fout); {open the file}
  For i:=ianf to iend do  {from "plotanf" to "plotend"}
  begin     
    If (i mod Abstd)=0 then
    begin
{     the first argument is the time:}
      Str(i*delT:10:5,Zahl);
      For j:=1 to Length(Zahl) do
      begin   {replace decimal-points by commata}
        If Zahl[j]<>'.' then write(fout,Zahl[j]);
        If Zahl[j]='.' then write(fout,',');
      end;
      Write(fout,chr(9));  {Tabulator for data-separation}
{     The first function is the Position of particle No. 1:}
      Str(x1[i]:10:5,Zahl);
      For j:=1 to Length(Zahl) do
      begin   { replace decimal-points by commata }
        If Zahl[j]<>'.' then write(fout,Zahl[j]);
        If Zahl[j]='.' then write(fout,',');
      end;
      Write(fout,chr(9));  {Tabulator for data-separation }
{     second column: Position of body 2:}
      Str(x2[i]:10:5,Zahl);
      For j:=1 to Length(Zahl) do
      begin   { replace decimal-points by commata }
        If Zahl[j]<>'.' then write(fout,Zahl[j]);
        If Zahl[j]='.' then write(fout,',');
      end;
      Write(fout,chr(9));  {Tabulator for data-separation }
{     third column: velocity of body 1:}
      Str(v1[i]:10:5,Zahl);
      For j:=1 to Length(Zahl) do
      begin   { replace decimal-points by commata }
        If Zahl[j]<>'.' then write(fout,Zahl[j]);
        If Zahl[j]='.' then write(fout,',');
      end;
      Write(fout,chr(9));  {Tabulator for data-separation }
{     fourth column: velocity of body 2:}
      Str(v2[i]:10:5,Zahl);
      For j:=1 to Length(Zahl) do
      begin   { replace decimal-points by commata }
        If Zahl[j]<>'.' then write(fout,Zahl[j]);
        If Zahl[j]='.' then write(fout,',');
      end;
      Writeln(fout,'');    {line-feed for data-separation}
    end;
  end;
  Close(fout);
end;    

Begin {Main program}
{ Initialisation: }
  D:=0; r:=0;         {Avoid Delphi-Messages}
  epo:=8.854187817E-12;{As/Vm}  {Magnetic field-constant }
  muo:=4*pi*1E-7;{Vs/Am}        {elektric field-constant }
  c:=Sqrt(1/muo/epo);{m/s}      {speed of light }
  m1:=1;{kg}                    {mass of body 1}
  m2:=1;{kg}                    {mass of body 2}
  delt:=1E-3;{sec.}             {Equidistant time-steps for the calculation of the motion}
  ianf:=0; iend:=100000;        {number of the first and the last time.-step }
  Abstd:=2;             {to plot every Abstd-th data-point}

  Writeln('Oscillator in DFEM with OVER-UNITY:');
  Writeln('epo=',epo:20,';  muo=',muo:20,';  c=',c:20);
  Writeln('m1,m2=',m1:15,', ',m2:15,';  D=',D:15);
  Writeln;

{ Begin of the Main Program}
{ Part 1 had been preparations for the program-development, not interesting any further}

{ Teil 2: Test -> anharmonic oscillation, with electrical charge or magnet: STATIC !}
  For i:=ianf to iend do
  begin
    x1[i]:=0;    x2[i]:=0;  {assign the positions to zero}
    v1[i]:=0;    v2[i]:=0;  {assign the velocities to zero}
  end;
  i:=0; {t:=i*delT;} {time in steps of delt.}
  Q1:=2.01E-5{C};  Q2:=2.01E-5{C}; {electrical charge of both bodies}
  D:=0.20;{N/m}                    {Hooke’s spring force constant }
  RLL:=6.0;{m}   {length of the spring without force} {rest-position of the bodies: +/-RLL/2}
  x1[0]:=-3.8;   x2[0]:=+3.8;   {starting-positions of the bodies}
  v1[0]:=00.00;  v2[0]:=00.00;  { starting-velocities of the bodies }
{Now we begin the determination of the motion, step-by-step:}
  Repeat
    i:=i+1;
    FL:=x2[i-1]-x1[i-1]; {length of the spring}
    FK1:=(FL-RLL)*D;  {spring-force, positive pulls to the right side, negative to the left}
    FK2:=(RLL-FL)*D;  {spring-force, positive pulls to the right side, negative to the left}
    FEL1:=0;  FEL2:=0;
    If FL<=1E-20 then
    begin
      Writeln;
      Writeln('Exception: Spring too much compressed in Part 2 at step ',i);
      Excel_Datenausgabe('XLS-Nr-02.DAT');
      Writeln('Data have been stored at "XLS-Nr-02.DAT", then termination of algorithm.');
      Wait; Halt;
    end;
    If FL>1E-20 then
    begin
      FEL1:=+Q1*Q2/4/pi/epo/FL/Abs(FL); {electrostatic force between Q1 & Q2}
      FEL2:=-Q1*Q2/4/pi/epo/FL/Abs(FL); {electrostatic force between Q1 & Q2}
    end;
{Check:} If i=1 then Writeln('El.-force:  ',FEL1,' and ',FEL2,' Newton');
{Check:} If i=1 then Writeln('Spring-force: ',FK1, ' and ',FK2,' Newton');
    a1:=(FK1+FEL1)/m1;  a2:=(FK2+FEL2)/m2; {acceleration of the bodies}
    v1[i]:=v1[i-1]+a1*delt; {alteration of the speed of body 1}
    v2[i]:=v2[i-1]+a2*delt; {alteration of the speed of body 2}
    x1[i]:=x1[i-1]+v1[i-1]*delt; {alteration of the position of body 1}
    x2[i]:=x2[i-1]+v2[i-1]*delt; {alteration of the position of body 2}
  Until i=iend;
  Excel_Datenausgabe('XLS-Nr-02.DAT'); {position and speed as a function of time}
  Writeln('Part 2 is ready.');
  
{ Part 3: Test -> Propagation of the fields with finite speed}
  P:=0; Pn:=0; {assign the machine-power to zero}
  For i:=ianf to iend do
  begin
    x1[i]:=0;    x2[i]:=0;  {assign the positions to zero}
    v1[i]:=0;    v2[i]:=0;  {assign the velocities to zero}
  end;
  i:=0;  {counter for the position and velocity}
  c:=1.4; {Sqrt(1/muo/epo);{m/s} {assign the speed of propagation of the fields here}
  Q1:=3E-5{C};  Q2:=3E-5{C};  {electrical charge of the bodies}
  D:=2.7;{N/m}                    { Hooke’s spring force constant }
  RLL:=8.0;{m}   {length of the spring without force} {rest-position of the bodies: +/-RLL/2}
  x1[0]:=-3.0;   x2[0]:=+3.0;   {starting-position of the bodies}
  v1[0]:=00.00;  v2[0]:=00.00;  {starting-velocity of the bodies }
  Ukp:=x2[0]; UkpAlt:=Ukp; unten:=true; neu:=true; {first reversal point}
  Writeln('reversal-point: ',Ukp:12:6,' m ');
  {Now we begin the determination of the motion, step-by-step:}
  Repeat
    i:=i+1;
    FL:=x2[i-1]-x1[i-1]; {length of the spring}
    FK1:=(FL-RLL)*D;  {spring-force, positive pulls to the right side, negative to the left}
    FK2:=(RLL-FL)*D;  {spring-force, positive pulls to the right side, negative to the left}
{   determination of the Field-motion-duration, Field-motion-distance, and Field-strength}
    FEL1:=0; FEL2:=0;
    tj:=i; ts:=i; {i mesures the time}
                  {Start the iteration with natural figures:}
{   Writeln('tj=',tj*delt:9:5,' ts=',ts*delt:9:5,'=>',x2[Round(tj)]-x1[Round(ts)]-c*(tj-ts)*delt:9:5); }
    Repeat
      ts:=ts-1;
      diff:=x2[Round(tj)]-x1[Round(ts)]-c*(tj-ts)*delt;
{     Writeln('tj=',tj*delt:9:5,' ts=',ts*delt:9:5,'=>',diff:9:5); }
    Until ((diff<0)or(ts<=0));
    If diff>=0 then  {before the motion begin at t=0, the bodies have been in rest.}
    begin
      r:=x2[Round(tj)]-x1[0];
{     Writeln('diff>=0; r=',r); }
    end;
    If diff<0 then   {linear interpolation to determine the fraction after the comma}
    begin
{     Writeln('diff<0 ==> tj=',tj,' ts=',ts);
      Write('x2[',Round(tj),']=',x2[Round(tj)]:13:9);
      Write(' und x1[',Round(ts),']=',x1[Round(ts)]:13:9);
      Write(' und x1[',Round(ts+1),']=',x1[Round(ts+1)]:13:9); Writeln; }
      ds:=x2[Round(tj)]-x1[Round(ts)]-c*(tj-ts)*delt;
      ds1:=x2[Round(tj)]-x1[Round(ts+1)]-c*(tj-(ts+1))*delt;
{     Writeln('ds1=',ds1:13:9,' und ds=',ds:13:9); }
      tr:=ts*delt+delt*(-ds)/(ds1-ds); {for linear interpolation}
      tj:=tj*delt;
{     Write('tj=',tj:13:9,' und tr_vor=',tr:13:9); }
      tr:=(tj-tr);  {interpolated moment of field-emission}
      r:=c*tr;           {interpolated real distance}
{     Writeln(' und tr=',tr:13:9,' und r=',r:13:9); } 
    end;  
    If r<=1E-10 then
    begin
      Writeln;
      Writeln('Exception: Spring too much compressed in Part 3 at step ',i);
      Excel_Datenausgabe('XV-03.DAT');
      Writeln('Data have been stored at "XV-03.DAT", then termination of algorithm.');
      Wait; Halt;
    end;
    If r>1E-10 then {Now insert data into Coulomb’s law:}
    begin
      FEL1:=+Q1*Q2/4/pi/epo/r/Abs(r); {electrostatic force between Q1 & Q2}
      FEL2:=-Q1*Q2/4/pi/epo/r/Abs(r); { electrostatic force between Q1 & Q2}
    end;
    Reib:=0.2;    {friction: computation begins here.}
    If i>=10000 then
    begin
      If FEL1>0 then FEL1:=FEL1-Reib;
      If FEL1<0 then FEL1:=FEL1+Reib;
      If FEL2>0 then FEL2:=FEL2-Reib;
      If FEL2<0 then FEL2:=FEL2+Reib;
      P:=P+Reib*Abs(x1[i]-x1[i-1])/delt;
      Pn:=Pn+1;
    end;           {Friction: computation ends here.}
    {Check:} If i=1 then Writeln('El.-force:  ',FEL1,' and ',FEL2,' Newton');
    {Check:} If i=1 then Writeln('spring-force: ',FK1, ' and ',FK2,' Newton');
    a1:=(FK1+FEL1)/m1;  a2:=(FK2+FEL2)/m2; {acceleration of the bodies}
    v1[i]:=v1[i-1]+a1*delt; {alteration of the speed of body 1}
    v2[i]:=v2[i-1]+a2*delt; {alteration of the speed of body 2}
    x1[i]:=x1[i-1]+v1[i-1]*delt; {alteration of the position of body 1}
    x2[i]:=x2[i-1]+v2[i-1]*delt; {alteration of the position of body 2}
{   If (i mod 1000)=0 then Writeln ('Feldstaerke= ',Q1/4/pi/epo/r/Abs(r),' N/C'); }
{   determination of the reversal-points, for determination of the amplitude’s-enhancement:}
    If unten then
    begin                                 
      If x2[i]>Ukp then begin Ukp:=x2[i]; end;
      If x2[i]<Ukp then
      begin
        Writeln('reversal-point: ',Ukp:12:6,' m , amplitude=',Abs(UkpAlt-Ukp));
        If Not(neu) then  AmplEnd:=Abs(UkpAlt-Ukp);
        If neu then begin AmplAnf:=Abs(UkpAlt-Ukp); neu:=false; end;
        unten:=Not(unten); UkpAlt:=Ukp;
      end;
    end;
    If Not(unten) then
    begin
      If x2[i]<Ukp then begin Ukp:=x2[i];  end;
      If x2[i]>Ukp then
      begin
        Writeln('reversal-point: ',Ukp:12:6,' m , amplitude=',Abs(UkpAlt-Ukp));
        If Not(neu) then  AmplEnd:=Abs(UkpAlt-Ukp);
        If neu then begin AmplAnf:=Abs(UkpAlt-Ukp); neu:=false; end;
        unten:=Not(unten); UkpAlt:=Ukp;
      end;
    end;
  Until i=iend;
  Writeln('enhancement of the amplitude: ',AmplEnd-AmplAnf:12:6,' Meter. ');
  Writeln('The machine-power is',P/Pn,' Watt.');
  Excel_Datenausgabe('XV-03.DAT'); {position and speed as a function of time}
  Wait; Wait;
End.