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;
  tEMFeldInhalt = (
    efA,efAX,efAY,efAZ,
    efDAXDT,efDAYDT,efDAZDT,
    efDPhiDX
  );
  tMaterieFeldInhalt = (
    mfN,mfDPsiDX,       // Teilchendichte, Geschwindigkeitsdichte
    mfP,mfPX,mfPY,mfPZ, // Impuls
    mfGamma,mfIGamma    // rel. Gamma-Faktor
  );

const
  erstesEMFMitAbleitung  =  efAX;
  letztesEMFMitAbleitung =  efDAZDT;
  erstesMatFMitAbleitung  = mfN;
  letztesMatFMitAbleitung = mfDPsiDX;

type

  tGitter = class;
  tFelder = class;
  tSimulation = 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;
    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;

  { tWertePunkt }

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

    constructor create(linkerNachbar: tWertePunkt; derBesitzer: tFelder);
    destructor destroy; override;
    procedure berechneGammaUndP;
    procedure berechneDPhiDX(dX: extended);
    procedure berechneNAbleitung(iDX,pDNMax: extended); overload;
    procedure berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean); overload;
    procedure berechneAbleitungen(iDX: extended);
    procedure liKo(in1,in2: tWertePunkt; fak2: extended); overload;                            // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
    procedure liKo(in1,in2,in3: tWertePunkt; fak2,fak3: extended); overload;
    procedure liKo(in1,in2,in3,in4: tWertePunkt; fak2,fak3,fak4: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6: tWertePunkt; fak2,fak3,fak4,fak5,fak6: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31,in32: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31,fak32: extended); overload;
    function maxDT: extended;
    function dump(mitNachbarn: boolean): string;
    procedure nichtnegativieren;
  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 tWertePunkt;
    par:    tGitter;

    constructor create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
    destructor destroy; override;
    procedure berechneAbleitungen(dX,iDX,pDNMax: extended);
    procedure liKo(in1,in2: tFelder; fak2: extended); overload;                            // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
    procedure liKo(in1,in2,in3: tFelder; fak2,fak3: extended); overload;
    procedure liKo(in1,in2,in3,in4: tFelder; fak2,fak3,fak4: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6: tFelder; fak2,fak3,fak4,fak5,fak6: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7: tFelder; fak2,fak3,fak4,fak5,fak6,fak7: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31: extended); overload;
    procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31,in32: tFelder; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31,fak32: extended); overload;
    function maxDT: extended;
  end;

  { tGitter }

  tGitter = class(tObject)
  private
    besitzer:                     tSimulation;
    prot:                         tProtokollant;
    zeitverfahren:                tZeitverfahren;
    pDNMaxDynamisch,abbruch:      boolean;
    dX,iDX,pDNMax,dTMaximum,xl,t: extended;
    kvs:                          tKnownValues;
    procedure diffusionsTermAnpassen(pro: tProtokollant);
    function pruefeMaxDT(aF: longint; var dTMax,dT: extended; dTMin: extended): boolean; // wurde verkleinert?
    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; 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; dTMin: 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,dTMin: 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[tEMFeldInhalt] of string = (
    'A','AX','AY','AZ','DAXDT','DAYDT','DAZDT','DPHIDX'
  );
  matFeldNamen: array[tMaterieFeldInhalt] of string = (
    'N','DPSIDX','P','PX','PY','PZ','GAMMA','IGAMMA'
  );
  minAusgabeBuffer = 1024*1024;

var
  sighupSimulationen: array of tSimulation;

