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.