unit Physikunit;

{$mode objfpc}{$H+}

{ $DEFINE Zeitschrittueberwachung}
{ $DEFINE Dichteueberwachung}

interface

uses
  Classes, SysUtils, Math;

type
  TZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaVier);
  TVerteilungsfunktion = function(x: extended): extended;
  TFeldInhalt = (
    fiA,fiAX,fiAY,fiAZ,
    fidAXDT,fidAYDT,fidAZDT,
    fiN,fidPhiDX,fidPsiDX,
    fiP,fiPX,fiPY,fiPZ,
    fiGamma,fiiGamma);

  TAusgabeDatei = record
    Ableitung: Boolean;
    Inhalt:    TFeldInhalt;
    Datei:     File;
  end;

  TProtokolldatei = record
    datei:       textfile;
    einrueckung: longint;
  end;

  { TProtokollant }

  TProtokollant = class(TObject)
  private
    sDat:   ^TProtokolldatei;
    Bes:    String;
    Kinder: array of TProtokollant;
    Elter:  TProtokollant;
  public
    constructor create(Dateiname: string); overload;
    constructor create(Elter_: TProtokollant; Besitzer: string); overload;
    destructor destroy; override;
    procedure schreibe(s: string);
    procedure schreibe(s: string; tee: boolean);
    procedure destroyall;
  end;

  { TFeld }

  pTFeld = ^TFeld;
  TFeld = class(TObject)
  public
    Groessen: array[Boolean,TFeldInhalt] of extended;
      // A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier!
    dX,iDX,x,pDNMax: extended;
    lN,rN:           pTFeld;
    constructor create(deltaX,pDNMa: extended);
    procedure berechneGammaUndP;
    procedure berechneNAbleitung;
    procedure berechneAbleitungen(dT: extended);
  end;

  { TGitter }

  TGitter = class(TObject)
  private
    groesse: longint;
    Prot:    TProtokollant;
    procedure berechneAbleitungen(Felder: longint; dt: extended); overload;
    procedure berechneAbleitungen(Felder: longint; dt: extended; out dTMax: extended); overload;
    procedure setzeRaender(Felder: longint);
  public
    aktuelleFelder: longint;
    Felders: array of array of TFeld;
    dTMaximum,dX,iDX,xl,xr,t,pDNMax: extended;
    pDNMaxDynamisch: boolean;
    Zeitverfahren: TZeitverfahren;
    Al,Ar: TVerteilungsfunktion;
    constructor create(size: longint; deltaT,deltaX,pDNMa: extended; n0: TVerteilungsfunktion; ZV: TZeitverfahren; Protokollant: TProtokollant);
    destructor destroy; override;
    procedure iteriereSchritt(var dT: extended);
    procedure macheAusgabe(AD: TAusgabedatei; sDX: extended);
    function gibErhaltungsgroessen: string;
  end;

function Nullfunktion(x: extended): extended;
function timetostr(t: extended): string;

const
  OptionNamen: array[TFeldInhalt] of string = (
    'A','AX','AY','AZ','dAXDT','dAYDT','dAZDT','N','dPhiDX','dPsiDX','P','PX','PY','PZ','Gamma','iGamma');
  AusgabeNamen: array[TFeldInhalt] of string = (
    'fiA','fiAX','fiAY','fiAZ','fidAXDT','fidAYDT','fidAZDT','fiN','fidPhiDX','fidPsiDX','fiP','fiPX','fiPY','fiPZ','fiGamma','fiiGamma');
  AusgabeHilfe = '--d?A[XYZ]? --d?dA[XYZ]dt --d?N --d?dPhiDX --d?dPsiDX --P[XYZ]? --Gamma --iGamma';

implementation

{ TProtokollant }

constructor TProtokollant.create(Dateiname: string);
begin
  inherited create;
  getmem(sDat,sizeof(TProtokolldatei));
  sDat^.einrueckung:=10;
  assignfile(sDat^.datei,Dateiname);
  rewrite(sDat^.datei);
  Bes:='';
  setlength(Kinder,0);
  Elter:=nil;
end;