{ 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;

  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<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;
  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: 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 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<0 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].emWerte[emInhalt,ableitung],(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);
        move(gitter.felders[gitter.aktuelleFelder].inhalt[i].matWerte[teilchen,matInhalt,ableitung],(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;

{ tWertePunkt }

constructor tWertePunkt.create(linkerNachbar: tWertePunkt; derBesitzer: tFelder);
var
  emF: tEMFeldInhalt;
  maF: tMaterieFeldInhalt;
  abl: boolean;
  i:   longint;
begin
  inherited create;
  besitzer:=derBesitzer;
  fillchar(matWerte,sizeof(matWerte),#0);
  setlength(matWerte,length(besitzer.spezLadungen));
  for i:=0 to length(matWerte)-1 do
    for maF:=low(tMaterieFeldInhalt) to high(tMaterieFeldInhalt) do
      for abl:=false to true do
        matWerte[i,maF,abl]:=0;
  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]:= besitzer.spezLadungen[i] * emWerte[efAX,false] + matWerte[i,mfDPsiDX,false];
    matWerte[i,mfPY,false]:= besitzer.spezLadungen[i] * emWerte[efAY,false];
    matWerte[i,mfPZ,false]:= besitzer.spezLadungen[i] * 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.berechneDPhiDX(dX: extended);
var
  i: longint;
begin
  // d2Phi/dx2 = rho  ... hier muss integriert werden:
  // dPhi/dx = dPhi/dx(x- deltax) + <rho> * deltax

  if assigned(lN) then begin
    emWerte[efDPhiDX,false]:=lN.emWerte[efDPhiDX,false];
    for i:=0 to length(matWerte)-1 do
      emWerte[efDPhiDX,false]:=
        emWerte[efDPhiDX,false]
        + besitzer.spezLadungen[i] * (lN.matWerte[i,mfN,false] + matWerte[i,mfN,false]) * dX/2;
  end
  else begin
    emWerte[efDPhiDX,false]:=0;
    for i:=0 to length(matWerte)-1 do
      emWerte[efDPhiDX,false]:=
        emWerte[efDPhiDX,false]
        + besitzer.spezLadungen[i] * matWerte[i,mfN,false] * dX/2;
  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] * iDX)
      + lN.matWerte[i,mfN,false] * lN.matWerte[i,mfIGamma,false] * (pDNMax + lN.matWerte[i,mfPX,false] * iDX)
      -    matWerte[i,mfN,false]  *   matWerte[i,mfIGamma,false] * 2 * pDNMax) / 2;
    // Der zweite Summand entspricht (in etwa) einer thermischen Verschmierung
    // und soll die numerische Stabilität bis zu p * d(ln n)/dx von pDNMax gewährleisten.
    // Es handelt sich um einen Diffusionsterm dn/dt = alpha * Laplace n mit
    // alpha = pDNMax * dX^2
end;

procedure tWertePunkt.berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean);
var
  i: longint;
begin
  if rechts then begin
    for i:=0 to length(matWerte)-1 do
      // rechter Randterm von dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n
      matWerte[i,mfN,true]:=
        ( lN.matWerte[i,mfN,false] * lN.matWerte[i,mfIGamma,false] * (pDNMax + lN.matWerte[i,mfPX,false] * iDX)
        -    matWerte[i,mfN,false]  *   matWerte[i,mfIGamma,false] * (pDNMax -    matWerte[i,mfPX,false] * iDX)) / 2;
  end
  else begin
    for i:=0 to length(matWerte)-1 do
      // linker Randterm von 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] * iDX)
        -    matWerte[i,mfN,false]  *   matWerte[i,mfIGamma,false] * (pDNMax +    matWerte[i,mfPX,false] * iDX)) / 2;
  end
end;

