path: root/Physikunit.pas

unit Physikunit;

{$mode objfpc}{$H+}

{$IFNDEF UNIX}
  {$ERROR This program can be compiled only on/for Unix/Linux based systems.}
{$ENDIF}

{$DEFINE Zeitschrittueberwachung}
{$DEFINE Dichteueberwachung}
{ $DEFINE negativeDichteueberwachung}

interface

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

type
  tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaDreiAchtel,zfRungeKuttaVier,zfRungeKuttaZehn,zfRungeKuttaZwoelf,zfRungeKuttaVierzehn);
  tVerteilungsfunktion = function(x: extended): extended;
  tEMFeldGroesze = (
    efAX,efAY,efAZ,
    efDAXDT,efDAYDT,efDAZDT,
    efEX,efBY,efBZ
  );
  tEMQuellGroesze = (
    eqRho,eqJX,eqJY,eqJZ
  );
  tEMQuellen = array[tEMQuellGroesze] of extended;

const
  erstesEMFMitAbleitung  =  efAX;
  letztesEMFMitAbleitung =  efDAZDT;

type

  tRaumPunkt = class;
  tFelder = class;
  tGitter = class;
  tSimulation = class;

  { tAusgabeDatei }

  tAusgabeDatei = class
  private
    ableitung:          boolean;
    emFelder:           tEMFeldGroesze;
    emQuellen:          tEMQuellGroesze;
    teilchen,           // -2: EM-Feld; -1: Gesamt-EM-Quelle; sonst: EM-Quelle der entsprechenden Teilchen
      nNum,tCnt:        longint;
    pre,suf:            string;
    datei:              file;
    pro:                tProtokollant;
    buf:                pointer;
    bufLen,bufPos:      longint;
    procedure schreibeBuffer(minBufLen: longint);
  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; sT, sDX: extended);
  end;

  { tTeilchenSpezies }

  tTeilchenSpezies = class
  private
    ortsGrenzen:   array of extended;
    dichteStuecke: array of string;
  public
    nMin,nMax,                   // minimale und 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;

  { tImpulsPunkt }

  tImpulsPunkt = class(tObject) // Besetzungsdichten verschiedener Teilchenspezies für ein Paar (r;p)
  private
    raumpunkt: tRaumPunkt;
    nachbarn:  array[0..3,boolean] of tImpulsPunkt; // x-, px-, py-, pz-Nachbarn (jeweil in Richtung - und +)
  public
    werte: array of array[boolean] of extended; // Besetzungsdichten und deren Ableitungen für verschiedene Teilchenspezies
    constructor create(rp: tRaumPunkt; linksRaum,linksImpuls,untenImpuls,vorneImpuls: tImpulsPunkt; anzTeilchen: longint);
    destructor destroy; override;
    {$DEFINE LiKotImpulsPunktHeader}
    {$INCLUDE linearkombinationen.inc}
    {$UNDEF LiKotImpulsPunktHeader}
    procedure nichtnegativieren;
    procedure akkumuliereEMQuellen(var emQuellen: tEMQuellen; vx,vy,vz,dVP: extended);
    function gibMirWas(var defizit: extended; entfernung, teilchen,richtung: longint; positiv: boolean): boolean; // ich habe jemanden erreicht
    procedure setzeNull;
    procedure berechneAbleitungen(pX,pY,pZ,iDX,iDPX,iDPY,iDPZ: extended);
  end;

  { tRaumPunkt }

  tRaumPunkt = class(tObject) // repräsentiert alle Werte eines Punktes im Raumgitter und deren Zeitableitungen
  private
    felder: tFelder;
  public
    phasenraum: array of tImpulsPunkt;
      // Phasenraumbesetzungsdichten
    emFelder:  array[tEMFeldGroesze,boolean] of extended;
      // EM-Felder und deren Zeitableitungen
      // dPhiDX und B[xyz] haben keine sinnvolle Ableitung hier!
    emQuellen: array[tEMQuellGroesze] of extended;
      // Quellen für die EM-Felder
    lN,rN:    tRaumPunkt;
      // Nachbarpunkte
    dPX,dPY,dPZ: extended;
      // Impulsdiskretisierung
    aPX,aPY,aPZ: longint;
      // Anzahl der Impulsbins

    constructor create(linkerNachbar: tRaumPunkt; derBesitzer: tFelder; teilchen,_aPX,_aPY,_aPZ: longint; _dPX,_dPY,_dPZ: extended);
    destructor destroy; override;
    procedure berechneEMFelder(dX: extended; rechts: boolean); overload;
    procedure berechneEMFelder(dX,iDX: extended); overload;
    procedure berechneEMAbleitungen(iDX: extended);
    procedure berechnePhasenraumAbleitungen(iDX: extended);
    {$DEFINE LiKotRaumPunktHeader}
    {$INCLUDE linearkombinationen.inc}
    {$UNDEF LiKotRaumPunktHeader}
    procedure nichtnegativieren;
    function impulsIntegral(teilchen: longint; emQ: tEMQuellGroesze): extended;
    procedure initialisiereDichte(teilchen: longint; n: extended);
    procedure reflektiereTeilchen(rechts: boolean);
    procedure setzeNull;
  end;

  { TFelder }

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

    constructor create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter; aPX,aPY,aPZ: longint; dPX,dPY,dPZ: extended);
    destructor destroy; override;
    procedure berechneAbleitungen(dX,iDX: extended);
    {$DEFINE LiKotFelderHeader}
    {$INCLUDE linearkombinationen.inc}
    {$UNDEF LiKotFelderHeader}
    procedure setzeNull;
  end;

  { tGitter }

  tGitter = class(tObject)
  private
    besitzer:                     tSimulation;
    prot:                         tProtokollant;
    zeitverfahren:                tZeitverfahren;
    abbruch:                      boolean;
    dX,iDX,xl,t:                  extended;
    kvs:                          tKnownValues;
    procedure abbrechen;
  public

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

    constructor create(derBesitzer: tSimulation; size: longint; deltaX: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string; aPX,aPY,aPZ: longint; dPX,dPY,dPZ: extended);
    destructor destroy; override;
    procedure iteriereSchritt(dT: extended);
  end;

  { tSimulation }

  tSimulation = class(tObject)
  private
    prot:                      tProtokollant;
    gitter:                    tGitter;
    dT,tEnde,sT,sDT,sDX:       extended;
    fortschrittsAnzeige:       boolean;
    ausgabeDateien:            array of tAusgabeDatei;
  public
    gotSighup:                 boolean;
    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;

