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;
  tEMFeldInhalt = (
    efA,efAX,efAY,efAZ,
    efDAXDT,efDAYDT,efDAZDT,
    efDPhiDX
  );
  tMaterieFeldInhalt = (
    mfN,mfDPsiDX,
    mfP,mfPX,mfPY,mfPZ,
    mfGamma,mfIGamma
  );

  tGitter = class;

  { tAusgabeDatei }

  tAusgabeDatei = class
  private
    ableitung:          boolean;
    emInhalt:           tEMFeldInhalt;
    matInhalt:          tMaterieFeldInhalt;
    teilchen,nNum,tCnt: longint; // wenn teilchen < 0, dann EM-Feld
    pre,suf:            string;
    datei:              file;
    pro:                tProtokollant;
  public
    zeitAnz:            longint;
    sDT:                extended;
    constructor create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string);
    destructor destroy; override;
    function dump: string;
    procedure schreibeKopf;
    procedure speichereWerte(gitter: tGitter; sDX: extended);
  end;

  { tTeilchenSpezies }

  tTeilchenSpezies = class
  private
    ortsGrenzen:   array of extended;
    dichteStuecke: array of string;
  public
    nMax,                        // maximale Massendichte in nc
    spezifischeLadung: extended; // q/m in e/me
    constructor create;
    destructor destroy; override;
    procedure dump(prot: tProtokollant; prefix: string);
    function gibDichte(x: extended; kvs: tKnownValues): extended; // Massendichte
    procedure liesDichteFunktion(rest: string; ifile: tMyStringList; kvs: tKnownValues; teilchen: array of tTeilchenSpezies; prot: tProtokollant);
  end;

  { tWertePunkt }

  tWertePunkt = class(tObject) // repräsentiert alle Werte eines Punktes im Gitter und deren Zeitableitungen
    matWerte: array of array[tMaterieFeldInhalt,boolean] of double;
      // Materiefelder (getrennt für verschiedene Teilchenspezies) und deren Zeitableitungen
    emWerte:  array[tEMFeldInhalt,boolean] of double;
      // EM-Felder und deren Zeitableitungen
      // A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier!
    lN,rN:    tWertePunkt;

    constructor create(linkerNachbar: tWertePunkt; materien: longint);
    destructor destroy; override;
    procedure berechneGammaUndP;
    procedure berechneNAbleitung(iDX,pDNMax: extended);
    procedure berechneAbleitungen(dX,iDX: extended);
    procedure liKo(in1,in2: tWertePunkt; fak2: extended); overload;                            // Werte werden auf (in1 + fak2*in2') gesetzt
    procedure liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
    function maxDT: extended;
  end;

  { TFelder }

  tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen
  private
    matAnz:   longint;
    lichters: array[boolean] of tMyStringlist;
    procedure setzeRaender(dT,iDX,iDT: extended);
  public
    inhalt: array of tWertePunkt;
    par:    tGitter;

    constructor create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
    destructor destroy; override;
    procedure berechneAbleitungen(dT,dX,iDT,iDX,pDNMax: extended);
    procedure liKo(in1,in2: tFelder; fak2: extended); overload;                            // Werte werden auf (in1 + fak2*in2') gesetzt
    procedure liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
    function maxDT: extended;
  end;

  { tGitter }

  tGitter = class(tObject)
  private
    prot:                         tProtokollant;
    zeitverfahren:                tZeitverfahren;
    pDNMaxDynamisch:              boolean;
    dX,iDX,pDNMax,dTMaximum,xl,t: extended;
    kvs:                          tKnownValues;
    procedure diffusionsTermAnpassen(pro: tProtokollant);
  public

    aktuelleFelder: longint;
    felders: array of tFelder; // mehrere komplette Simulationsboxen von Feldern, benötigt um Zwischenschritte für die Zeitentwicklung zu speichern

    constructor create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string);
    destructor destroy; override;
    procedure iteriereSchritt(var dT: extended);
    function dumpErhaltungsgroessen: boolean;
    procedure berechne(was: char; teilchen: longint);
  end;

  { tSimulation }

  tSimulation = class(tObject)
  private
    prot:                      tProtokollant;
    gitter:                    tGitter;
    dT,tEnde,sT,sDT,sDX:       extended;
    fortschrittsAnzeige:       boolean;
    ausgabeDateien:            array of tAusgabeDatei;
  public
    constructor create(inName: string; protokollant: tProtokollant; name: string);
    destructor destroy; override;
    function iteriereSchritt(start: extended; var zeitPhysik,zeitDatei: extended): boolean; // noch nicht zu Ende?
  end;

implementation

const
  emFeldNamen: array[tEMFeldInhalt] of string = (
    'A','AX','AY','AZ','DAXDT','DAYDT','DAZDT','DPHIDX'
  );
  matFeldNamen: array[tMaterieFeldInhalt] of string = (
    'N','DPSIDX','P','PX','PY','PZ','GAMMA','IGAMMA'
  );

{ tAusgabeDatei }

constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string);
var
  emF:     tEMFeldInhalt;
  maF:     tMaterieFeldInhalt;
  num:     longint;
  abl,gef: boolean;