procedure tWertePunkt.berechneAbleitungen(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]:=
      besitzer.spezLadungen[i] * emWerte[efDPhiDX,false]
        - (rN.matWerte[i,mfGamma,false] - lN.matWerte[i,mfGamma,false]) * iDX/2;

  // d2A/dt2 = Laplace(A) ...    ( - Nabla(dPhi/dt) wird auf dA/dt aufgeschlagen !!! )
  emWerte[efDAXDT,true]:=
    ( rN.emWerte[efAX,false] - 2*emWerte[efAX,false] + lN.emWerte[efAX,false] )*sqr(iDX);
  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);
  // ... - Rho/gamma p
  for i:=0 to length(matWerte)-1 do begin
    emWerte[efDAXDT,true]:=
      emWerte[efDAXDT,true]
      - (besitzer.spezLadungen[i] * matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPX,false]);
    emWerte[efDAYDT,true]:=
      emWerte[efDAYDT,true]
      - (besitzer.spezLadungen[i] * matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPY,false]);
    emWerte[efDAZDT,true]:=
      emWerte[efDAZDT,true]
      - (besitzer.spezLadungen[i] * matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * matWerte[i,mfPZ,false]);
  end;

  // dA/dt = dA/dt - Nabla(Phi)
  emWerte[efAX,true]:=emWerte[efDAXDT,false] - emWerte[efDPhiDX,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;

function tWertePunkt.maxDT: extended;
var
  i: longint;
begin
  result:=infinity;
  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;

function tWertePunkt.dump(mitNachbarn: boolean): string;
var
  i: longint;
begin
  if mitNachbarn and assigned(lN) then result:=lN.dump(false)+' << '
  else result:='';

  for i:=0 to length(matWerte)-1 do
    result:=
      result+
      'N['+  inttostr(i)+'] = '+floattostr(matWerte[i,mfN,false])+' '+
      'N''['+inttostr(i)+'] = '+floattostr(matWerte[i,mfN,true])+' ';

  result:=result+'maxDT: '+floattostr(maxDT);

  if mitNachbarn and assigned(rN) then result:=result+' >> '+rN.dump(false);
end;

procedure tWertePunkt.nichtnegativieren; // Dichten nicht negativ machen
var
  i:       longint;
  defizit: extended;
  li,re:   tWertePunkt;
  pro:     tProtokollant;
begin
  for i:=0 to length(matWerte)-1 do
    if matWerte[i,mfN,false]<0 then begin
      defizit:=-matWerte[i,mfN,false];
      matWerte[i,mfN,false]:=0;
      li:=self.lN;
      re:=self.rN;
      repeat
        if assigned(li) then begin
          if li.matWerte[i,mfN,false]>0 then begin
            if defizit>=li.matWerte[i,mfN,false] then begin
              defizit:=defizit-li.matWerte[i,mfN,false];
              li.matWerte[i,mfN,false]:=0;
            end
            else begin
              li.matWerte[i,mfN,false]:=li.matWerte[i,mfN,false]-defizit;
              defizit:=0;;
            end;
          end;
          li:=li.lN;
        end;
        if assigned(re) then begin
          if re.matWerte[i,mfN,false]>0 then begin
            if defizit>=re.matWerte[i,mfN,false] then begin
              defizit:=defizit-re.matWerte[i,mfN,false];
              re.matWerte[i,mfN,false]:=0;
            end
            else begin
              re.matWerte[i,mfN,false]:=re.matWerte[i,mfN,false]-defizit;
              defizit:=0;;
            end;
          end;
          re:=re.rN;
        end;
      until (defizit<=0) or not (assigned(li) or assigned(re));
      if defizit>0 then begin
        pro:=tProtokollant.create(besitzer.par.prot,'WertePunkt.nichtnegativieren');
        pro.schreibe('Ich konnte ein Teilchendefizit nicht auffüllen!',true);
        besitzer.par.abbrechen;
      end;
    end;
end;

{ linearkombinationsspezifische Methoden von tWertePunkt und tFelder }

{$INCLUDE linearkombination.inc}
{$DEFINE lkA3}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA4}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA5}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA6}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA7}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA8}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA9}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA10}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA11}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA12}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA14}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA15}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA16}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA17}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA18}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA19}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA22}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA23}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA31}
{$INCLUDE linearkombination.inc}
{$DEFINE lkA32}
{$INCLUDE linearkombination.inc}

{ 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(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]:=tWertePunkt.create(nil,self);
  for i:=1 to length(inhalt)-1 do
    inhalt[i]:=tWertePunkt.create(inhalt[i-1],self);
  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
      inhalt[i].matWerte[j,mfN,false]:=teilchen[j].gibDichte(par.xl+(i-2)*par.dX,parent.kvs);
      inhalt[i].matWerte[j,mfDPsiDX,false]:=0;
      massen[j]:=massen[j]+inhalt[i].matWerte[j,mfN,false]*par.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:    tEMFeldInhalt;
  rechts: boolean;
  i:      longint;
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
      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);
end;

procedure tFelder.berechneAbleitungen(dX,iDX,pDNMax: extended);
var
  i: longint;
