unit Physikunit;

{$mode objfpc}{$H+}

{ $DEFINE Zeitschrittueberwachung}
{ $DEFINE Dichteueberwachung}

interface

uses
  Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit;

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;
    name:      string;
    datei:     file;
  end;

  tFelder = class;

  { tWerteGitter }

  tWerteGitter = class(tObject) // repräsentiert ein Gitter von Werten und deren Zeitableitungen
    werte: array[boolean] of array of extended;
    par:   tFelder;
  end;

  { TFelder }

  tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen
  public
    materieFelder:            array of tWerteGitter; // Materiefelder (getrennt für verschiedene Teilchenspezies) und deren Zeitableitungen
    elektroMagnetischeFelder: tWerteGitter;          // EM-Felder und deren Zeitableitungen
      // A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier!
    dX,iDX,x,pDNMax: extended;

    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 tFelder; // mehrere komplette Simulationsboxen von Feldern, benötigt um Zwischenschritte für die Zeitentwicklung zu speichern
    dTMaximum,dX,iDX,xl,xr,t,pDNMax: extended;
    pDNMaxDynamisch: boolean;
    Zeitverfahren: TZeitverfahren;
    Al,Ar: TVerteilungsfunktion;
    constructor create(size: longint; deltaT,deltaX,pDNMa: extended; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant);
    destructor destroy; override;
    procedure iteriereSchritt(var dT: extended);
    procedure macheAusgabe(AD: TAusgabedatei; sDX: extended);
    function gibErhaltungsgroessen: string;
  end;

  { tSimulation }

  tSimulation = class(tObject)
  private
    prot:   tProtokollant;
    gitter: tGitter;
    kvs:    tKnownValues;
  public
    constructor create(inName: string; Protokollant: tProtokollant);
    destructor destroy; override;
  end;

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

{ 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; dichten, lichter: tMyStringlist; 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;
      knownValues.add('x',felders[i,j].x);
      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 length(Felders)-1 do begin
    for j:=0 to length(Felders[i])-1 do
      Felders[i,j].free;
    setlength(Felders[i],0);
  end;
  setlength(Felders,0);
  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;

{ tSimulation }

constructor tSimulation.create(inName: string; Protokollant: tProtokollant);
var
  ifile:                                   tMyStringlist;
  zeitverfahren:                           tZeitverfahren;
  s,t,aSuffix,aPrefix:                     string;
  deltaX,deltaT,endzeit,breite,sDT,pDNMax: extended;
  ausgabeDateien:                          array of tAusgabeDatei;
  abl:                                     boolean;
  fi:                                      tFeldinhalt;
  i,j:                                     longint;
begin
  inherited create;
  prot:=tProtokollant.create(Protokollant,'tSimulation');
  kvs:=tKnownValues.create;

  ifile:=tMyStringlist.create(prot);
  ifile.loadfromfile(inName);
  if not ifile.unfoldMacros then begin
    prot.schreibe('Fehlerhafte Macros in Parameterdatei '''+inName+'''!',true);
    halt(1);
  end;

  // Standardeinstellungen Bereich 'allgemein'
  zeitverfahren:=zfRungeKuttaVier;
  deltaX:=1e-2;
  kvs.add('λ',1/deltaX);
  deltaT:=-1;
  sDT:=-1;
  endzeit:=100;
  breite:=10.0;
  pDNMax:=-1;

  // Standardeinstellungen Bereich 'ausgaben'
  aPrefix:=extractfilepath(paramstr(0));
  aSuffix:='.dat';
  setlength(ausgabeDateien,0);

  while ifile.readln(s) do begin
    if s='allgemein' then begin
      repeat
        if not ifile.readln(s) then begin
          prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich allgemein!',true);
          halt(1);
        end;
        if s='allgemeinEnde' then break;
        if s='runge-Kutta-4' then begin
          Zeitverfahren:=zfRungeKuttaVier;
          continue;
        end;
        if s='euler-Vorwärts' then begin
          Zeitverfahren:=zfEulerVorwaerts;
          continue;
        end;
        if startetMit('ortsschritt',s) then begin
          deltaX:=exprtofloat(false,s,kvs,nil);
          kvs.add('λ',1/deltaX);
          continue;
        end;
        if startetMit('zeitschritt',s) then begin
          deltaT:=exprtofloat(false,s,kvs,nil);
          kvs.add('dT',deltaT);
          continue;
        end;
        if startetMit('diffusionsterm',s) then begin
          pDNMax:=exprtofloat(false,s,kvs,nil);
          continue;
        end;
        if startetMit('zeit',s) then begin
          endzeit:=exprtofloat(false,s,kvs,nil);
          continue;
        end;
        if startetMit('breite',s) then begin
          breite:=exprtofloat(false,s,kvs,nil);
          continue;
        end;
        prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich allgemein!',true);
        halt(1);
      until false;
      continue;
    end;

    if s='ausgaben' then begin
      repeat
        if not ifile.readln(s) then begin
          prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich ausgaben!',true);
          halt(1);
        end;
        if s='ausgabenEnde' then break;
        if startetMit('suffix',s) then begin
          aSuffix:=s;
          continue;
        end;
        if startetMit('prefix',s) then begin
          aPrefix:=s;
          continue;
        end;
        if startetMit('zeitschritt',s) then begin
          sDT:=exprtofloat(false,s,kvs,nil);
          continue;
        end;
        if startetMit('felder',s) then begin
          s:=s+',';
          while s<>'' do begin
            t:=erstesArgument(s,',');
            for abl:=false to true do
              for fI:=low(tFeldInhalt) to high(tFeldInhalt) do
                if not (abl and (fI in [fiGamma,fiP])) then
                  if t=copy('d',1,byte(abl))+optionNamen[fI] then begin
                    setlength(ausgabedateien,length(ausgabedateien)+1);
                    ausgabedateien[length(ausgabedateien)-1].inhalt:=fI;
                    ausgabedateien[length(ausgabedateien)-1].ableitung:=abl;
                    ausgabedateien[length(ausgabedateien)-1].name:=aPrefix+t+aSuffix;
                    assignFile(ausgabedateien[length(ausgabedateien)-1].datei,ausgabedateien[length(ausgabedateien)-1].name);
                  end;
          end;
          continue;
        end;
        prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich ausgaben!',true);
        halt(1);
      until false;
      continue;
    end;

    prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+'''!',true);
    halt(1);
  end;

  if length(ausgabedateien)=0 then begin
    prot.schreibe('Du solltest irgendetwas abspeichern lassen!',true);
    halt(1);
  end;

  if deltaT<0 then
    deltaT:=deltaX/10;
  if sDT<0 then
    sDT:=deltaT;

  for i:=0 to length(ausgabedateien)-1 do begin
    rewrite(ausgabedateien[i].datei,1);
    j:=0;
    BlockWrite(Ausgabedateien[i].Datei,j,sizeof(longint));
    j:=floor(Endzeit/sDT+1);
    BlockWrite(Ausgabedateien[i].Datei,j,sizeof(longint));
  end;
  ifile.free;

  gitter:=tGitter.create(Breite,deltaT,deltaX,pDNMax,@Anfangsdichte,Zeitverfahren,Prot);
end;

destructor tSimulation.destroy;
begin
  gitter.free;
  kvs.free;
  prot.free;
  inherited destroy;
end;

end.