begin
  inherited create;
  pro:=tProtokollant.create(prot,nam);
  num:=length(feldName);
  while (num>0) and (feldName[num] in ['0'..'9']) do
    dec(num);
  inc(num);
  if num<=length(feldName) then begin
    teilchen:=strtoint(copy(feldName,num,length(feldName)))-1;
    delete(feldName,num,length(feldName));
  end
  else
    teilchen:=0;

  emInhalt:=efA;
  matInhalt:=mfN;
  feldName:=ansiUpperCase(feldName);
  nNum:=0; // Header 0 wurde also noch nicht geschrieben
  zeitAnz:=-1;
  tCnt:=-1;

  gef:=false;
  for abl:=false to true do begin
    if not gef then
      for emF:=low(tEMFeldInhalt) to high(tEMFeldInhalt) do
        if copy('D',1,byte(abl))+emFeldNamen[emF] = feldName then begin
          emInhalt:=emF;
          teilchen:=-1;
          gef:=true;
          break;
        end;
    if not gef then
      for maF:=low(tMaterieFeldInhalt) to high(tMaterieFeldInhalt) do
        if copy('D',1,byte(abl))+matFeldNamen[maF] = feldName then begin
          matInhalt:=maF;
          gef:=true;
          break;
        end;
  end;

  if not gef then begin
    pro.schreibe('tAusgabeDatei.create kennt Feld '''+feldName+''' nicht!',true);
    pro.destroyAll;
    halt(1);
  end;

  if teilchen>=0 then
    pre:=matFeldNamen[maF]+inttostr(teilchen+1)
  else
    pre:=emFeldNamen[emF];

  if ableitung then
    pre:='d'+pre;

  pre:=prefix+pre;
  suf:=suffix;
end;

destructor tAusgabeDatei.destroy;
begin
  if (tCnt<>zeitAnz) and ((nNum<>0) or (tCnt<>-1)) then begin
    pro.schreibe('Falsche Anzahl an Zeitschritten in '+pre+'-'+inttostr(nNum-1)+suf+' geschrieben ('+inttostr(tCnt)+' statt '+inttostr(zeitAnz)+')!',true);
    pro.destroyall;
    halt(1);
  end;
  inherited destroy;
end;

function tAusgabeDatei.dump: string;
begin
  if teilchen<0 then
    result:=emFeldNamen[emInhalt]
  else
    result:=matFeldNamen[matInhalt]+inttostr(teilchen);
  if ableitung then
    result:=result+'''';
  result:=result+' -> '+pre+'-$i'+suf;
end;

procedure tAusgabeDatei.schreibeKopf;
begin
  if (tCnt<>zeitAnz) and ((nNum<>0) or (tCnt<>-1)) then begin
    pro.schreibe('Falsche Anzahl an Zeitschritten in '+pre+'-'+inttostr(nNum-1)+suf+' geschrieben ('+inttostr(tCnt)+' statt '+inttostr(zeitAnz)+')!',true);
    pro.destroyall;
    halt(1);
  end;
  tCnt:=0;
  assignfile(datei,pre+'-'+inttostr(nNum)+suf);
  inc(nNum);
  rewrite(datei,1);
  blockwrite(datei,nNum,sizeof(integer));
  blockwrite(datei,zeitAnz,sizeof(integer));
  closefile(datei);
end;

procedure tAusgabeDatei.speichereWerte(gitter: tGitter; sDX: extended);
var i,cnt: longint;
    sX,cX: extended;
begin
  if (teilchen>=gitter.felders[gitter.aktuelleFelder].matAnz) then begin
    pro.schreibe('Teilchen '+inttostr(teilchen+1)+' gibt es nicht, da kann ich auch nichts speichern!',true);
    pro.destroyall;
    halt(1);
  end;
  if (teilchen<0) and (emInhalt=efA) then
    gitter.berechne('A',teilchen);
  if (teilchen>=0) and (matInhalt=mfP) then
    gitter.berechne('P',teilchen);

  if gitter.t>=nNum+sDT/2 then
    schreibeKopf;

  cnt:=floor((length(gitter.felders[gitter.aktuelleFelder].inhalt)-1)/sDX*gitter.dX+1);
  cX:=gitter.xl;
  sX:=cX+(cnt-1)*sDX;

  inc(tCnt);
  reset(datei,1);
  seek(datei,fileSize(datei));
  blockWrite(datei,cX,sizeof(extended));
  blockWrite(datei,sX,sizeof(extended));
  blockWrite(datei,cnt,sizeof(integer));
  sX:=cX-Min(gitter.dX,sDX)/2;
  if teilchen<0 then begin
    for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
      while cX>=sX do begin
        dec(cnt);
        blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].emWerte[emInhalt,ableitung],sizeof(double));
        sX:=sX+sDX;
      end;
      cX:=cX+gitter.dX;
    end;
  end
  else begin
    for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
      while cX>=sX do begin
        dec(cnt);
        blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].matWerte[teilchen,matInhalt,ableitung],sizeof(double));
        sX:=sX+sDX;
      end;
      cX:=cX+gitter.dX;
    end;
  end;
  closeFile(datei);
  if cnt<>0 then begin
    pro.schreibe('Falsche Anzahl an Ortsschritten geschrieben ('+inttostr(cnt)+')!',true);
    pro.destroyall;
    halt(1);
  end;
end;

{ tTeilchenSpezies }