procedure SignalCapture(signal : longint); cdecl;

implementation

const
  emFeldNamen: array[tEMFeldGroesze] of string = (
    'AX','AY','AZ','DAXDT','DAYDT','DAZDT','EX','BY','BZ'
  );
  emQuellNamen: array[tEMQuellGroesze] of string = (
    'RHO','JX','JY','JZ'
  );
  minAusgabeBuffer = 1024*1024;

var
  sighupSimulationen: array of tSimulation;

// tAusgabeDatei ***************************************************************

constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string);
var
  emF:     tEMFeldGroesze;
  emQ:     tEMQuellGroesze;
  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:=-1;

  emFelder:=efAX;
  emQuellen:=eqRho;
  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(tEMFeldGroesze) to high(tEMFeldGroesze) do
        if copy('D',1,byte(abl))+emFeldNamen[emF] = feldName then begin
          emFelder:=emF;
          ableitung:=abl;
          teilchen:=-2;
          gef:=true;
          break;
        end;
  end;
  if not gef then
    for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do
      if emQuellNamen[emQ] = feldName then begin
        ableitung:=false;
        emQuellen:=emQ;
        gef:=true;
        break;
      end;

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

  if teilchen=-2 then
    pre:=emFeldNamen[emF]
  else begin
    pre:=emQuellNamen[emQ];
    if teilchen>=0 then
      pre:=pre+inttostr(teilchen+1);
  end;

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

  bufLen:=0;
  buf:=nil;
  bufPos:=0;

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

destructor tAusgabeDatei.destroy;
begin
  if bufPos<>0 then
    schreibeBuffer(0);
  freeMem(buf,bufLen);
  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;

procedure tAusgabeDatei.schreibeBuffer(minBufLen: longint);
begin
  if buf=nil then begin
    bufLen:=max(minAusgabeBuffer,minBufLen);
    getmem(buf,bufLen);
  end
  else begin
    reset(datei,1);
    seek(datei,fileSize(datei));
    blockWrite(datei,buf^,bufPos);
  end;
  bufPos:=0;
end;

function tAusgabeDatei.dump: string;
begin
  if teilchen=-2 then
    result:=emFeldNamen[emFelder]
  else begin
    result:=emQuellNamen[emQuellen];
    if teilchen>=0 then
      result:=result+inttostr(teilchen);
  end;
  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;
  if bufPos<>0 then
    schreibeBuffer(0);
  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; sT, sDX: extended);
var i,cnt:     longint;
    sX,cX,val: 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 sT>=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);
  if bufLen-bufPos < sizeof(integer) + (2+cnt)*sizeof(extended) then
    schreibeBuffer(sizeof(integer) + (2+cnt)*sizeof(extended));

  if bufLen-bufPos < sizeof(integer) + (2+cnt)*sizeof(extended) then begin
    pro.schreibe('Schreibbuffer ist zu klein ('+inttostr(bufLen)+' statt mindestens '+inttostr(sizeof(integer) + (2+cnt)*sizeof(extended))+' Bytes)!');
    pro.destroyall;
    halt(1);
  end;

  move(cX,(buf+bufPos)^,sizeof(extended));
  inc(bufPos,sizeof(extended));
  move(sX,(buf+bufPos)^,sizeof(extended));
  inc(bufPos,sizeof(extended));
  move(cnt,(buf+bufPos)^,sizeof(integer));
  inc(bufPos,sizeof(integer));

  sX:=cX-Min(gitter.dX,sDX)/2;
  if teilchen=-2 then begin
    for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
      while cX>=sX do begin
        dec(cnt);
        move(gitter.felders[gitter.aktuelleFelder].inhalt[i].emFelder[emFelder,ableitung],(buf+bufPos)^,sizeof(extended));
        inc(bufPos,sizeof(extended));
        sX:=sX+sDX;
      end;
      cX:=cX+gitter.dX;
    end;
  end
  else if teilchen=-1 then begin
    for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
      while cX>=sX do begin
        dec(cnt);
        move(gitter.felders[gitter.aktuelleFelder].inhalt[i].emQuellen[emQuellen],(buf+bufPos)^,sizeof(extended));
        inc(bufPos,sizeof(extended));
        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);
        val:=gitter.felders[gitter.aktuelleFelder].inhalt[i].impulsIntegral(teilchen,emQuellen);
        move(val,(buf+bufPos)^,sizeof(extended));
        inc(bufPos,sizeof(extended));
        sX:=sX+sDX;
      end;
      cX:=cX+gitter.dX;
    end;
  end;
  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);
  nMin:=1E-3;
  nMax:=nMin;
  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-nMin)*exprToFloat(false,dichteStuecke[i],kvs,nil)+nMin;
  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;