constructor TProtokollant.create(Elter_: TProtokollant; Besitzer: string);
begin
  inherited create;
  sDat:=Elter_.sDat;
  setlength(Elter_.Kinder,length(Elter_.Kinder)+1);
  Elter_.Kinder[length(Elter_.Kinder)-1]:=self;
  Elter:=Elter_;
  setlength(Kinder,0);
  Bes:=Elter.Bes+'.'+Besitzer;
  if pos('.',Bes)=1 then delete(Bes,1,1);
  if length(Bes)+4 > sDat^.einrueckung then
    sDat^.einrueckung:=length(Bes)+4;
end;

destructor TProtokollant.destroy;
var i: longint;
begin
  while length(Kinder)>0 do
    Kinder[length(Kinder)-1].free;
  if assigned(Elter) then begin
    flush(sDat^.datei);
    i:=0;
    while (i<length(Elter.Kinder)) and (Elter.Kinder[i]<>self) do
      inc(i);
    if i>=length(Elter.Kinder) then
      schreibe('destroy fehlgeschlagen, ich kann mich nicht sehen.');
    inc(i);
    while i<length(Elter.Kinder) do begin
      Elter.Kinder[i-1]:=Elter.Kinder[i];
      inc(i);
    end;
    setlength(Elter.Kinder,length(Elter.Kinder)-1);
  end
  else begin
    closefile(sDat^.datei);
    freemem(sDat,sizeof(TProtokolldatei));
  end;
  inherited destroy;
end;

procedure TProtokollant.destroyall;
begin
  if Assigned(Elter) then Elter.destroyall
  else destroy;
end;

procedure TProtokollant.schreibe(s: string);
begin
  schreibe(s,false);
end;

procedure TProtokollant.schreibe(s: string; tee: boolean);
var i: longint;
begin
  for i:=length(Bes)+1 to sDat^.einrueckung do
    s:=' '+s;
  writeln(sDat^.datei,Bes+s);
  if tee then
    writeln(Bes+s);
end;

{ TFeld }

constructor TFeld.create(deltaX,pDNMa: extended);
var FI:  TFeldinhalt;
    Abl: Boolean;
begin
  inherited create;
  lN:=nil;
  rN:=nil;
  for Abl:=False to True do
    for FI:=low(TFeldinhalt) to high(TFeldinhalt) do
      Groessen[Abl,FI]:=0;
  dX:=deltaX;
  iDX:=1/dX;
  x:=0;
  pDNMax:=pDNMa;
  Groessen[false,fiN]:=1;
  Groessen[false,fiGamma]:=1;
  Groessen[false,fiiGamma]:=1;
end;

procedure TFeld.berechneGammaUndP; // aus aktuellem dPsi/dX und A
var b:   byte;
    tmp: extended;
begin
  tmp:=0;
  for b:=0 to 2 do begin
    Groessen[false,TFeldinhalt(Byte(fiPX)+b)]:=
      Groessen[false,TFeldinhalt(Byte(fiAX)+b)] + Groessen[false,fidPsiDX]*Byte(b=0);
    tmp:=tmp+sqr(Groessen[false,TFeldinhalt(Byte(fiPX)+b)]);
  end;
  Groessen[false,fiGamma]:=sqrt(1+tmp);
  Groessen[false,fiiGamma]:=1/Groessen[false,fiGamma];
end;