constructor tTeilchenSpezies.create;
begin
  inherited create;
  fillchar(ortsGrenzen,sizeof(ortsGrenzen),#0);
  setlength(ortsGrenzen,0);
  fillchar(dichteStuecke,sizeof(dichteStuecke),#0);
  setlength(dichteStuecke,0);
  nMax:=0;
  spezifischeLadung:=0;
end;

destructor tTeilchenSpezies.destroy;
begin
  setlength(ortsGrenzen,0);
  setlength(dichteStuecke,0);
  inherited destroy;
end;

procedure tTeilchenSpezies.dump(prot: tProtokollant; prefix: string);
var
  i: longint;
begin
  prot.schreibe(prefix+'nMax = '+floattostr(nMax)+' * nc');
  prot.schreibe(prefix+'q/m  = '+floattostr(spezifischeLadung)+' e/me');
  prot.schreibe(prefix+'n(x):');
  for i:=0 to length(dichteStuecke)-1 do begin
    prot.schreibe(prefix+' '+dichteStuecke[i]);
    if i<length(ortsGrenzen) then
      prot.schreibe(prefix+'  '+floattostr(ortsGrenzen[i]));
  end;
end;

function tTeilchenSpezies.gibDichte(x: extended; kvs: tKnownValues): extended;
var
  i: longint;
begin
  i:=0;
  while (i<length(ortsGrenzen)) and (ortsGrenzen[i]<x) do
    inc(i);
  kvs.add('x',x);
  result:=nMax*exprToFloat(false,dichteStuecke[i],kvs,nil);
  kvs.rem('x');
end;

procedure tTeilchenSpezies.liesDichteFunktion(rest: string; ifile: tMyStringList; kvs: tKnownValues; teilchen: array of tTeilchenSpezies; prot: tProtokollant);
var
  s:     string;
  ori,i: longint;
begin
  if rest='stückweise' then begin
    setlength(ortsGrenzen,0);
    setlength(dichteStuecke,0);
    repeat
      if not ifile.readln(s) then begin
        prot.schreibe('Unerwartetes Dateiende in Parameterdatei im Bereich teilchen -> verteilung stückweise!',true);
        prot.destroyall;
        halt(1);
      end;

      if length(dichteStuecke)=length(ortsGrenzen) then begin // neues Funktionsstück wird definiert
        setlength(dichteStuecke,length(dichteStuecke)+1);
        dichteStuecke[length(dichteStuecke)-1]:=s;
        continue;
      end;

      // kein neues Funktionsstück =>
      if s='stückweiseEnde' then break; // Ende oder
      // neue Ortsgrenze

      setlength(ortsGrenzen,length(ortsGrenzen)+1);
      ortsGrenzen[length(ortsGrenzen)-1]:=exprToFloat(false,s,kvs,nil);
    until false;
  end
  else if startetMit('wie teilchen',rest) then begin
    ori:=strtoint(rest)-1;
    if (ori<0) or (ori>=length(teilchen)-1) then begin
      prot.schreibe('Kann mich nicht auf die Dichte von Teilchen '+inttostr(ori+1)+' beziehen, weil es noch nicht definiert wurde!',true);
      prot.destroyall;
      halt(1);
    end;
    setlength(ortsGrenzen,length(teilchen[ori].ortsGrenzen));
    for i:=0 to length(ortsGrenzen)-1 do
      ortsGrenzen[i]:=teilchen[ori].ortsGrenzen[i];
    setlength(dichteStuecke,length(teilchen[ori].dichteStuecke));
    for i:=0 to length(dichteStuecke)-1 do
      dichteStuecke[i]:=teilchen[ori].dichteStuecke[i];
  end
  else begin
    setlength(ortsGrenzen,0);
    setlength(dichteStuecke,1);
    dichteStuecke[0]:=rest;
  end;
end;

{ tWertePunkt }

constructor tWertePunkt.create(linkerNachbar: tWertePunkt; materien: longint);
var
  emF: tEMFeldInhalt;
  maF: tMaterieFeldInhalt;
  abl: boolean;
  i:   longint;
begin
  inherited create;
  fillchar(matWerte,sizeof(matWerte),#0);
  setlength(matWerte,materien);
  for i:=0 to length(matWerte)-1 do begin
    for maF:=low(tMaterieFeldInhalt) to high(tMaterieFeldInhalt) do
      for abl:=false to true do
        matWerte[i,maF,abl]:=0;
    matWerte[i,mfGamma,false]:=0;
    matWerte[i,mfIGamma,false]:=0;
  end;
  for emF:=low(tEMFeldInhalt) to high(tEMFeldInhalt) do
    for abl:=false to true do
      emWerte[emF,abl]:=0;
  lN:=linkerNachbar;
  rN:=nil;
  if assigned(lN) then
    lN.rN:=self;
end;

destructor tWertePunkt.destroy;
begin
  setlength(matWerte,0);
  if assigned(lN) then
    lN.rN:=nil;
  if assigned(rN) then
    rN.lN:=nil;
  inherited destroy;
end;

procedure tWertePunkt.berechneGammaUndP; // aus aktuellem dPsi/dX und A
var
  i: longint;
begin
  for i:=0 to length(matWerte)-1 do begin
    matWerte[i,mfPX,false]:= emWerte[efAX,false] + matWerte[i,mfDPsiDX,false];
    matWerte[i,mfPY,false]:= emWerte[efAY,false];
    matWerte[i,mfPZ,false]:= emWerte[efAZ,false];

    matWerte[i,mfGamma,false]:=
      sqrt(
        1 +
        sqr(matWerte[i,mfPX,false]) +
        sqr(matWerte[i,mfPY,false]) +
        sqr(matWerte[i,mfPZ,false])
      );
    matWerte[i,mfIGamma,false]:=
      1/matWerte[i,mfGamma,false];
  end;
end;

procedure tWertePunkt.berechneNAbleitung(iDX,pDNMax: extended); // Zeitableitung von n berechnen
var
  i: longint;
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   *)

  for i:=0 to length(matWerte)-1 do
    // dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n
    matWerte[i,mfN,true]:=
      ( rN.matWerte[i,mfN,false] * rN.matWerte[i,mfIGamma,false] *  (pDNMax - rN.matWerte[i,mfPX,false])
      - ln.matWerte[i,mfN,false] * ln.matWerte[i,mfIGamma,false] * (-pDNMax - rN.matWerte[i,mfPX,false])
      - matWerte[i,mfN,false]  *  matWerte[i,mfIGamma,false] * 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 tWertePunkt.berechneAbleitungen(dX,iDX: extended); // Zeitableitungen berechnen
var
  i: longint;
begin
  // d2Psi/dxdt = dPhi/dx - dGamma/dx
  for i:=0 to length(matWerte)-1 do
    matWerte[i,mfDPsiDX,true]:=
      emWerte[efDPhiDX,false]
        - (rN.matWerte[i,mfGamma,false] - lN.matWerte[i,mfGamma,false]) * iDX/2;

  // d3Phi/dx2dt = dn/dt  ... hier muss integriert werden:
  // d2Phi/dxdt = d2Phi/dxdt(x- deltax) + <dn/dt> * deltax
  emWerte[efDPhiDX,true]:=
    lN.emWerte[efDPhiDX,true];
  for i:=0 to length(matWerte)-1 do
    emWerte[efDPhiDX,true]:=
      emWerte[efDPhiDX,true]
      + (lN.matWerte[i,mfN,true] + matWerte[i,mfN,true])*dX/2;

  // d2A/dt2 = Laplace(A) - Nabla(phi) ...
  emWerte[efDAXDT,true]:=
    ( rN.emWerte[efAX,false] - 2*emWerte[efAX,false] + lN.emWerte[efAX,false] )*sqr(iDX)
     - emWerte[efDPhiDX,false];
  emWerte[efDAYDT,true]:=
    ( rN.emWerte[efAY,false] - 2*emWerte[efAY,false] + lN.emWerte[efAY,false] )*sqr(iDX);
  emWerte[efDAZDT,true]:=
    ( rN.emWerte[efAZ,false] - 2*emWerte[efAZ,false] + lN.emWerte[efAZ,false] )*sqr(iDX);
  // ... - n/gamma p
  for i:=0 to length(matWerte)-1 do begin
    emWerte[efDAXDT,true]:=
      emWerte[efDAXDT,true]
      - (matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPX,false]);
    emWerte[efDAYDT,true]:=
      emWerte[efDAYDT,true]
      - (matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPY,false]);
    emWerte[efDAZDT,true]:=
      emWerte[efDAZDT,true]
      - (matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPZ,false]);
  end;

  // dA/dt = dA/dt
  emWerte[efAX,true]:=emWerte[efDAXDT,false];// + emWerte[efDAXDT,true]*dT;
  emWerte[efAY,true]:=emWerte[efDAYDT,false];// + emWerte[efDAYDT,true]*dT;
  emWerte[efAZ,true]:=emWerte[efDAZDT,false];// + emWerte[efDAZDT,true]*dT;