// tImpulsPunkt ****************************************************************

constructor tImpulsPunkt.create(rp: tRaumPunkt; linksRaum,linksImpuls,untenImpuls,vorneImpuls: tImpulsPunkt; anzTeilchen: longint);
var
  i: longint;
  b: boolean;
begin
  inherited create;
  raumPunkt:=rp;
  fillchar(werte,sizeof(werte),#0);
  setlength(werte,anzTeilchen);
  for i:=0 to length(werte)-1 do
    for b:=false to true do
      werte[i,b]:=0;
  nachbarn[0,false]:=linksRaum;
  nachbarn[0,true]:=nil;
  nachbarn[1,false]:=linksImpuls;
  nachbarn[1,true]:=nil;
  nachbarn[2,false]:=untenImpuls;
  nachbarn[2,true]:=nil;
  nachbarn[3,false]:=vorneImpuls;
  nachbarn[3,true]:=nil;
  for i:=0 to 3 do
    if assigned(nachbarn[i,false]) then
      nachbarn[i,false].nachbarn[i,true]:=self;
end;

destructor tImpulsPunkt.destroy;
var
  i: longint;
  b: boolean;
begin
  for i:=0 to 3 do
    for b:=false to true do
      if assigned(nachbarn[i,b]) then
        nachbarn[i,b].nachbarn[i,not b]:=nil;
  setlength(werte,0);
  inherited destroy;
end;

procedure tImpulsPunkt.nichtnegativieren; // Dichten nicht negativ machen
var
  i,j,entfernung:   longint;
  defizit:          extended;
  pro:              tProtokollant;
  jemandErreicht,b: boolean;
begin
  for i:=0 to length(werte)-1 do
    if werte[i,false]<0 then begin
      defizit:=-werte[i,false];
      werte[i,false]:=0;
      entfernung:=1;
      jemandErreicht:=true;
      while (defizit>0) and jemandErreicht do begin
        jemandErreicht:=false;
        for j:=0 to 3 do
          for b:=false to true do
            if assigned(nachbarn[j,b]) then
              jemandErreicht:=nachbarn[j,b].gibMirWas(defizit,i,entfernung,j,b) or jemandErreicht;

        inc(entfernung);
      end;
      if defizit>0 then begin
        pro:=tProtokollant.create(raumpunkt.felder.gitter.prot,'WertePunkt.nichtnegativieren');
        pro.schreibe('Ich konnte ein Teilchendefizit der Sorte '+inttostr(i+1)+' nicht auffüllen!',true);
        raumpunkt.felder.gitter.abbrechen;
      end;
    end;
end;

procedure tImpulsPunkt.akkumuliereEMQuellen(var emQuellen: tEMQuellen; vx,vy,vz,dVP: extended);
var
  i:   longint;
  tmp: extended;
begin
  for i:=0 to length(werte)-1 do begin
    tmp:=werte[i,false] * raumpunkt.felder.spezLadungen[i] * dVP;
    emQuellen[eqRho]:=emQuellen[eqRho] + tmp;
    emQuellen[eqJX]:= emQuellen[eqJX]  + tmp * vx;
    emQuellen[eqJY]:= emQuellen[eqJY]  + tmp * vy;
    emQuellen[eqJZ]:= emQuellen[eqJZ]  + tmp * vz;
  end;
end;

function tImpulsPunkt.gibMirWas(var defizit: extended; entfernung, teilchen,richtung: longint; positiv: boolean): boolean; // ich habe jemanden erreicht
var
  i: longint;
begin
  result:=entfernung<=0;
  if result then begin
    if (defizit>0) and (werte[teilchen,false]>0) then begin
      if defizit>werte[teilchen,false] then begin
        defizit:=defizit-werte[teilchen,false];
        werte[teilchen,false]:=0;
      end
      else begin
        werte[teilchen,false]:=werte[teilchen,false]-defizit;
        defizit:=0;;
      end;
    end;
  end
  else
    if defizit>0 then begin
      result:=false;
      for i:=richtung to 4 do
        if assigned(nachbarn[i,positiv]) then
          result:=nachbarn[i,positiv].gibMirWas(defizit,entfernung-1,teilchen,i,positiv) or result;
    end;
end;

procedure tImpulsPunkt.setzeNull;
var
  i:   longint;
  abl: boolean;
begin
  for i:=0 to length(werte)-1 do
    for abl:=false to true do
      werte[i,abl]:=0;
end;

procedure tImpulsPunkt.berechneAbleitungen(pX,pY,pZ,iDX,iDPX,iDPY,iDPZ: extended);
var
  igam: extended;
  i:    longint;
begin
  igam:=1/sqrt(1+sqr(pX)+sqr(pY)+sqr(pZ));
  pX:=pX*igam;
  pY:=pY*igam;
  pZ:=pZ*igam;

  // df/dt = - v Nabla f - q(E + v/c x B) Nabla_p f
  for i:=0 to length(werte)-1 do
    werte[i,true]:=
      - pX * (nachbarn[0,true].werte[i,false] - nachbarn[0,false].werte[i,false]) * iDX/2
      - raumpunkt.felder.spezLadungen[i] * (
        (raumpunkt.emFelder[efEX,false] + pY*raumpunkt.emFelder[efBZ,false] - pZ*raumpunkt.emFelder[efBY,false])
         * (nachbarn[1,true].werte[i,false] - nachbarn[1,false].werte[i,false]) * iDPX/2
        -pX*raumpunkt.emFelder[efBZ,false]
         * (nachbarn[2,true].werte[i,false] - nachbarn[2,false].werte[i,false]) * iDPY/2
        +pX*raumpunkt.emFelder[efBY,false]
         * (nachbarn[3,true].werte[i,false] - nachbarn[3,false].werte[i,false]) * iDPZ/2
      );
end;

// tRaumPunkt ******************************************************************

constructor tRaumPunkt.create(linkerNachbar: tRaumPunkt; derBesitzer: tFelder; teilchen,_aPX,_aPY,_aPZ: longint; _dPX,_dPY,_dPZ: extended);
var
  emF:            tEMFeldGroesze;
  emQ:            tEMQuellGroesze;
  abl:            boolean;
  i,j,k:          longint;
  xN,pxN,pyN,pzN: tImpulsPunkt;
begin
  inherited create;
  felder:=derBesitzer;
  for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do
    emQuellen[emQ]:=0;
  for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do
    for abl:=false to true do
      emFelder[emF,abl]:=0;
  lN:=linkerNachbar;
  rN:=nil;
  if assigned(lN) then
    lN.rN:=self;
  aPX:=_aPX+byte(not odd(_aPX));
  aPY:=_aPY+byte(not odd(_aPY));
  aPZ:=_aPZ+byte(not odd(_aPZ));
  dPX:=_dPX;
  dPY:=_dPY;
  dPZ:=_dPZ;
  fillchar(phasenraum,sizeof(phasenraum),#0);
  setlength(phasenraum,aPX*aPY*aPZ);
  for i:=0 to aPX-1 do
    for j:=0 to aPY-1 do
      for k:=0 to aPZ-1 do begin
        if assigned(linkerNachbar) then
          xN:=linkerNachbar.phasenraum[i+aPX*(j+k*aPY)]
        else
          xN:=nil;
        if i>0 then
          pxN:=phasenraum[i-1+aPX*(j+k*aPY)]
        else
          pxN:=nil;
        if j>0 then
          pyN:=phasenraum[i+aPX*(j-1+k*aPY)]
        else
          pyN:=nil;
        if k>0 then
          pzN:=phasenraum[i+aPX*(j+(k-1)*aPY)]
        else
          pzN:=nil;
        phasenraum[i+aPX*(j+k*aPY)]:=tImpulsPunkt.create(self,xN,pxN,pyN,pzN,teilchen);
      end;
end;

destructor tRaumPunkt.destroy;
var
  i: longint;
begin
  for i:=0 to length(phasenraum)-1 do
    phasenraum[i].free;
  setlength(phasenraum,0);
  if assigned(lN) then
    lN.rN:=nil;
  if assigned(rN) then
    rN.lN:=nil;
  inherited destroy;
end;

procedure tRaumPunkt.berechneEMFelder(dX: extended; rechts: boolean);
begin
  // dEX/dx = rho  ... hier muss integriert werden:
  // EX = EX(x- deltax) + <rho> * deltax

  if rechts then
    emFelder[efEX,false]:=
      lN.emFelder[efEX,false] +
      (lN.emQuellen[eqRho] + emQuellen[eqRho]) * dX/2
  else
    emFelder[efEX,false]:=
      emQuellen[eqRho] * dX;

  // B = rot A
  emFelder[efBY,false]:=0;  // kann man vllt auch hübscher machen, ich bin aber gerade zu faul dafür
  emFelder[efBZ,false]:=0;
end;

procedure tRaumPunkt.berechneEMFelder(dX,iDX: extended);
begin
  // dEX/dx = rho  ... hier muss integriert werden:
  // EX = EX(x- deltax) + <rho> * deltax

  emFelder[efEX,false]:=
    lN.emFelder[efEX,false] +
    (lN.emQuellen[eqRho] + emQuellen[eqRho]) * dX/2;

  // B = rot A
  emFelder[efBY,false]:=
    (lN.emFelder[efAZ,false]-rN.emFelder[efAZ,false])*iDX/2;
  emFelder[efBZ,false]:=
    (rN.emFelder[efAY,false]-lN.emFelder[efAY,false])*iDX/2;
end;

procedure tRaumPunkt.berechneEMAbleitungen(iDX: extended); // Zeitableitungen der EM-Potentiale berechnen
var
  i,j,k: longint;
  emQ:   tEMQuellGroesze;
begin
  for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do
    emQuellen[emQ]:=0;
  for i:=0 to aPX-1 do
    for j:=0 to aPY-1 do
      for k:=0 to aPZ-1 do
        phasenraum[i+aPX*(j+aPY*k)].akkumuliereEMQuellen(emQuellen,(i-aPX/2)*dPX,(j-aPY/2)*dPY,(k-aPZ/2)*dPZ,dPX*dPY*dPZ);

  // d2A/dt2 = Laplace(A) - j    ( - dE/dt wird auf dA/dt direkt aufgeschlagen !!! )
  emFelder[efDAXDT,true]:=
    ( rN.emFelder[efAX,false] - 2*emFelder[efAX,false] + lN.emFelder[efAX,false] )*sqr(iDX)
    - emQuellen[eqJX];
  emFelder[efDAYDT,true]:=
    ( rN.emFelder[efAY,false] - 2*emFelder[efAY,false] + lN.emFelder[efAY,false] )*sqr(iDX)
    - emQuellen[eqJY];
  emFelder[efDAZDT,true]:=
    ( rN.emFelder[efAZ,false] - 2*emFelder[efAZ,false] + lN.emFelder[efAZ,false] )*sqr(iDX)
    - emQuellen[eqJZ];

  // dA/dt = dA/dt - E
  emFelder[efAX,true]:=emFelder[efDAXDT,false] - emFelder[efEX,false];
  emFelder[efAY,true]:=emFelder[efDAYDT,false];
  emFelder[efAZ,true]:=emFelder[efDAZDT,false];
end;

procedure tRaumPunkt.berechnePhasenraumAbleitungen(iDX: extended);
var
  i,j,k: longint;
begin
  // df/dt = - v Nabla f - q(E + v/c x B) Nabla_p f
  for i:=1 to aPX-2 do
    for j:=1 to aPY-2 do
      for k:=1 to aPZ-2 do
        phasenraum[i+aPX*(j+aPY*k)].berechneAbleitungen((i-(aPX div 2))*dPX,(j-(aPY div 2))*dPY,(k-(aPZ div 2))*dPZ,iDX,1/dPX,1/dPY,1/dPZ);
end;

procedure tRaumPunkt.nichtnegativieren; // Dichten nicht negativ machen
var
  i:       longint;
begin
  for i:=0 to length(phasenraum)-1 do
    phasenraum[i].nichtnegativieren;
end;

function tRaumPunkt.impulsIntegral(teilchen: longint; emQ: tEMQuellGroesze): extended;
var
  i,j,k:    longint;
begin
  if teilchen<0 then begin
    result:=emQuellen[emQ]; // das ist leicht :-)
    exit;
  end;

  result:=0;

  case emQ of
    eqRho:
      for i:=0 to aPX-1 do
        for j:=0 to aPY-1 do
          for k:=0 to aPZ-1 do
            result:=
              result +
              phasenraum[i+aPX*(j+aPY*k)].werte[teilchen,false]*dPX*dPY*dPZ;
    eqJX:
      for i:=0 to aPX-1 do
        for j:=0 to aPY-1 do
          for k:=0 to aPZ-1 do
            result:=
              result +
              phasenraum[i+aPX*(j+aPY*k)].werte[teilchen,false]*dPX*dPY*dPZ * (i-aPX/2)*dPX;
    eqJY:
      for i:=0 to aPX-1 do
        for j:=0 to aPY-1 do
          for k:=0 to aPZ-1 do
            result:=
              result +
              phasenraum[i+aPX*(j+aPY*k)].werte[teilchen,false]*dPX*dPY*dPZ * (j-aPY/2)*dPY;
    eqJZ:
      for i:=0 to aPX-1 do
        for j:=0 to aPY-1 do
          for k:=0 to aPZ-1 do
            result:=
              result +
              phasenraum[i+aPX*(j+aPY*k)].werte[teilchen,false]*dPX*dPY*dPZ * (k-aPZ/2)*dPZ;
  end;
end;

procedure tRaumPunkt.initialisiereDichte(teilchen: longint; n: extended);
var
  i: longint;
begin
  for i:=0 to length(phasenraum)-1 do
    if (aPX div 2) + aPX*((aPY div 2) + aPY*(aPZ div 2)) = i then
      phasenraum[i].werte[teilchen,false]:=n
    else
      phasenraum[i].werte[teilchen,false]:=0;
end;

procedure tRaumPunkt.reflektiereTeilchen(rechts: boolean);
var
  i,j,k: longint;
begin
  for i:=(aPX div 2 + 1)*byte(rechts) to (aPX div 2)-1 + (apx-(aPX div 2))*byte(rechts) do
    for j:=0 to aPY*aPZ-1 do
      for k:=0 to felder.matAnz-1 do
        if phasenraum[i+aPX*j].werte[k,false]>0 then begin
          phasenraum[aPX-1-i+aPX*j].werte[k,false]:=phasenraum[i+aPX*j].werte[k,false];
          phasenraum[i+aPX*j].werte[k,false]:=0;
        end;
end;

procedure tRaumPunkt.setzeNull;
var
  i:   longint;
  emF: tEMFeldGroesze;
  emQ: tEMQuellGroesze;
  abl: boolean;
begin
  for i:=0 to length(phasenraum)-1 do
    phasenraum[i].setzeNull;
  for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do
    emQuellen[emQ]:=0;
  for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do
    for abl:=false to true do
      emFelder[emF,abl]:=0;
end;

// linearkombinationsspezifische Methoden von tImpulsPunkt, tRaumPunkt und tFelder

{$DEFINE LiKoimplementation}
{$INCLUDE linearkombinationen.inc}
{$UNDEF LiKoimplementation}

// tFelder *********************************************************************

constructor tFelder.create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter; aPX,aPY,aPZ: longint; dPX,dPY,dPZ: extended);
var
  i,j:    longint;
  rechts: boolean;
  dens:   extended;
begin
  inherited create;
  gitter:=parent;
  matAnz:=length(teilchen);
  fillchar(spezLadungen,sizeof(spezLadungen),#0);
  setlength(spezLadungen,matAnz);
  for i:=0 to length(spezLadungen)-1 do
    spezLadungen[i]:=teilchen[i].spezifischeLadung;
  fillchar(inhalt,sizeof(inhalt),#0);
  setlength(inhalt,groesse+4); // zwei Felder links und rechts extra für Randbedingungen
  inhalt[0]:=tRaumPunkt.create(nil,self,length(teilchen),aPX,aPY,aPZ,dPX,dPY,dPZ);
  for i:=1 to length(inhalt)-1 do
    inhalt[i]:=tRaumPunkt.create(inhalt[i-1],self,length(teilchen),aPX,aPY,aPZ,dPX,dPY,dPZ);
  fillchar(massen,sizeof(massen),#0);
  setlength(massen,length(teilchen));
  for i:=0 to length(massen)-1 do
    massen[i]:=0;
  for i:=0 to length(inhalt)-1 do
    for j:=0 to length(teilchen)-1 do begin
      dens:=teilchen[j].gibDichte(gitter.xl+(i-2)*gitter.dX,parent.kvs);
      inhalt[i].initialisiereDichte(j,dens);
      massen[j]:=massen[j]+dens*gitter.dX;
    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;
var
  i: longint;
begin
  setlength(spezLadungen,0);
  for i:=0 to length(inhalt)-1 do
    inhalt[i].free;
  setlength(inhalt,0);
  lichters[false].free;
  lichters[true].free;
  inherited destroy;
end;

procedure tFelder.setzeRaender(iDX: extended);
var
  emF:    tEMFeldGroesze;
  rechts: boolean;
  i:      longint;
begin
  for emF:=efAX to efAZ do begin // Vakuumrandbedingungen für das A-Feld
    inhalt[0].emFelder[emF,true]:=
      (inhalt[1].emFelder[emF,false] -
       inhalt[0].emFelder[emF,false])*iDX;
    inhalt[length(inhalt)-1].emFelder[emF,true]:=
      (inhalt[length(inhalt)-2].emFelder[emF,false] -
       inhalt[length(inhalt)-1].emFelder[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
      inhalt[(length(inhalt)-1)*byte(rechts)].emFelder[efAY,true]:=
        inhalt[(length(inhalt)-1)*byte(rechts)].emFelder[efAY,true] +
        exprToFloat(false,lichters[rechts][i],gitter.kvs,nil);
end;

procedure tFelder.berechneAbleitungen(dX,iDX: extended);
var
  i: longint;
  b: boolean;
begin
  for b:=false to true do
    inhalt[(length(inhalt)-1)*byte(b)].reflektiereTeilchen(b);

  inhalt[0].berechneEMFelder(dX,false);
  for i:=1 to length(inhalt)-2 do
    inhalt[i].berechneEMFelder(dX,iDX);
  inhalt[length(inhalt)-1].berechneEMFelder(dX,false);

  setzeRaender(iDX);

  for i:=1 to length(inhalt)-2 do
    inhalt[i].berechneEMAbleitungen(iDX);
  for i:=1 to length(inhalt)-2 do
    inhalt[i].berechnePhasenraumAbleitungen(iDX);
end;

procedure tFelder.setzeNull;
var
  i: longint;
begin
  for i:=0 to length(inhalt)-1 do
    inhalt[i].setzeNull;
end;

// tGitter *********************************************************************

constructor tGitter.create(derBesitzer: tSimulation; size: longint; deltaX: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string; aPX,aPY,aPZ: longint; dPX,dPY,dPZ: extended);
var
  i: longint;
begin
  inherited create;
  abbruch:=false;
  besitzer:=derBesitzer;
  zeitverfahren:=zv;
  kvs:=bekannteWerte;
  prot:=tProtokollant.create(protokollant,name);
  dX:=deltaX;
  iDX:=1/dX;
  case Zeitverfahren of
    zfEulerVorwaerts:
      Setlength(Felders,2);
    zfRungeKuttaDreiAchtel,zfRungeKuttaVier:
      Setlength(Felders,5);
    zfRungeKuttaZehn:
      Setlength(Felders,18);
    zfRungeKuttaZwoelf:
      Setlength(Felders,26);
    zfRungeKuttaVierzehn:
      Setlength(Felders,36);
  end{of Case};
  xl:=dX/2;

  for i:=0 to length(felders)-1 do
    felders[i]:=tFelder.create(size,teilchen,lichter,self,aPX,aPY,aPZ,dPX,dPY,dPZ);

  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.abbrechen;
var
  i: longint;
begin
  abbruch:=true;
  for i:=0 to length(felders)-1 do
    felders[i].setzeNull;
end;

procedure tGitter.iteriereSchritt(dT: extended);
    {$IFDEF negativeDichteueberwachung}
var i,j:     longint;
    pro:     tProtokollant;
    {$ENDIF}
begin
  if abbruch then begin
    t:=t+dT;
    exit;
  end;
  kvs.add('t',t);

  repeat
    felders[aktuelleFelder].berechneAbleitungen(dX,iDX); // y' = y'(t,y(t))

    case zeitverfahren of
      zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt
        felders[1-aktuelleFelder].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT);
      zfRungeKuttaDreiAchtel: begin
        {$INCLUDE rk3_8.inc}
      end;
      zfRungeKuttaVier: begin
        {$INCLUDE rk4.inc}
      end;
      zfRungeKuttaZehn: begin      // Quelle:  http://sce.uhcl.edu/rungekutta/rk108.txt
        {$INCLUDE rk108.inc}
      end;
      zfRungeKuttaZwoelf: begin    // Quelle:  http://sce.uhcl.edu/rungekutta/rk1210.txt
        {$INCLUDE rk1210.inc}
      end;
      zfRungeKuttaVierzehn: begin  // Quelle:  http://sce.uhcl.edu/rungekutta/rk1412.txt
        {$INCLUDE rk1412.inc}
      end;
    end{of case};

    break;
  until abbruch;

  aktuelleFelder:=1-aktuelleFelder;
  {$IFDEF negativeDichteueberwachung}
  for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do
    for j:=0 to felders[aktuelleFelder].matAnz-1 do
        if felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false]<0 then begin
          pro:=tProtokollant.create(prot,'iteriereSchritt');
          pro.schreibe('n<0 bei:',true);
          pro.schreibe(' t    = ' +floattostr(t)+',',true);
          pro.schreibe(' i    = ' +inttostr(i)+' (x = '+floattostr(xl+i*dX)+'),',true);
          pro.schreibe(' alte Werte:',true);
          with felders[1-aktuelleFelder] do begin
            pro.schreibe(' v_l  = ' +floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]*
                                                inhalt[i-1].matWerte[j,mfIGamma,false])+',',true);
            pro.schreibe(' n_l  = ' +floattostr(inhalt[i-1].matWerte[j,mfN,false])+',',true);
            pro.schreibe(' n_l'' = '+floattostr(inhalt[i-1].matWerte[j,mfN,true])+',',true);
            pro.schreibe(' v    = ' +floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]*
                                                inhalt[i].matWerte[j,mfIGamma,false])+',',true);
            pro.schreibe(' n    = ' +floattostr(inhalt[i].matWerte[j,mfN,false])+',',true);
            pro.schreibe(' n''   = '+floattostr(inhalt[i].matWerte[j,mfN,true])+',',true);
            pro.schreibe(' v_r  = ' +floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]*
                                                inhalt[i+1].matWerte[j,mfIGamma,false])+',',true);
            pro.schreibe(' n_r  = ' +floattostr(inhalt[i+1].matWerte[j,mfN,false])+',',true);
            pro.schreibe(' n_r'' = '+floattostr(inhalt[i+1].matWerte[j,mfN,true]),true);
          end;
          pro.schreibe(' neue Werte:',true);
          with felders[aktuelleFelder] do begin
            pro.schreibe(' v_l  = ' +floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]*
                                                inhalt[i-1].matWerte[j,mfIGamma,false])+',',true);
            pro.schreibe(' n_l  = ' +floattostr(inhalt[i-1].matWerte[j,mfN,false])+',',true);
            pro.schreibe(' n_l'' = '+floattostr(inhalt[i-1].matWerte[j,mfN,true])+',',true);
            pro.schreibe(' v    = ' +floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]*
                                                inhalt[i].matWerte[j,mfIGamma,false])+',',true);
            pro.schreibe(' n    = ' +floattostr(inhalt[i].matWerte[j,mfN,false])+',',true);
            pro.schreibe(' n''   = '+floattostr(inhalt[i].matWerte[j,mfN,true])+',',true);
            pro.schreibe(' v_r  = ' +floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]*
                                                inhalt[i+1].matWerte[j,mfIGamma,false])+',',true);
            pro.schreibe(' n_r  = ' +floattostr(inhalt[i+1].matWerte[j,mfN,false])+',',true);
            pro.schreibe(' n_r'' = '+floattostr(inhalt[i+1].matWerte[j,mfN,true]),true);
          end;
          pro.free;
          abbrechen;
        end;
  {$ENDIF}
  t:=t+dT;