begin
  for i:=0 to length(inhalt)-1 do
    inhalt[i].berechneGammaUndP;
  for i:=0 to length(inhalt)-1 do
    inhalt[i].berechneDPhiDX(dX);
  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;

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

  inhalt[1].berechneNAbleitung(iDX,pDNMax,false);
  inhalt[length(inhalt)-2].berechneNAbleitung(iDX,pDNMax,true);

  setzeRaender(iDX);

  for i:=1 to length(inhalt)-2 do
    inhalt[i].berechneAbleitungen(iDX);
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(derBesitzer: tSimulation; 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;
  abbruch:=false;
  besitzer:=derBesitzer;
  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);
    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);

  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.diffusionsTermAnpassen(pro: tProtokollant);
var
  i,j,k:      longint;
  curMax,tmp: extended;
begin
  curMax:=0;
  k:=0;
  for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do
    for j:=0 to felders[aktuelleFelder].matAnz-1 do begin
      tmp:=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-100));
      if tmp>curMax then begin
        curMax:=tmp;
        k:=i;
      end;
    end;
  if curMax > pDNMax then begin
    if pDNMaxDynamisch then begin
      pDNMax:=2*curMax;
      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, x='+floattostr(xL+k*dX)+'), die numerische Stabilität ist nicht mehr gewährleistet, Simulation wird abgebrochen!',true);
      pro.schreibe(' Lösung: größeren Diffusionsterm wählen',true);
      pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern',true);
      abbrechen;
    end;
  end;
end;

function tGitter.pruefeMaxDT(aF: longint; var dTMax,dT: extended; dTMin: extended): boolean; // wurde verkleinert?
begin
  dTMax:=min(dTMax,felders[aF].maxDT);
    // das maximale dT, welches bei den momentanen Ableitungen nicht zu
    // unphysikalischen Effekten (z.B. negativen Dichten) führt

  result:=(dTMax<dT) and (dT>1.1*dTMin);
  if result then begin // beabsichtigter Zeitschritt ist zu groß,
    dT:=dTMax/4;       // dann machen wir ihn doch kleiner!
    if dT<dTMin then
      dT:=dTMin
    {$IFDEF Zeitschrittueberwachung}
    else begin
      prot.schreibe('pruefeMaxDT: dT -> '+floattostr(dT)+' (t = '+floattostr(t)+')');
      prot.spuelen;
    end
    {$ENDIF};
  end;
end;

procedure tGitter.abbrechen;
var
  i,j,k: longint;
  maF:   tMaterieFeldInhalt;
  emF:   tEMFeldInhalt;
  abl:   boolean;
begin
  abbruch:=true;
  for i:=0 to length(felders)-1 do
    for j:=0 to length(felders[i].inhalt)-1 do begin
      for emF:=low(tEMFeldInhalt) to high(tEMFeldInhalt) do
        for abl:=false to true do
          felders[i].inhalt[j].emWerte[emF,false]:=0;
      for k:=0 to felders[i].matAnz-1 do
        for maF:=low(tMaterieFeldInhalt) to high(tMaterieFeldInhalt) do
          for abl:=false to true do
            felders[i].inhalt[j].matWerte[k,maF,false]:=0;
    end;
end;

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

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

    if pruefeMaxDT(aktuelleFelder,dTMax,dT,dTMin) then
      continue;

    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;
  dT:=max(dT,dTMax/100);
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 (not abbruch) and (abs(ns[i]-felders[aktuelleFelder].massen[i])>0.5) then begin
      result:=false;
      pro.schreibe(' n['+inttostr(i)+'] ('+floattostr(ns[i])+') unterscheidet sich stark vom Anfangswert ('+floattostr(felders[aktuelleFelder].massen[i])+'), es scheinen sehr viele Teilchen entstanden oder verloren gegangen 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.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;
  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;
  pDNMax:=-1;
  fortschrittsAnzeige:=false;
  dTMin:=-1;
  gotSighup:=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-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('minZeitschritt ',s) then begin
          dTMin:=exprtofloat(false,s,kvs,nil);
          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('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 dTMin<0 then
    dTMin:=min(1E-30,dT*1E-5);

  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('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;

  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,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;
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,dTMin);
  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);
    if errorCode<2 then
      if not gitter.dumpErhaltungsgroessen then begin
        errorCode:=2;
        exit;
      end;
  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.