end;

procedure tWertePunkt.liKo(in1,in2: tWertePunkt; fak2: extended); // Werte werden auf (in1 + fak2*in2') gesetzt
var
  emF: tEMFeldInhalt;
  maF: tMaterieFeldInhalt;
  i:   longint;
begin
(*  tEMFeldInhalt = (
      efA,efAX,efAY,efAZ,
      efDAXDT,efDAYDT,efDAZDT,
      efDPhiDX
    ); *)
  for emF:=efAX to efDPhiDX do // alles außer efA, welchen Ableitung ja nicht berechnet wurde
    emWerte[emF,false]:= in1.emWerte[emF,false] + fak2 * in2.emWerte[emF,true];

(*  tMaterieFeldInhalt = (
      mfN,mfDPsiDX,
      mfP,mfPX,mfPY,mfPZ,
      mfGamma,mfIGamma
    ); *)
  for i:=0 to length(matWerte)-1 do // siehe oben
    for maF:=mfN to mfDPsiDX do
      matWerte[i,maF,false]:= in1.matWerte[i,maF,false] + fak2 * in2.matWerte[i,maF,true];
end;

procedure tWertePunkt.liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
var
  emF: tEMFeldInhalt;
  maF: tMaterieFeldInhalt;
  i:   longint;
begin
(*  tEMFeldInhalt = (
      efA,efAX,efAY,efAZ,
      efDAXDT,efDAYDT,efDAZDT,
      efDPhiDX
    ); *)
  for emF:=efAX to efDPhiDX do // alles außer efA, welchen Ableitung ja nicht berechnet wurde
    emWerte[emF,false]:=
      in1.emWerte[emF,false]
      + fak2 * in2.emWerte[emF,true]
      + fak3 * in3.emWerte[emF,true]
      + fak4 * in4.emWerte[emF,true]
      + fak5 * in5.emWerte[emF,true];