end;

// tSimulation *****************************************************************

constructor tSimulation.create(inName: string; protokollant: tProtokollant; name: string);
var
  ifile:                                       tMyStringlist;
  zeitverfahren:                               tZeitverfahren;
  s,t,aSuffix,aPrefix:                         string;
  deltaX,breite,dPX,dPY,dPZ,pXMax,pYMax,pZMax: extended;
  i:                                           longint;
  kvs:                                         tKnownValues;
  teilchen:                                    array of tTeilchenSpezies;
  lichter:                                     tMyStringlist;
  pro:                                         tProtokollant;
  na:                                          pSigActionRec;
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;
  fortschrittsAnzeige:=false;
  gotSighup:=false;
  dPX:=-1;
  dPY:=-1;
  dPZ:=-1;
  pXMax:=3;
  pYMax:=3;
  pZMax:=3;

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

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

  // Standardeinstellungen Bereich '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-3/8' then begin
          Zeitverfahren:=zfRungeKuttaDreiAchtel;
          continue;
        end;
        if s='runge-Kutta-4' then begin
          Zeitverfahren:=zfRungeKuttaVier;
          continue;
        end;
        if s='runge-Kutta-10' then begin
          Zeitverfahren:=zfRungeKuttaZehn;
          continue;
        end;
        if s='runge-Kutta-12' then begin
          Zeitverfahren:=zfRungeKuttaZwoelf;
          continue;
        end;
        if s='runge-Kutta-14' then begin
          Zeitverfahren:=zfRungeKuttaVierzehn;
          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('impulsschritt ',s) then begin
          if startetMit('x ',s) then begin
            dPX:=exprtofloat(false,s,kvs,nil);
            continue;
          end;
          if startetMit('y ',s) then begin
            dPY:=exprtofloat(false,s,kvs,nil);
            continue;
          end;
          if startetMit('z ',s) then begin
            dPZ:=exprtofloat(false,s,kvs,nil);
            continue;
          end;
          dPX:=exprtofloat(false,s,kvs,nil);
          dPY:=dPX;
          dPZ:=dPX;
          continue;
        end;
        if startetMit('maximalimpuls ',s) then begin
          if startetMit('x ',s) then begin
            pXMax:=exprtofloat(false,s,kvs,nil);
            continue;
          end;
          if startetMit('y ',s) then begin
            pYMax:=exprtofloat(false,s,kvs,nil);
            continue;
          end;
          if startetMit('z ',s) then begin
            pZMax:=exprtofloat(false,s,kvs,nil);
            continue;
          end;
          pXMax:=exprtofloat(false,s,kvs,nil);
          pYMax:=pXMax;
          pZMax:=pXMax;
          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('minimaldichte ',s) then begin
          teilchen[length(teilchen)-1].nMin:=exprToFloat(false,s,kvs,nil);
          kvs.add('nmin'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].nMin);
          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;

  if dPX<0 then
    dPX:=deltaX;
  if dPY<0 then
    dPY:=sqrt(deltaX);
  if dPZ<0 then
    dPZ:=sqrt(deltaX);

  pro:=tProtokollant.create(prot,'create');
  case zeitverfahren of
    zfEulerVorwaerts:
      pro.schreibe('Iteration mittels Euler-Vorwärts');
    zfRungeKuttaDreiAchtel:
      pro.schreibe('Iteration mittels Runge-Kutta-3/8');
    zfRungeKuttaVier:
      pro.schreibe('Iteration mittels Runge-Kutta-4');
    zfRungeKuttaZehn:
      pro.schreibe('Iteration mittels Runge-Kutta-10');
    zfRungeKuttaZwoelf:
      pro.schreibe('Iteration mittels Runge-Kutta-12');
    zfRungeKuttaVierzehn:
      pro.schreibe('Iteration mittels Runge-Kutta-14');
    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('dPX  =  '+floattostr(dPX));
  pro.schreibe('pXMax = '+floattostr(pXMax));
  pro.schreibe('dPY  =  '+floattostr(dPY));
  pro.schreibe('pYMax = '+floattostr(pYMax));
  pro.schreibe('dPZ  =  '+floattostr(dPZ));
  pro.schreibe('pZMax = '+floattostr(pZMax));
  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;

  if length(sighupSimulationen)=0 then begin
    new(na);
    na^.sa_Handler := sigActionHandler(@signalCapture);
    fillchar(na^.sa_Mask,sizeof(na^.sa_mask),#0);
    na^.sa_Flags := 0;
    {$ifdef Linux}               // Linux specific
    na^.sa_Restorer := Nil;
    {$endif}
    if (fPSigaction(SIGHUP, na, nil) <> 0) then begin
      writeln(stderr, 'Error: ', fpgeterrno, '.');
      halt(1);
    end;
    dispose(na);
  end;
  setlength(sighupSimulationen,length(sighupSimulationen)+1);
  sighupSimulationen[length(sighupSimulationen)-1]:=self;

  gitter:=tGitter.create(self,round(breite/deltaX)+1,deltaX,kvs,teilchen,lichter,zeitverfahren,prot,'gitter',2*round(pXMax/dPX)+1,2*round(pYMax/dPY)+1,2*round(pZMax/dPZ)+1,dPX,dPY,dPZ);
  for i:=0 to length(teilchen)-1 do
    teilchen[i].free;
  setlength(teilchen,0);
  lichter.free;
end;

destructor tSimulation.destroy;
var
  i,j: longint;
begin
  for i:=0 to length(ausgabeDateien)-1 do
    ausgabeDateien[i].free;
  setlength(ausgabeDateien,0);
  gitter.free;
  prot.free;
  for i:=length(sighupSimulationen)-1 downto 0 do
    if sighupSimulationen[i]=self then begin
      for j:=i+1 to length(sighupSimulationen)-1 do
        sighupSimulationen[j-1]:=sighupSimulationen[j];
      setlength(sighupSimulationen,length(sighupSimulationen)-1);
    end;
  inherited destroy;
end;

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

  if gitter.abbruch or (dT>sDT) then
    dT:=sDT;

  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
    for i:=0 to length(ausgabeDateien)-1 do
      ausgabeDateien[i].speichereWerte(gitter,sT,sDX);
    sT:=sT+sDT;
  end;
  zeitDatei:=zeitDatei+now;

  if gotSighup or
     (fortschrittsAnzeige and (floor(100*Gitter.t/tEnde) < floor(100*(Gitter.t+dT)/tEnde))) then begin
    gotSighup:=false;
    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);
  end;

  result:=gitter.t<tEnde;
end;

// allgemeine Funktionen *******************************************************

procedure SignalCapture(signal : longint); cdecl;
var
  i: longint;
begin
  if signal=SIGHUP then
    for i:=0 to length(sighupSimulationen)-1 do
      sighupSimulationen[i].gotSighup:=true;
end;

initialization
  setlength(sighupSimulationen,0);
finalization
  setlength(sighupSimulationen,0);
end.