procedure TFeld.berechneNAbleitung; // Zeitableitung von n berechnen
begin
(*  // dn/dt = -d(n/Gamma p_x)/dx
    Groessen[true,fiN]:=
      -(   rN^.rN^.Groessen[false,fiN]*rN^.rN^.Groessen[false,fiiGamma]*rN^.rN^.Groessen[false,fiPX]
         + rN^.Groessen[false,fiN]*rN^.Groessen[false,fiiGamma]*rN^.Groessen[false,fiPX]
         - lN^.Groessen[false,fiN]*lN^.Groessen[false,fiiGamma]*lN^.Groessen[false,fiPX]
         - lN^.lN^.Groessen[false,fiN]*lN^.lN^.Groessen[false,fiiGamma]*lN^.lN^.Groessen[false,fiPX]
       )*iDX/6;
    // es wird über die beiden nächsten linken bzw. rechten Nachbarn gemittelt   *)

  // dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n
  Groessen[true,fiN]:=
    ( rN^.Groessen[false,fiN]*rN^.Groessen[false,fiiGamma]*(pDNMax - rN^.Groessen[false,fiPX])
    - lN^.Groessen[false,fiN]*lN^.Groessen[false,fiiGamma]*(-pDNMax - lN^.Groessen[false,fiPX])
    - Groessen[false,fiN]*Groessen[false,fiiGamma]*2*pDNMax)*iDX/2
  // Der zweite Summand entspricht (in etwa) einer thermischen Verschmierung
  // und soll die numerische Stabilität bis zu p * dn/dx / n von 2*pDNMax gewährleisten.
  // Es handelt sich um einen Diffusionsterm dn/dt = alpha * Laplace n mit
  // alpha = pDNMax * dX
end;

procedure TFeld.berechneAbleitungen(dT: extended); // Zeitableitungen berechnen
var b:  byte;
begin
  // d2Psi/dxdt = dPhi/dx - dGamma/dx
  Groessen[true,fidPsiDX]:=
    Groessen[false,fidPhiDX]
      - (rN^.Groessen[false,fiGamma] - lN^.Groessen[false,fiGamma])*iDX/2;
  // d3Phi/dx2dt = dn/dt  ... hier muss integriert werden:
  // d2Phi/dxdt = d2Phi/dxdt(x- deltax) + <dn/dt> * deltax
  Groessen[true,fidPhiDX]:=
    lN^.Groessen[true,fidPhiDX]
      + (lN^.Groessen[true,fiN] + Groessen[true,fiN])*dX/2;
  // d2A/dt2 = - n/gamma p + Laplace(A) ...
  for b:=0 to 2 do
    Groessen[true,TFeldinhalt(Byte(fidAXDT)+b)]:=
        - (Groessen[false,fiN]*Groessen[false,fiiGamma]*Groessen[false,TFeldinhalt(Byte(fiPX)+b)])
        + (rN^.Groessen[false,TFeldinhalt(Byte(fiAX)+b)]
           - 2*Groessen[false,TFeldinhalt(Byte(fiAX)+b)]
           + lN^.Groessen[false,TFeldinhalt(Byte(fiAX)+b)])*sqr(iDX);
  Groessen[true,fidAXDT]:=Groessen[true,fidAXDT] - Groessen[false,fidPhiDX]; // ... - Nabla(phi)
  // dA/dt = dA/dt
  for b:=0 to 2 do
    Groessen[true,TFeldInhalt(Byte(fiAX)+b)]:=
      Groessen[false,TFeldInhalt(Byte(fidAXDT)+b)] +
      Groessen[true,TFeldInhalt(Byte(fidAXDT)+b)]*dT;
end;

{ TGitter }

procedure TGitter.berechneAbleitungen(Felder: longint; dT: extended; out dTMax: extended);
var i:   longint;
begin
  berechneAbleitungen(Felder,dT);
  dTMax:=dTMaximum;
  for i:=0 to groesse+3 do
    if Felders[Felder,i].Groessen[true,fiN]<0 then begin
      dTMax:=min(dTMax,Abs(Felders[Felder,i].Groessen[false,fiN]/Felders[Felder,i].Groessen[true,fiN]));
    end;
end;

procedure TGitter.berechneAbleitungen(Felder: longint; dT: extended);
var i:         longint;
begin
  for i:=0 to groesse+3 do begin
(*    Felders[Felder,i].Groessen[false,fiN]:=
      max(Felders[Felder,i].Groessen[false,fiN],0); // n >= 0 *)
    Felders[Felder,i].berechneGammaUndP;
  end;
  Felders[Felder,1].Groessen[false,fiPX]:= // Teilchen werden reflektiert
    Abs(Felders[Felder,1].Groessen[false,fiPX]);
  Felders[Felder,groesse+2].Groessen[false,fiPX]:=
    -Abs(Felders[Felder,groesse+2].Groessen[false,fiPX]);
  setzeRaender(Felder);
  for i:=1 to groesse+2 do
    Felders[Felder,i].berechneNAbleitung;
  for i:=1 to groesse+2 do
    Felders[Felder,i].berechneAbleitungen(dT);
end;

procedure TGitter.setzeRaender(Felder: longint);
var FI: TFeldInhalt;
begin
  for FI:=fiAX to fiAZ do begin // Vakuumrandbedingungen für das A-Feld
    Felders[Felder,0].Groessen[true,FI]:=
      (Felders[Felder,1].Groessen[false,FI] -
       Felders[Felder,0].Groessen[false,FI])*iDX;
    Felders[Felder,groesse+3].Groessen[true,FI]:=
      (Felders[Felder,groesse+2].Groessen[false,FI] -
       Felders[Felder,groesse+3].Groessen[false,FI])*iDX;
  end; // (ein bisschen wird noch reflektiert, vmtl. durch numerische Fehler)
  Felders[Felder,0].Groessen[true,fiAy]:=
    Felders[Felder,0].Groessen[true,fiAy]
     + (Al(t+dTMaximum)-Al(t))/dTMaximum;
  Felders[Felder,groesse+3].Groessen[true,fiAy]:=
    Felders[Felder,groesse+3].Groessen[true,fiAy]
     + (Ar(t+dTMaximum)-Ar(t))/dTMaximum;
end;

constructor TGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; n0: TVerteilungsfunktion; ZV: TZeitverfahren; Protokollant: TProtokollant);
var i,j: longint;
begin
  inherited create;
  Ar:=@Nullfunktion;
  Al:=@Nullfunktion;
  Zeitverfahren:=ZV;
  groesse:=size;
  Prot:=TProtokollant.create(Protokollant,'TGitter');
  dTMaximum:=deltaT;
  dX:=deltaX;
  iDX:=1/dX;
  if pDNMa < 0 then begin
    pDNMax:=0;
    pDNMaxDynamisch:=true;
  end
  else pDNMax:=pDNMa;
  case Zeitverfahren of
    zfEulerVorwaerts: Setlength(Felders,2);
    zfRungeKuttaVier: Setlength(Felders,5);
  end{of Case};
  xl:=dX/2;
  for i:=0 to length(Felders)-1 do begin
    setlength(Felders[i],size+4); // auf jeder Seite je zwei Felder Rand extra
    for j:=0 to length(felders[i])-1 do begin
      felders[i,j]:=TFeld.create(dX,pDNMax);
      felders[i,j].x:=xl+j*dX;
      felders[i,j].Groessen[false,fiN]:=n0(felders[i,j].x);
      felders[i,j].Groessen[false,fidPsiDX]:=0.0001*Cos(j/(length(felders[i])-1)*pi*15);
    end;
    felders[i,0].lN:=@felders[i,0];
    for j:=1 to length(felders[i])-1 do begin
      felders[i,j].lN:=@felders[i,j-1];
      felders[i,j-1].rN:=@felders[i,j];
    end;
    felders[i,length(felders[i])-1].rN:=@felders[i,length(felders[i])-1];
  end;
  xr:=xl+dX*(length(Felders[0])-1);
  aktuelleFelder:=0;
  t:=0;
end;

destructor TGitter.destroy;
var i,j: longint;
begin
  for i:=0 to 1 do
    for j:=0 to length(Felders[i])-1 do
      Felders[i,j].free;
  inherited destroy;
end;

procedure TGitter.iteriereSchritt(var dT: extended);
var i:     longint;
    FI:    TFeldInhalt;
    dTMax: extended;
    {$IFDEF Zeitschrittueberwachung}
    Pro:   TProtokollant;
    {$ENDIF}
const zeitObergrenze  = 1/4;
      zeitSollwert    = 1/8;
      zeitUntergrenze = 1/16;
begin
  berechneAbleitungen(aktuelleFelder,dT,dTMax); // y' = y'(t,y(t))
  {$IFDEF Zeitschrittueberwachung}
  Pro:=TProtokollant.create(Prot,'iteriereSchritt');
  {$ENDIF}
  while dT>=dTMax*zeitObergrenze do begin
    dT:=dTMax*zeitSollwert;
    berechneAbleitungen(aktuelleFelder,dT,dTMax); // y' = y'(t,y(t))
    {$IFDEF Zeitschrittueberwachung}
    Pro.schreibe('+ neues dT: '+floattostr(dT)+' (t='+floattostr(t)+')');
    {$ENDIF}
  end;
  {$IFDEF Zeitschrittueberwachung}
  if dT<1e-30 then begin
    for i:=1 to groesse+2 do
      if (Felders[aktuelleFelder,i].Groessen[true,fiN]<0) and
         (Abs(Felders[aktuelleFelder,i].Groessen[false,fiN]/
              Felders[aktuelleFelder,i].Groessen[true,fiN])<1e-15) then
        Prot.schreibe(
          floattostr(t)+' '+
          inttostr(i)+' '+
          floattostr(Felders[aktuelleFelder,i-1].Groessen[false,fidPsiDX]*
                     Felders[aktuelleFelder,i-1].Groessen[false,fiiGamma])+' '+
          floattostr(Felders[aktuelleFelder,i-1].Groessen[false,fiN])+' '+
          floattostr(Felders[aktuelleFelder,i-1].Groessen[true,fiN])+' '+
          floattostr(Felders[aktuelleFelder,i].Groessen[false,fidPsiDX]*
                     Felders[aktuelleFelder,i].Groessen[false,fiiGamma])+' '+
          floattostr(Felders[aktuelleFelder,i].Groessen[false,fiN])+' '+
          floattostr(Felders[aktuelleFelder,i].Groessen[true,fiN])+' '+
          floattostr(Felders[aktuelleFelder,i+1].Groessen[false,fidPsiDX]*
                     Felders[aktuelleFelder,i+1].Groessen[false,fiiGamma])+' '+
          floattostr(Felders[aktuelleFelder,i+1].Groessen[false,fiN])+' '+
          floattostr(Felders[aktuelleFelder,i+1].Groessen[true,fiN]),
          true);
    Pro.destroyall;
    halt(1);
  end;
  {$ENDIF}

  case Zeitverfahren of
    zfEulerVorwaerts:
      begin // y(t+dt) = y(t) + y' dt
        for i:=0 to groesse+3 do
          for FI:=low(TFeldinhalt) to fidPsiDX do
            Felders[1-aktuelleFelder,i].Groessen[false,FI]:=
              Felders[aktuelleFelder,i].Groessen[false,FI]
               + Felders[aktuelleFelder,i].Groessen[true,FI]*dT;
      end;
    zfRungeKuttaVier:
      begin
        for i:=0 to groesse+3 do // ya = y(t) + y' dt/2
          for FI:=low(TFeldinhalt) to fidPsiDX do
            Felders[2,i].Groessen[false,FI]:=
              Felders[aktuelleFelder,i].Groessen[false,FI]
               + Felders[aktuelleFelder,i].Groessen[true,FI]*dT/2;
        berechneAbleitungen(2,dT/2,dTMax);  // ya' = y'(t+dt/2,ya)
        if dT/2>dTMax*zeitObergrenze then begin
          {$IFDEF Zeitschrittueberwachung}
          Pro.schreibe('1. '+floattostr(dT/2)+'>'+floattostr(dTMax*zeitObergrenze)+' (t='+floattostr(t)+')',true);
          Pro.free;
          {$ENDIF}
          dT:=dTMax*zeitSollwert;
          exit;
        end;
        for i:=0 to groesse+3 do // yb = y(t) + ya' dt/2
          for FI:=low(TFeldinhalt) to fidPsiDX do
            Felders[3,i].Groessen[false,FI]:=
              Felders[aktuelleFelder,i].Groessen[false,FI]
               + Felders[2,i].Groessen[true,FI]*dT/2;
        berechneAbleitungen(3,dT/2,dTMax);  // yb' = y'(t+dt/2,yb)
        if dT/2>dTMax*zeitObergrenze then begin
          {$IFDEF Zeitschrittueberwachung}
          Pro.schreibe('2. '+floattostr(dT/2)+'>'+floattostr(dTMax*zeitObergrenze)+' (t='+floattostr(t)+')',true);
          Pro.free;
          {$ENDIF}
          dT:=dTMax*zeitSollwert;
          exit;
        end;
        for i:=0 to groesse+3 do // yc = y(t) + yb' dt
          for FI:=low(TFeldinhalt) to fidPsiDX do
            Felders[4,i].Groessen[false,FI]:=
              Felders[aktuelleFelder,i].Groessen[false,FI]
               + Felders[3,i].Groessen[true,FI]*dT;
        berechneAbleitungen(4,dT/2,dTMax);  // yc' = y'(t+dt,yc)
        if dT/2>dTMax*zeitObergrenze then begin
          {$IFDEF Zeitschrittueberwachung}
          Pro.schreibe('3. '+floattostr(dT/2)+'>'+floattostr(dTMax*zeitObergrenze)+' (t='+floattostr(t)+')',true);
          Pro.free;
          {$ENDIF}
          dT:=dTMax*zeitSollwert;
          exit;
        end;
        for i:=0 to groesse+3 do // y(t+dt) = y(t) + (y' + 2(ya' + yb') + yc') dt/6
          for FI:=low(TFeldinhalt) to fidPsiDX do
            Felders[1-aktuelleFelder,i].Groessen[false,FI]:=
              Felders[aktuelleFelder,i].Groessen[false,FI]
               + (Felders[aktuelleFelder,i].Groessen[true,FI] +
                  2*(Felders[2,i].Groessen[true,FI] +
                     Felders[3,i].Groessen[true,FI]) +
                  Felders[4,i].Groessen[true,FI])*dT/6;
      end;
  end{of case};

  t:=t+dT;
  aktuelleFelder:=1-aktuelleFelder;
  if dT<dTMax*zeitUntergrenze then begin
    dT:=dTMax*zeitSollwert;
    {$IFDEF Zeitschrittueberwachung}
      Pro.schreibe('- neues dT: '+floattostr(dT)+' (t='+floattostr(t)+')');
    {$ENDIF}
  end;
  {$IFDEF Zeitschrittueberwachung}
  for i:=0 to groesse+3 do
    if Felders[aktuelleFelder,i].Groessen[false,fiN]<0 then begin
      Pro.schreibe(
        'n<=0 bei '+
        floattostr(t)+' '+
        inttostr(i)+' '+
        floattostr(Felders[aktuelleFelder,i-1].Groessen[false,fidPsiDX]*
                   Felders[aktuelleFelder,i-1].Groessen[false,fiiGamma])+' '+
        floattostr(Felders[aktuelleFelder,i-1].Groessen[false,fiN])+' '+
        floattostr(Felders[aktuelleFelder,i-1].Groessen[true,fiN])+' '+
        floattostr(Felders[aktuelleFelder,i].Groessen[false,fidPsiDX]*
                   Felders[aktuelleFelder,i].Groessen[false,fiiGamma])+' '+
        floattostr(Felders[aktuelleFelder,i].Groessen[false,fiN])+' '+
        floattostr(Felders[aktuelleFelder,i].Groessen[true,fiN])+' '+
        floattostr(Felders[aktuelleFelder,i+1].Groessen[false,fidPsiDX]*
                   Felders[aktuelleFelder,i+1].Groessen[false,fiiGamma])+' '+
        floattostr(Felders[aktuelleFelder,i+1].Groessen[false,fiN])+' '+
        floattostr(Felders[aktuelleFelder,i+1].Groessen[true,fiN]),true);
      Pro.destroyall;
      halt(1);
    end;
  {$ENDIF}
end;

procedure TGitter.macheAusgabe(AD: TAusgabedatei; sDX: extended);
var i:      longint;
    tmp,sX: extended;
    b:      byte;
begin
  if AD.Inhalt in [fiA,fiP] then
    for i:=0 to length(Felders[aktuelleFelder])-1 do begin
      tmp:=0;
      for b:=1 to 3 do
        tmp:=tmp + sqr(Felders[aktuelleFelder,i].Groessen[AD.Ableitung,TFeldinhalt(Byte(AD.Inhalt)+b)]);
      Felders[aktuelleFelder,i].Groessen[AD.Ableitung,AD.Inhalt]:=sqrt(tmp);
    end;
  i:=floor((xr-xl)/sDX+1);
  tmp:=xl+(i-1)*sDX;
  BlockWrite(AD.Datei,xl,sizeof(extended));
  BlockWrite(AD.Datei,tmp,sizeof(extended));
  BlockWrite(AD.Datei,i,sizeof(longint));
  sX:=xl-Min(dX,sDX)/2;
  for i:=0 to length(Felders[aktuelleFelder])-1 do
    while Felders[aktuelleFelder,i].x>=sX do begin
      BlockWrite(AD.Datei,Felders[aktuelleFelder,i].Groessen[AD.Ableitung,AD.Inhalt],sizeof(extended));
      sX:=sX+sDX;
    end;
end;

function TGitter.gibErhaltungsgroessen: string;
var i,j,k: integer;
    n:     extended;
    Pro:   TProtokollant;
begin
  Pro:=TProtokollant.create(Prot,'gibErhaltungsgroessen');
  n:=0;
  for i:=0 to groesse+3 do
    n:=n+Felders[aktuelleFelder,i].Groessen[false,fiN];
  for i:=1 to groesse+2 do
    if abs(Felders[aktuelleFelder,i].Groessen[false,fiPx]*
            (Felders[aktuelleFelder,i+1].Groessen[false,fiN]-Felders[aktuelleFelder,i-1].Groessen[false,fiN])/
            Max(Felders[aktuelleFelder,i+1].Groessen[false,fiN]+Felders[aktuelleFelder,i-1].Groessen[false,fiN],1e-10))
        > pDNMax then begin
      if pDNMaxDynamisch then begin
        pDNMax:=
          2*
          abs(Felders[aktuelleFelder,i].Groessen[false,fiPx]*
              (Felders[aktuelleFelder,i+1].Groessen[false,fiN]-Felders[aktuelleFelder,i-1].Groessen[false,fiN])/
              Max(Felders[aktuelleFelder,i+1].Groessen[false,fiN]+Felders[aktuelleFelder,i-1].Groessen[false,fiN],1e-10));
        for j:=0 to length(felders)-1 do
          for k:=0 to length(felders[j])-1 do
            felders[j,k].pDNMax:=pDNMax;
        Pro.schreibe('Neues maximales px * dn/dx / n: '+floattostr(pDNMax)+' (t='+floattostr(t)+')');
      end
      else
      if pDNMax>0 then begin
        Pro.schreibe('Warnung: maximaler Impuls * dn/dx / n von '+floattostr(pDNMax)+' wurde überschritten (t='+floattostr(t)+
          'T), die numerische Stabilität ist nicht mehr gewährleistet!',true);
        Pro.schreibe(' Lösung: größeren Diffusionsterm wählen (-D)',true);
        Pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern (-x)',true);
        pDNMax:=-1;
      end;
    end;
  result:='n='+floattostr(n);
  {$IFDEF Dichteueberwachung}
  if n>1000 then begin
    errorCode:=2;
    Pro.schreibe(' n > 1000, es scheinen sehr viele neue Teilchen entstanden zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')');
  end;
  {$ENDIF}
  Pro.free;
end;

function Nullfunktion(x: extended): extended;
begin
  result:=0*x;
end;

function timetostr(t: extended): string;
var b: boolean;
begin
  result:='';
  b:=t>=1;
  if b then begin
    result:=result+inttostr(floor(t))+' Tage ';
    t:=t-floor(t);
  end;
  t:=t*24;
  b:=b or (t>=1);
  if b then begin
    result:=result+inttostr(floor(t))+'h ';
    t:=t-floor(t);
  end;
  t:=t*60;
  b:=b or (t>=1);
  if b then begin
    result:=result+inttostr(floor(t))+'min ';
    t:=t-floor(t);
  end;
  t:=t*60;
  result:=result+inttostr(round(t))+'s';
end;

end.