(*  tMaterieFeldInhalt = (
      mfN,mfDPsiDX,
      mfP,mfPX,mfPY,mfPZ,
      mfGamma,mfIGamma
    ); *)
  for i:=0 to length(matWerte)-1 do // siehe oben
    for maF:=mfN to mfDPsiDX do
      matWerte[i,maF,false]:=
        in1.matWerte[i,maF,false]
        + fak2 * in2.matWerte[i,maF,true]
        + fak3 * in3.matWerte[i,maF,true]
        + fak4 * in4.matWerte[i,maF,true]
        + fak5 * in5.matWerte[i,maF,true];
end;

function tWertePunkt.maxDT: extended;
var
  i: longint;
begin
  result:=1;
  for i:=0 to length(matWerte)-1 do
    if (matWerte[i,mfN,true]<0) and
       (matWerte[i,mfN,false]>0) then
      result:=min(result,matWerte[i,mfN,false]/matWerte[i,mfN,true]);
end;

{ tFelder }

constructor tFelder.create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
var
  i,j:    longint;
  rechts: boolean;
begin
  inherited create;
  par:=parent;
  matAnz:=length(teilchen);
  fillchar(inhalt,sizeof(inhalt),#0);
  setlength(inhalt,groesse+4); // zwei Felder links und rechts extra für Randbedingungen
  inhalt[0]:=tWertePunkt.create(nil,matAnz);
  for i:=1 to length(inhalt)-1 do
    inhalt[i]:=tWertePunkt.create(inhalt[i-1],matAnz);
  for i:=0 to length(inhalt)-1 do
    for j:=0 to length(teilchen)-1 do begin
      inhalt[i].matWerte[j,mfN,false]:=teilchen[j].gibDichte(par.xl+(i-2)*par.dX,parent.kvs);
      inhalt[i].matWerte[j,mfDPsiDX,false]:=0;
    end;
  for rechts:=false to true do begin
    lichters[rechts]:=tMyStringlist.create(nil,'');
    lichters[rechts].text:=lichter.text;
  end;
  lichters[false].grep('^links .*');
  lichters[false].subst('^links *','');
  lichters[true].grep('^rechts .*');
  lichters[true].subst('^rechts *','');
end;

destructor tFelder.destroy;
begin
  lichters[false].free;
  lichters[true].free;
  inherited destroy;
end;

procedure tFelder.setzeRaender(dT,iDX,iDT: extended);
var
  emF:    tEMFeldInhalt;
  rechts: boolean;
  i:      longint;
  oVal:   extended;
begin
  for emF:=efAX to efAZ do begin // Vakuumrandbedingungen für das A-Feld
    inhalt[0].emWerte[emF,true]:=
      (inhalt[1].emWerte[emF,false] -
       inhalt[0].emWerte[emF,false])*iDX;
    inhalt[length(inhalt)-1].emWerte[emF,true]:=
      (inhalt[length(inhalt)-2].emWerte[emF,false] -
       inhalt[length(inhalt)-1].emWerte[emF,false])*iDX;
  end; // (ein bisschen wird noch reflektiert, vmtl. durch numerische Fehler bzw. nicht beachtete numerische Dispersion)

  for rechts:=false to true do
    for i:=0 to lichters[rechts].count-1 do begin
      par.kvs.add('t',par.t-dT);
      oVal:=exprToFloat(false,lichters[rechts][i],par.kvs,nil);
      par.kvs.add('t',par.t);
      inhalt[(length(inhalt)-1)*byte(rechts)].emWerte[efAY,true]:=
        inhalt[(length(inhalt)-1)*byte(rechts)].emWerte[efAY,true] +
        (exprToFloat(false,lichters[rechts][i],par.kvs,nil) - oVal)*iDT;
    end;
end;

procedure tFelder.berechneAbleitungen(dT,dX,iDT,iDX,pDNMax: extended);
var
  i: longint;
begin
  for i:=0 to length(inhalt)-1 do
    inhalt[i].berechneGammaUndP;
  for i:=0 to matAnz-1 do begin // Teilchen werden reflektiert
    inhalt[1].matWerte[i,mfPX,false]:=
       abs(inhalt[1].matWerte[i,mfPX,false]);
    inhalt[length(inhalt)-2].matWerte[i,mfPX,false]:=
      -abs(inhalt[length(inhalt)-2].matWerte[i,mfPX,false]);
  end;

  setzeRaender(dT,iDX,iDT);

  for i:=1 to length(inhalt)-2 do
    inhalt[i].berechneNAbleitung(iDX,pDNMax);
  for i:=1 to length(inhalt)-2 do
    inhalt[i].berechneAbleitungen(dX,iDX);
end;

procedure tFelder.liKo(in1,in2: tFelder; fak2: extended); // Werte werden auf (in1 + fak2*in2') gesetzt
var
  i: longint;
begin
  for i:=0 to length(inhalt)-1 do
    inhalt[i].liKo(in1.inhalt[i],in2.inhalt[i],fak2);
end;

procedure tFelder.liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
var
  i: longint;
begin
  for i:=0 to length(inhalt)-1 do
    inhalt[i].liKo(in1.inhalt[i],in2.inhalt[i],in3.inhalt[i],in4.inhalt[i],in5.inhalt[i],fak2,fak3,fak4,fak5);
end;

function tFelder.maxDT: extended;
var
  i: longint;
begin
  result:=inhalt[0].maxDT;
  for i:=1 to length(inhalt)-1 do
    result:=min(result,inhalt[i].maxDT);
end;

{ tGitter }

constructor tGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string);
var
  i: longint;
begin
  inherited create;
  zeitverfahren:=ZV;
  kvs:=bekannteWerte;
  prot:=tProtokollant.create(protokollant,name);
  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
    felders[i]:=tFelder.create(size,teilchen,lichter,self);

  aktuelleFelder:=0;
  t:=0;
  kvs.add('t',t);
end;

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

procedure tGitter.iteriereSchritt(var dT: extended);
var i:       longint;
    {$IFDEF Zeitschrittueberwachung}
    pro:     tProtokollant;
    j:       longint;
    {$ENDIF}
    iDT,mDT: extended;
begin
  kvs.add('t',t+dT);
  diffusionsTermAnpassen(prot);

  repeat
    iDT:=1/dT;
    felders[aktuelleFelder].berechneAbleitungen(dT,dX,iDT,iDX,pDNMax); // y' = y'(t,y(t))

    mDT:=felders[aktuelleFelder].maxDT;
      // das maximale dT, welches bei den Momentanen Ableitungen nicht zu
      // unphysikalischen Effekten (z.B. negativen Dichten) führt

    if mDT<dT*2 then begin // beabsichtigter Zeitschritt ist zu groß
      dT:=mDT/4;           // machen wir ihn kleiner!
      {$IFDEF Zeitschrittueberwachung}
      pro:=tProtokollant.create(prot,'iteriereSchritt');
      pro.schreibe('dT -> '+floattostr(dT));
      if dT<1E-30 then begin
        pro.schreibe('Zeitschritt geht gegen Null (ist bereits '+floattostr(dT)+' < 10^-30) - irgendwas scheint grundsätzlich kaputt zu sein!',true);
        pro.destroyall;
        halt(1);
      end;
      pro.free;
      {$ENDIF}
      continue;
    end;

    break;
  until false;

  if dT<mDT/1000 then
    dT:=min(1,mDT)/1000;

  case zeitverfahren of
    zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt
      felders[1-aktuelleFelder].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT);
    zfRungeKuttaVier:
      begin
        felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT/2); // ya  = y(t) + y' dt/2
        felders[2].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax);                // ya' = y'(t+dt/2,ya)

        felders[3].liKo(felders[aktuelleFelder],felders[2],dT/2);              // yb  = y(t) + ya' dt/2
        felders[3].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax);                // yb' = y'(t+dt/2,yb)

        felders[4].liKo(felders[aktuelleFelder],felders[3],dT);                // yc  = y(t) + yb' dt
        felders[4].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax);                // yc' = y'(t+dt,yc)

        felders[1-aktuelleFelder].liKo(
          felders[aktuelleFelder],
          felders[aktuelleFelder],
          felders[2],
          felders[3],
          felders[4],
          dT/6,
          dT/3,
          dT/3,
          dT/6
        ); // y(t+dt) = y(t) + (y' + 2(ya' + yb') + yc') dt/6
      end;
  end{of case};

  aktuelleFelder:=1-aktuelleFelder;
  {$IFDEF Dichteueberwachung}
  with felders[aktuelleFelder] do
    for i:=0 to length(inhalt)-1 do
      for j:=0 to matAnz-1 do
        if inhalt[i].matWerte[j,mfN,false]<0 then begin
          pro:=tProtokollant.create(prot,'iteriereSchritt');
          pro.schreibe(
            'n<=0 bei '+
            floattostr(t)+' '+
            inttostr(i)+' '+
            floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]*
                       inhalt[i-1].matWerte[j,mfIGamma,false])+' '+
            floattostr(inhalt[i-1].matWerte[j,mfN,false])+' '+
            floattostr(inhalt[i-1].matWerte[j,mfN,true])+' '+
            floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]*
                       inhalt[i].matWerte[j,mfIGamma,false])+' '+
            floattostr(inhalt[i].matWerte[j,mfN,false])+' '+
            floattostr(inhalt[i].matWerte[j,mfN,true])+' '+
            floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]*
                       inhalt[i+1].matWerte[j,mfIGamma,false])+' '+
            floattostr(inhalt[i+1].matWerte[j,mfN,false])+' '+
            floattostr(inhalt[i+1].matWerte[j,mfN,true]),true);
          pro.destroyall;
          halt(1);
        end;
  {$ENDIF}
  t:=t+dT;
end;

function tGitter.dumpErhaltungsgroessen: boolean;
var i,j: integer;
    ns:  tExtendedArray;
    pro: tProtokollant;
    s:   string;
begin
  pro:=tProtokollant.create(prot,'dumpErhaltungsgroessen');
  setlength(ns,felders[aktuelleFelder].matAnz);

  result:=true;

  s:='';
  for i:=0 to length(ns)-1 do
    ns[i]:=0;
  for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do
    for j:=0 to length(ns)-1 do
      ns[j]:=ns[j] + felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false];
  for i:=0 to length(ns)-1 do begin
    ns[i]:=ns[i]*dX;
    s:=s+ ' n['+inttostr(i)+']='+floattostr(ns[i]);
    {$IFDEF Dichteueberwachung}
    if ns[i]>1000 then begin
      result:=false;
      pro.schreibe(' n['+inttostr(i)+'] > 1000, es scheinen sehr viele neue Teilchen entstanden zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')',true);
    end;
    {$ENDIF}
  end;
  delete(s,1,1);
  prot.schreibe(s);
  pro.free;
end;

procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant);
var
  i,j: longint;
begin
  for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do
    for j:=0 to felders[aktuelleFelder].matAnz-1 do
      if abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPx,false]*
              (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
              max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10))
          > pDNMax then begin
        if pDNMaxDynamisch then begin
          pDNMax:=
            2*
            abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPX,false]*
                (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
                Max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10));
          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;
end;

procedure tGitter.berechne(was: char; teilchen: longint);
var
  i:   longint;
  tmp: extended;
  emF: tEMFeldInhalt;
  maF: tMaterieFeldInhalt;
begin
  if was='A' then begin
    for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do begin
      tmp:=0;
      for emF:=efAX to efAZ do
        tmp:=tmp+sqr(felders[aktuelleFelder].inhalt[i].emWerte[emF,false]);
      felders[aktuelleFelder].inhalt[i].emWerte[efA,false]:=sqrt(tmp);
    end;
    exit;
  end;
  if was='P' then begin
    for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do begin
      tmp:=0;
      for maF:=mfPX to mfPZ do
        tmp:=tmp+sqr(felders[aktuelleFelder].inhalt[i].matWerte[teilchen,maF,false]);
      felders[aktuelleFelder].inhalt[i].matWerte[teilchen,maF,false]:=sqrt(tmp);
    end;
    exit;
  end;
  prot.schreibe('Kann '''+was+''' nicht berechnen, weil ich es nicht verstehe!',true);
  prot.destroyall;
  halt;
end;

{ tSimulation }

constructor tSimulation.create(inName: string; protokollant: tProtokollant; name: string);
var
  ifile:                                   tMyStringlist;
  zeitverfahren:                           tZeitverfahren;
  s,t,aSuffix,aPrefix:                     string;
  deltaX,breite,pDNMax:                    extended;
  i:                                       longint;
  kvs:                                     tKnownValues;
  teilchen:                                array of tTeilchenSpezies;
  lichter:                                 tMyStringlist;
  pro:                                     tProtokollant;
begin
  inherited create;
  prot:=tProtokollant.create(protokollant,name);
  kvs:=tKnownValues.create;

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

  // Standardeinstellungen Bereich 'allgemein'
  zeitverfahren:=zfRungeKuttaVier;
  deltaX:=1e-2;
  kvs.add('λ',1);
  kvs.add('T',1);
  kvs.add('ω',2*pi);
  kvs.add('dX',deltaX);
  kvs.add('q',1);
  kvs.add('me',1);
  dT:=-1;
  sDT:=-1;
  sDX:=-1;
  tEnde:=100;
  breite:=10.0;
  pDNMax:=-1;
  fortschrittsAnzeige:=false;

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

  // Standardeinstellungen Breich 'lichtVon...'
  lichter:=tMyStringlist.create(prot,'lichter');

  // Standardeinstellungen Breich 'teilchen...'
  setlength(teilchen,0);

  repeat
    if not ifile.readln(s) then begin
      prot.schreibe('Unerwartetes Dateiende in '''+inName+'''!',true);
      prot.destroyall;
      halt(1);
    end;

    if s='Dateiende' then
      break;

    if s='allgemein' then begin
      repeat
        if not ifile.readln(s) then begin
          prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich allgemein!',true);
          prot.destroyall;
          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('dX',deltaX);
          continue;
        end;
        if startetMit('zeitschritt ',s) then begin
          dT:=exprtofloat(false,s,kvs,nil);
          kvs.add('dT',dT);
          continue;
        end;
        if startetMit('diffusionsterm ',s) then begin
          pDNMax:=exprtofloat(false,s,kvs,nil);
          continue;
        end;
        if startetMit('zeit ',s) then begin
          tEnde:=exprtofloat(false,s,kvs,nil);
          continue;
        end;
        if startetMit('breite ',s) then begin
          breite:=exprtofloat(false,s,kvs,nil);
          continue;
        end;
        if s='mit Fortschrittsanzeige' then begin
          fortschrittsAnzeige:=true;
          continue;
        end;
        if s='ohne Fortschrittsanzeige' then begin
          fortschrittsAnzeige:=false;
          continue;
        end;
        prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich allgemein!',true);
        prot.destroyall;
        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);
          prot.destroyall;
          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('ortsschritt ',s) then begin
          sDX:=exprtofloat(false,s,kvs,nil);
          continue;
        end;
        if startetMit('felder ',s) then begin
          s:=s+',';
          while s<>'' do begin
            t:=erstesArgument(s,',');
            setlength(ausgabeDateien,length(ausgabeDateien)+1);
            ausgabeDateien[length(ausgabeDateien)-1]:=tAusgabeDatei.create(t,aPrefix,aSuffix,prot,'ausgabeDateien['+inttostr(length(ausgabeDateien)-1)+']');
          end;
          continue;
        end;
        prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich ausgaben!',true);
        prot.destroyall;
        halt(1);
      until false;
      continue;
    end;

    if startetMit('licht von ',s) then begin
      if startetMit('links ',s) then begin
        lichter.add('links '+s);
        continue;
      end;
      if startetMit('rechts ',s) then begin
        lichter.add('rechts '+s);
        continue;
      end;
      prot.schreibe('Licht kann nicht von '''+erstesArgument(s)+''' kommen!',true);
      prot.destroyall;
      halt(1);
    end;

    if startetMit('teilchen',s) then begin
      if (s<>'') and (strtoint(s)<>length(teilchen)+1) then begin
        prot.schreibe('Ich erwarte die Teilchen beginnend bei 1 der Reihe nach in Parameterdatei '''+inName+'''!',true);
        prot.destroyall;
        halt(1);
      end;
      setlength(teilchen,length(teilchen)+1);
      teilchen[length(teilchen)-1]:=tTeilchenSpezies.create;
      repeat
        if not ifile.readln(s) then begin
          prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich teilchen!',true);
          prot.destroyall;
          halt(1);
        end;
        if s='teilchen'+inttostr(length(teilchen))+'Ende' then break;
        if startetMit('spezifische Ladung ',s) then begin
          teilchen[length(teilchen)-1].spezifischeLadung:=exprToFloat(false,s,kvs,nil);
          kvs.add('qm'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].spezifischeLadung);
          continue;
        end;
        if startetMit('maximaldichte ',s) then begin
          teilchen[length(teilchen)-1].nMax:=exprToFloat(false,s,kvs,nil);
          kvs.add('nmax'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].nMax);
          continue;
        end;
        if startetMit('verteilung ',s) then begin
          teilchen[length(teilchen)-1].liesDichteFunktion(s,ifile,kvs,teilchen,prot);
          continue;
        end;

        prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich teilchen!',true);
        prot.destroyall;
        halt(1);
      until false;
      continue;
    end;

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

  ifile.free;

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

  if dT<0 then
    dT:=deltaX/10;
  if sDT<0 then
    sDT:=dT;
  sDT:=1/round(1/sDT);
  sT:=sDT/2;
  if sDX<0 then
    sDX:=deltaX;

  pro:=tProtokollant.create(prot,'create');
  case zeitverfahren of
    zfEulerVorwaerts:
      pro.schreibe('Iteration mittels Euler-Vorwärts');
    zfRungeKuttaVier:
      pro.schreibe('Iteration mittels Runge-Kutta-4');
    else
      pro.schreibe('Iteration mittels unbekanntem Verfahren');
  end{of case};
  pro.schreibe('dX = '+floattostr(deltaX));
  pro.schreibe('dT = '+floattostr(dT));
  pro.schreibe('breite = '+floattostr(breite)+' (= '+inttostr(round(breite/deltaX)+1)+' Punkte)');
  pro.schreibe('pDNMax = '+floattostr(pDNMax));
  pro.schreibe('bekannte Werte:');
  kvs.dump(pro,'  ');
  if length(teilchen)>0 then begin
    pro.schreibe('teilchen:');
    for i:=0 to length(teilchen)-1 do
      teilchen[i].dump(pro,'  ');
  end;
  if lichter.count>0 then begin
    pro.schreibe('lichter:');
    lichter.dump(pro,'  ');
  end;
  if length(ausgabeDateien)>0 then begin
    pro.schreibe('ausgaben:');
    for i:=0 to length(ausgabeDateien)-1 do begin
      ausgabeDateien[i].zeitAnz:=round(1/sDT);
      ausgabeDateien[i].sDT:=sDT;
      pro.schreibe('  '+ausgabeDateien[i].dump);
    end;
  end;
  pro.free;

  gitter:=tGitter.create(round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,teilchen,lichter,zeitverfahren,prot,'gitter');
  for i:=0 to length(teilchen)-1 do
    teilchen[i].free;
  setlength(teilchen,0);
  lichter.free;
end;

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

function tSimulation.iteriereSchritt(start: extended; var zeitPhysik,zeitDatei: extended): boolean; // noch nicht zu Ende?
var
  i: longint;
begin
  result:=false;

  zeitPhysik:=zeitPhysik-now;
  if errorCode<2 then
    gitter.iteriereSchritt(dT);
  zeitPhysik:=zeitPhysik+now;
  zeitDatei:=zeitDatei-now;
  while (gitter.t>=sT) and (sT<tEnde) do begin
    sT:=sT+sDT;
    for i:=0 to length(ausgabeDateien)-1 do
      ausgabeDateien[i].speichereWerte(gitter,sDX);
  end;
  zeitDatei:=zeitDatei+now;

  if fortschrittsAnzeige then begin
    if (floor(100*Gitter.t/tEnde) < floor(100*(Gitter.t+dT)/tEnde)) then begin
      prot.schreibe(inttostr(round(100*Gitter.t/tEnde))+'% (t='+floattostr(Gitter.t)+'T)',true);
      prot.schreibe(timetostr(now-start)+' ('+floattostr(zeitPhysik/max(1e-11,zeitPhysik+zeitDatei))+')',true);
      prot.schreibe('ETA: '+timetostr((now-start)*(tEnde-Gitter.t)/max(Gitter.t,dT)),true);
      prot.schreibe('aktueller Zeitschritt: '+floattostr(dT)+'T',true);
      if errorCode<2 then
        if not gitter.dumpErhaltungsgroessen then begin
          errorCode:=2;
          exit;
        end;
    end;
  end;

  result:=gitter.t<tEnde;
end;

end.