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, fftw;

type
  tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaDreiAchtel,zfRungeKuttaVier,zfRungeKuttaZehn,zfRungeKuttaZwoelf,zfRungeKuttaVierzehn);
  tVerteilungsfunktion = function(x: double): double;
  tEMFeldGroesze = (
    efAX,efAY,
    efDAXDT,efDAYDT,
    efEX,efBZ
  );
  tEMQuellGroesze = (
    eqRho,eqJX,eqJY
  );
  tMaterieFeldGroesze = (mfPY);
  tMaterieSpeicherGroesze = (msPX,msPXSqr,msPY,msVX,msVY,msN,msPXRipple);
  tTeilchenSpeziesGroeszen = (tsgMasse,tsgIQdrMasse,tsgLadung);
  tEMQuellen = array[tEMQuellGroesze] of double;

const
  erstesEMFMitAbleitung  =  efAX;
  letztesEMFMitAbleitung =  efDAYDT;

type

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

  { tAusgabeDatei }

  tAusgabeDatei = class
  private
    ableitung:          boolean;
    emFeld:             tEMFeldGroesze;
    emQuelle:           tEMQuellGroesze;
    matFeld:            tMaterieSpeicherGroesze;
    wasSpeichern,       // 0: EM-Feld; 1: EM-Quelle; 2: Mat-Feld
    teilchen,           // -1: Gesamt-Feld; sonst: Feld des 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:                double;
    constructor create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string);
    destructor destroy; override;
    function dump: string;
    procedure schreibeKopf;
    procedure speichereWerte(gitter: tGitter; sT, sDX: double);
  end;

  { tTeilchenSpezies }

  tTeilchenSpezies = class
  private
    ortsGrenzen:     array of double;
    dichteStuecke:   array of string;
  public
    nMin,nMax,breite: double;                                    // minimale und maximale Massendichte in nc; Breite der Impulsverteilung
    eigenschaften:    array[tTeilchenSpeziesGroeszen] of double; // q in e; m in me
    diffusion:        array[0..1] of double;                     // Diffusionsterm in x- und p-Richtung (= maximales d(ln n)/d[xp] * iD[XP])
    constructor create; overload;
    constructor create(original: tTeilchenSpezies); overload;
    destructor destroy; override;
    procedure dump(prot: tProtokollant; prefix: string);
    function gibDichte(x: double; kvs: tKnownValues): double; // Massendichte
    procedure liesDichteFunktion(rest: string; ifile: tMyStringList; kvs: tKnownValues; teilchen: array of tTeilchenSpezies; prot: tProtokollant);
  end;

  { TFelder }
(*
  procedure berechneAbleitungen(pX: double); inline;
  function nichtnegativieren: double; inline;
*)

  tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen
  private

    // 2d-Arrays (x & px) sind wie folgt indiziert: [iP + iX*aP]

    teilchen:           array of tTeilchenSpezies;
    lichters:           array[boolean] of tMyStringlist;
    massen:             array of double;                         // Gesamtmassen je Teilchenspezies zur Überwachung der Erhaltungsgrößen
    gesamtDefizit:      double;
    impulsRaumGradient: array of array[0..1] of pDouble;         // Ableitung des Impulsraumes nach x und px
    iMGamma:            array of pDouble;                        // 1/m/gamma (= v/p)
    _aX,_aP:            longint;
    _dX,_dP,iDX:        double;
    ffts:               array of array[0..1,boolean] of fftw_plan_double; // fft-Algorithmen
    fftTmp,lFftTmp:     Pcomplex_double;                         // Zwischenspeicher für ffts
    procedure initialisiereDichte(ort,tlc: longint; breite,n: double); inline;
    procedure pruefeArrayEnden(fehler: string);
    procedure nichtnegativieren;
  public
    emFelder:  array[tEMFeldGroesze,boolean] of pDouble;         // EM-Felder und deren Zeitableitungen
                                                                   // dPhiDX und B[xyz] haben keine sinnvolle Ableitung hier!
    emQuellen: array[tEMQuellGroesze] of pDouble;                // Quellen für die EM-Felder
    matFelder: array of array[tMaterieFeldGroesze] of pDouble;   // px-unabhängige Felder, nach Teilchenspezies sortiert
    impulsraum: array of array[boolean] of pDouble;              // px-abhängige Felder, nach Teilchenspezies sortiert, sowie deren Ableitung

    gitter: tGitter;

    constructor create(anzX,anzP: longint; deltaX,deltaP: double; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
    destructor destroy; override;
    procedure berechneAbleitungen;
    procedure setzeNull; inline;
    procedure dumpErhaltungsgroeszen;
    property aX: longint read _aX;
    property aP: longint read _aP;
    property dX: double read _dX;
    property dP: double read _dP;
    function impulsIntegral(ort,tlc: longint; emQ: tEMQuellGroesze): double; overload; inline;
    function impulsIntegral(ort,tlc: longint; maF: tMaterieSpeicherGroesze): double; overload; inline;
    {$DEFINE LiKoInterface}
    {$INCLUDE linearkombinationen.inc}
    {$UNDEF LiKoInterface}
    procedure schreibePhasenraum(tlc: longint; datNam: string);
  end;

  { tGitter }

  tGitter = class(tObject)
  private
    simulation:                   tSimulation;
    prot:                         tProtokollant;
    zeitverfahren:                tZeitverfahren;
    abbruch:                      boolean;
    dX,iDX,xl,t:                  double;
    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(besitzer: tSimulation; aX,aP: longint; deltaX,deltaP: double; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; name: string);
    destructor destroy; override;
    procedure iteriereSchritt(dT: double);
    procedure dumpErhaltungsgroeszen;
  end;

  { tSimulation }

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

procedure SignalCapture(signal : longint); cdecl;

implementation

const
  emFeldNamen: array[tEMFeldGroesze] of string = (
    'AX','AY','DAXDT','DAYDT','EX','BZ'
  );
  emQuellNamen: array[tEMQuellGroesze] of string = (
    'RHO','JX','JY'
  );
  matSpeicherNamen: array[tMaterieSpeicherGroesze] of string = (
    'PX','PXSQR','PY','VX','VY','N','PXRIPPLE'
  );
  minAusgabeBuffer = 1024*1024;

var
  simulationen: array of tSimulation;

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

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

  emFeld:=efAX;
  emQuelle:=eqRho;
  matFeld:=msN;
  wasSpeichern:=0;
  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
          emFeld:=emF;
          wasSpeichern:=0;
          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;
        emQuelle:=emQ;
        wasSpeichern:=1;
        gef:=true;
        break;
      end;
  if not gef then
    for maF:=low(tMaterieSpeicherGroesze) to high(tMaterieSpeicherGroesze) do
      if matSpeicherNamen[maF] = feldName then begin
        ableitung:=false;
        matFeld:=maF;
        wasSpeichern:=2;
        gef:=true;
        break;
      end;

  if not gef then begin
    pro.schreibe('tAusgabeDatei.create kennt Größe '''+feldName+''' nicht!',true);
    nam:='';
    for abl:=false to true do
      for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do
        nam:=nam+', '+copy('D',1,byte(abl))+emFeldNamen[emF];
    for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do
      nam:=nam+', '+emQuellNamen[emQ];
    for maF:=low(tMaterieSpeicherGroesze) to high(tMaterieSpeicherGroesze) do
      nam:=nam+', '+matSpeicherNamen[maF];
    delete(nam,1,2);
    pro.schreibe('Ich kenne nur: '+nam,true);
    raise exception.create('Fehler in tAusgabeDatei.create!');
  end;

  case wasSpeichern of
    0: pre:=emFeldNamen[emFeld];
    1: pre:=emQuellNamen[emQuelle];
    2: pre:=matSpeicherNamen[matFeld];
  end;
  if teilchen>=0 then
    pre:=pre+inttostr(teilchen+1);

  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);
    raise exception.create('Fehler in tAusgabeDatei.destroy!');
  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
  case wasSpeichern of
    0: result:=emFeldNamen[emFeld];
    1: result:=emQuellNamen[emQuelle];
    2: result:=matSpeicherNamen[matFeld];
  end;
  if teilchen>=0 then
    result:=result+'['+inttostr(teilchen+1)+']';
  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);
    raise exception.create('Fehler in tAusgabeDatei.schreibeKopf!');
  end;
  if bufPos<>0 then
    schreibeBuffer(0);
  tCnt:=0;
  assignfile(datei,pre+'-'+inttostr(nNum)+suf);
  inc(nNum);
  rewrite(datei,1);
  blockwrite(datei,nNum,sizeof(longint));
  blockwrite(datei,zeitAnz,sizeof(longint));
  closefile(datei);
end;

procedure tAusgabeDatei.speichereWerte(gitter: tGitter; sT, sDX: double);
var i,cnt:     longint;
    sX,cX,val: double;
begin
  if (teilchen>=length(gitter.felders[gitter.aktuelleFelder].teilchen)) then begin
    pro.schreibe('Teilchen '+inttostr(teilchen+1)+' gibt es nicht, da kann ich auch nichts speichern!',true);
    raise exception.create('Fehler in tAusgabeDatei.speichereWerte!');
  end;

  if sT>=nNum+sDT/2 then
    schreibeKopf;

  cnt:=floor(((gitter.felders[gitter.aktuelleFelder].aX-1)*gitter.dX+Min(gitter.dX,sDX)/2)/sDX+1);
  cX:=gitter.xl;
  sX:=cX+(cnt-1)*sDX;

  inc(tCnt);
  if bufLen-bufPos < sizeof(longint) + (2+cnt)*sizeof(double) then
    schreibeBuffer(sizeof(longint) + (2+cnt)*sizeof(double));

  if bufLen-bufPos < sizeof(longint) + (2+cnt)*sizeof(double) then begin
    pro.schreibe('Schreibbuffer ist zu klein ('+inttostr(bufLen)+' statt mindestens '+inttostr(sizeof(longint) + (2+cnt)*sizeof(double))+' Bytes)!');
    raise exception.create('Fehler in tAusgabeDatei.speichereWerte!');
  end;

  move(cX,(buf+bufPos)^,sizeof(double));
  inc(bufPos,sizeof(double));
  move(sX,(buf+bufPos)^,sizeof(double));
  inc(bufPos,sizeof(double));
  move(cnt,(buf+bufPos)^,sizeof(longint));
  inc(bufPos,sizeof(longint));

  sX:=cX-Min(gitter.dX,sDX)/2;
  case wasSpeichern of
    0: // em-Feld speichern
      for i:=0 to gitter.felders[gitter.aktuelleFelder].aX-1 do begin
        while cX>=sX do begin
          dec(cnt);
          move((gitter.felders[gitter.aktuelleFelder].emFelder[emFeld,ableitung]+i)^,(buf+bufPos)^,sizeof(double));
          inc(bufPos,sizeof(double));
          sX:=sX+sDX;
        end;
        cX:=cX+gitter.dX;
      end;
    1:  // em-Quelle
      for i:=0 to gitter.felders[gitter.aktuelleFelder].aX-1 do begin
        while cX>=sX do begin
          dec(cnt);
          val:=gitter.felders[gitter.aktuelleFelder].impulsIntegral(i,teilchen,emQuelle);
          move(val,(buf+bufPos)^,sizeof(double));
          inc(bufPos,sizeof(double));
          sX:=sX+sDX;
        end;
        cX:=cX+gitter.dX;
      end;
    2:  // Materiefeld
      for i:=0 to gitter.felders[gitter.aktuelleFelder].aX-1 do begin
        while cX>=sX do begin
          dec(cnt);
          val:=gitter.felders[gitter.aktuelleFelder].impulsIntegral(i,teilchen,matFeld);
          move(val,(buf+bufPos)^,sizeof(double));
          inc(bufPos,sizeof(double));
          sX:=sX+sDX;
        end;
        cX:=cX+gitter.dX;
      end;
  end{of case};
  if cnt<>0 then begin
    pro.schreibe('Falsche Anzahl an Ortsschritten geschrieben ('+inttostr(cnt)+')!',true);
    raise exception.create('Fehler in tAusgabeDatei.speichereWerte!');
  end;
end;

// tTeilchenSpezies ************************************************************

constructor tTeilchenSpezies.create;
var
  egr: tTeilchenSpeziesGroeszen;
begin
  inherited create;
  fillchar(ortsGrenzen,sizeof(ortsGrenzen),#0);
  setlength(ortsGrenzen,0);
  fillchar(dichteStuecke,sizeof(dichteStuecke),#0);
  setlength(dichteStuecke,0);
  nMin:=1E-3;
  nMax:=nMin;
  for egr:=low(tTeilchenSpeziesGroeszen) to high(tTeilchenSpeziesGroeszen) do
    eigenschaften[egr]:=byte(egr in [tsgMasse,tsgIQdrMasse]);
  diffusion[0]:=0;
  diffusion[1]:=0;
  breite:=1;
end;

constructor tTeilchenSpezies.create(original: tTeilchenSpezies);
var
  egr: tTeilchenSpeziesGroeszen;
  i:   longint;
begin
  create;
  nMin:=original.nMin;
  nMax:=original.nMax;
  for egr:=low(tTeilchenSpeziesGroeszen) to high(tTeilchenSpeziesGroeszen) do
    eigenschaften[egr]:=original.eigenschaften[egr];
  for i:=0 to 1 do
    diffusion[i]:=original.diffusion[i];
  breite:=original.breite;
  fillchar(ortsGrenzen,sizeof(ortsGrenzen),#0);
  setlength(ortsGrenzen,length(original.ortsGrenzen));
  for i:=0 to length(ortsGrenzen)-1 do
    ortsGrenzen[i]:=original.ortsGrenzen[i];
  fillchar(dichteStuecke,sizeof(dichteStuecke),#0);
  setlength(dichteStuecke,length(original.dichteStuecke));
  for i:=0 to length(dichteStuecke)-1 do
    dichteStuecke[i]:=original.dichteStuecke[i];
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+'m    =    '+floattostr(eigenschaften[tsgMasse])+' me');
  prot.schreibe(prefix+'q    =    '+floattostr(eigenschaften[tsgLadung])+' e');
  prot.schreibe(prefix+'dLnN/dX < '+floattostr(diffusion[0])+' / dX');
  prot.schreibe(prefix+'dLnN/dP < '+floattostr(diffusion[1])+' / dP');
  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: double; kvs: tKnownValues): double;
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);
        raise exception.create('Fehler in tTeilchenSpezies.liesDichteFunktion!');
      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);
      raise exception.create('Fehler in tTeilchenSpezies.liesDichteFunktion!');
    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;

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

constructor tFelder.create(anzX,anzP: longint; deltaX,deltaP: double; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
var
  i,j:        longint;
  rechts,abl: boolean;
  dens:       double;
  emF:        tEMFeldGroesze;
  emQ:        tEMQuellGroesze;
  maF:        tMaterieFeldGroesze;
begin
  inherited create;
  gitter:=parent;
  gesamtDefizit:=0;
  _aX:=anzX;
  _aP:=anzP;
  _dX:=deltaX;
  _dP:=deltaP;
  iDX:=1/dX;
  fillchar(teilchen,sizeof(teilchen),#0);
  setlength(teilchen,length(_teilchen));
  for i:=0 to length(teilchen)-1 do
    teilchen[i]:=tTeilchenSpezies.create(_teilchen[i]);

  fillchar(impulsRaumGradient,sizeof(impulsRaumGradient),#0);
  setlength(impulsRaumGradient,length(teilchen));
  for i:=0 to length(impulsRaumGradient)-1 do
    for j:=0 to length(impulsRaumGradient[i])-1 do begin
      fftw_getmem(impulsRaumGradient[i,j],(aP*aX+1)*sizeof(double));
      fillchar(impulsRaumGradient[i,j]^,aP*aX*sizeof(double),#0);
      (impulsRaumGradient[i,j]+aP*aX)^:=pi;
    end;
  fillchar(iMGamma,sizeof(iMGamma),#0);
  setlength(iMGamma,length(teilchen));
  for i:=0 to length(iMGamma)-1 do begin
    fftw_getmem(iMGamma[i],aP*aX*sizeof(double));
    fillchar(iMGamma[i]^,aP*aX*sizeof(double),#0);
  end;
  for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do begin
    emFelder[emF,true]:=nil;
    for abl:=false to not(emF in [efEX,efBZ]) do begin
      fftw_getmem(emFelder[emF,abl],aX*sizeof(double));
      fillchar(emFelder[emF,abl]^,aX*sizeof(double),#0);
    end;
  end;
  for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do begin
    fftw_getmem(emQuellen[emQ],aX*sizeof(double));
    fillchar(emQuellen[emQ]^,aX*sizeof(double),#0);
  end;
  fillchar(matFelder,sizeof(matFelder),#0);
  setlength(matFelder,length(teilchen));
  for i:=0 to length(matFelder)-1 do
    for maF:=low(tMaterieFeldGroesze) to high(tMaterieFeldGroesze) do begin
      fftw_getmem(matFelder[i,maF],aX*sizeof(double));
      fillchar(matFelder[i,maF]^,aX*sizeof(double),#0);
    end;
  fillchar(impulsraum,sizeof(impulsraum),#0);
  setlength(impulsraum,length(teilchen));
  for i:=0 to length(impulsraum)-1 do
    for abl:=false to true do begin
      fftw_getmem(impulsraum[i,abl],(aX*aP+1)*sizeof(double));
      fillchar(impulsraum[i,abl]^,aX*aP*sizeof(double),#0);
      (impulsraum[i,abl]+aX*aP)^:=pi;
    end;

  i:=max((aP div 2 + 1)*aX,(aX div 2 + 1)*aP);
  fftw_getmem(fftTmp,(i+1)*sizeof(complex_double));
  fillchar(fftTmp^,(i+1)*sizeof(complex_double),#0);
  lFftTmp:=fftTmp+i;

  lFftTmp^.re:=pi;
  lFftTmp^.im:=sqr(pi);

  fillchar(ffts,sizeof(ffts),#0);     // Vorsicht, die Planung der FFTs überschreibt den Input, auch wenn die eigentlichen FFTs das nicht tun!
  setlength(ffts,length(teilchen));
  for i:=0 to length(ffts)-1 do begin
    ffts[i,0,false]:=  // Planung der Hintransformationen über x
      fftw_plan_many_dft_r2c(1,@_aX,aP,impulsraum[i,false],nil,aP,1,fftTmp,nil,aP,1,[fftw_preserve_input,fftw_exhaustive]);
    ffts[i,0,true]:=  // Planung der Rücktransformationen über x
      fftw_plan_many_dft_c2r(1,@_aX,aP,fftTmp,nil,aP,1,impulsRaumGradient[i,0],nil,aP,1,[fftw_destroy_input,fftw_exhaustive]);

    ffts[i,1,false]:=  // Planung der Hintransformationen über p
      fftw_plan_many_dft_r2c(1,@_aP,aX,impulsraum[i,false],nil,1,aP,fftTmp,nil,1,aP div 2 + 1,[fftw_preserve_input,fftw_exhaustive]);
    ffts[i,1,true]:=  // Planung der Rücktransformationen über p
      fftw_plan_many_dft_c2r(1,@_aP,aX,fftTmp,nil,1,aP div 2 + 1,impulsRaumGradient[i,1],nil,1,aP,[fftw_destroy_input,fftw_exhaustive]);
  end;

  pruefeArrayEnden('FFT schreibt beim Planen über das Ende hinaus!');

  (*
     fftw_plan_many_dft_r2c(int rank, const int *n, int howmany,
                                      double *in, const int *inembed,
                                      int istride, int idist,
                                      fftw_complex *out, const int *onembed,
                                      int ostride, int odist,
                                      unsigned flags);
     fftw_plan fftw_plan_many_dft_c2r(int rank, const int *n, int howmany,
                                      fftw_complex *in, const int *inembed,
                                      int istride, int idist,
                                      double *out, const int *onembed,
                                      int ostride, int odist,
                                      unsigned flags);
  *)

  fillchar(massen,sizeof(massen),#0);
  setlength(massen,length(teilchen));
  for i:=0 to length(massen)-1 do
    massen[i]:=0;
  for i:=0 to aX-1 do
    for j:=0 to length(teilchen)-1 do begin
      dens:=teilchen[j].gibDichte(gitter.xl+(i-2)*gitter.dX,parent.kvs);
      initialisiereDichte(i,j,teilchen[j].breite,dens);
      massen[j]:=massen[j]+dens*gitter.dX*teilchen[j].eigenschaften[tsgMasse];
    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,j: longint;
  abl: boolean;
  emF: tEMFeldGroesze;
  emQ: tEMQuellGroesze;
  maF: tMaterieFeldGroesze;
begin
  setlength(massen,0);
  for i:=0 to length(teilchen)-1 do
    teilchen[i].free;
  setlength(teilchen,0);
  for i:=0 to length(impulsRaumGradient)-1 do
    for j:=0 to length(impulsRaumGradient[i])-1 do
      fftw_freemem(impulsRaumGradient[i,j]);
  setlength(iMGamma,0);
  for i:=0 to length(iMGamma)-1 do
    fftw_freemem(iMGamma[i]);
  setlength(iMGamma,0);
  for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do
    for abl:=false to not(emF in [efEX,efBZ]) do
      fftw_freemem(emFelder[emF,abl]);
  for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do
    fftw_freemem(emQuellen[emQ]);
  for i:=0 to length(matFelder)-1 do
    for maF:=low(tMaterieFeldGroesze) to high(tMaterieFeldGroesze) do
      fftw_freemem(matFelder[i,maF]);
  setlength(matFelder,0);
  for i:=0 to length(impulsraum)-1 do
    for abl:=false to true do
      fftw_freemem(impulsraum[i,abl]);
  setlength(impulsraum,0);
  lichters[false].free;
  lichters[true].free;
  for i:=0 to length(ffts)-1 do
    for j:=0 to 1 do
      for abl:=false to true do
        fftw_destroy_plan(ffts[i,j,abl]);
  setlength(ffts,0);
  fftw_freemem(fftTmp);
  inherited destroy;
end;

procedure tFelder.initialisiereDichte(ort,tlc: longint; breite,n: double);
var
  i:   longint;
  ges: double;
begin
  ges:=0;
  for i:=0 to aP-1 do begin
    (impulsraum[tlc,false]+i+ort*aP)^:=power(1/2,sqr((i-aP/2)*dP*2/breite));
    ges:=ges + (impulsraum[tlc,false]+i+ort*aP)^;
  end;
  ges:=n/ges/dP;
  for i:=0 to aP-1 do
    (impulsraum[tlc,false]+i+ort*aP)^:=(impulsraum[tlc,false]+i+ort*aP)^ * ges;
end;

procedure tFelder.pruefeArrayEnden(fehler: string);
var
  i,j: longint;
  abl: boolean;
begin
  if (abs(lFftTmp^.re-pi)>1e-10) or (abs(lFftTmp^.im-sqr(pi))>1e-10) then begin
    gitter.prot.schreibe(fehler);
    gitter.prot.schreibe(floattostr(lFftTmp^.re)+' vs. '+floattostr(pi)+' und '+floattostr(lFftTmp^.im)+' vs. '+floattostr(sqr(pi)));
    raise exception.create('Fehler in tFelder.pruefeArrayEnden!');
  end;
  for i:=0 to length(impulsRaumGradient)-1 do
    for j:=0 to 1 do
      if abs((impulsRaumGradient[i,j]+aP*aX)^-pi)>1e-10 then begin
        gitter.prot.schreibe(fehler);
        gitter.prot.schreibe(floattostr((impulsRaumGradient[i,j]+aP*aX)^)+' vs. '+floattostr(pi));
        raise exception.create('Fehler in tFelder.pruefeArrayEnden!');
      end;
  for i:=0 to length(impulsRaum)-1 do
    for abl:=false to true do
      if abs((impulsRaum[i,abl]+aP*aX)^-pi)>1e-10 then begin
        gitter.prot.schreibe(fehler);
        gitter.prot.schreibe(floattostr((impulsRaum[i,abl]+aP*aX)^)+' vs. '+floattostr(pi));
        raise exception.create('Fehler in tFelder.pruefeArrayEnden!');
      end;
end;

procedure tFelder.berechneAbleitungen;
var
  i,j,k:  longint;
  emQ:    tEMQuellGroesze;
  emF:    tEMFeldGroesze;
  rechts: boolean;
  tmp:    double;
begin

  // PX ggf. korrigieren
  for i:=0 to length(teilchen)-1 do      // linker Rand
    for j:=0 to (aP div 2)-1 do          // negativer Impuls
      if (impulsraum[i,false]+j)^>0 then begin
        (impulsraum[i,false]+aP-j-byte(j=0))^:=(impulsraum[i,false]+aP-j-byte(j=0))^ + (impulsraum[i,false]+j)^;
        (impulsraum[i,false]+j)^:=0;
      end;

  for i:=0 to length(teilchen)-1 do      // rechter Rand
    for j:=(aP div 2)+1 to aP-1 do       // positiver Impuls
      if (impulsraum[i,false]+j+(aX-1)*aP)^>0 then begin
        (impulsraum[i,false]+aP-j+(aX-1)*aP)^:=(impulsraum[i,false]+aP-j+(aX-1)*aP)^ + (impulsraum[i,false]+j+(aX-1)*aP)^;
        (impulsraum[i,false]+j+(aX-1)*aP)^:=0;
      end;

  // PY und iMGamma berechnen
  for i:=0 to length(matFelder)-1 do
    for j:=0 to aX-1 do begin
      (matFelder[i,mfPY]+j)^:=-teilchen[i].eigenschaften[tsgLadung] * (emFelder[efAY,false]+j)^;
      for k:=0 to aP-1 do
        (iMGamma[i]+k+j*aP)^:=
          1/(sqrt(1 + (sqr((k-(aP div 2))*dP) + sqr((matFelder[i,mfPY]+j)^))*teilchen[i].eigenschaften[tsgIQdrMasse])*teilchen[i].eigenschaften[tsgMasse]);
        // 1/(m Gamma) = 1/sqrt(1 + (pX/m)^2 + (pY/m)^2)*m
    end;

  // rho und j berechnen
  for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do
    for i:=0 to aX-1 do
      (emQuellen[emQ]+i)^:=0;
  for i:=0 to length(matFelder)-1 do
    for j:=0 to aX-1 do
      for k:=0 to aP-1 do begin
        tmp:=(impulsraum[i,false]+k+j*aP)^ * teilchen[i].eigenschaften[tsgLadung] * dP;
        (emQuellen[eqRho]+j)^:=(emQuellen[eqRho]+j)^ + tmp;    // * iGamma[i] ??? TODO ... falls, dann aber ohne M!
        (emQuellen[eqJX]+j)^:= (emQuellen[eqJX]+j)^  + tmp * (k-(aP div 2))*dP * (iMGamma[i]+k+j*aP)^;
        (emQuellen[eqJY]+j)^:= (emQuellen[eqJY]+j)^  + tmp * (matFelder[i,mfPY]+j)^ * (iMGamma[i]+k+j*aP)^;
      end;

  // E und B aus Rho und A berechnen

  // dEX/dx = rho  ... hier muss integriert werden:
  // EX = EX(x- deltax) + <rho> * deltax

  emFelder[efEX,false]^:=emQuellen[eqRho]^ * dX; // linker Rand
  for i:=1 to aX-1 do                            // Rest
    (emFelder[efEX,false]+i)^:=
      (emFelder[efEX,false]+i-1)^ +
      ((emQuellen[eqRho]+i-1)^ + (emQuellen[eqRho]+i)^) * dX/2;

  // B = rot A
  emFelder[efBZ,false]^:=0;        // linker Rand
  for i:=1 to aX-2 do              // Mitte
    (emFelder[efBZ,false]+i)^:=
      ((emFelder[efAY,false]+i+1)^ - (emFelder[efAY,false]+i-1)^)*iDX/2;
  (emFelder[efBZ,false]+aX-1)^:=0; // rechter Rand

  // Randbedingungen für die Ableitungen des A-Felds setzen
  for emF:=efAX to efAY do begin // Vakuumrandbedingungen für das A-Feld
    emFelder[emF,true]^:=
      ((emFelder[emF,false]+1)^ -
       emFelder[emF,false]^)*iDX;
    (emFelder[emF,true]+aX-1)^:=
      ((emFelder[emF,false]+aX-2)^ -
       (emFelder[emF,false]+aX-1)^)*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
      (emFelder[efAY,true]+(aX-1)*byte(rechts))^:=
        (emFelder[efAY,true]+(aX-1)*byte(rechts))^ +
        exprToFloat(false,lichters[rechts][i],gitter.kvs,nil);

  // sonstige Ableitungen des A-Feldes berechnen

  // d2A/dt2 = Laplace(A) - j    ( - dE/dt wird auf dA/dt direkt aufgeschlagen !!! )
  for i:=1 to aX-2 do
    (emFelder[efDAXDT,true]+i)^:=
      ( (emFelder[efAX,false]+i+1)^ - 2*(emFelder[efAX,false]+i)^ + (emFelder[efAX,false]+i-1)^ )*sqr(iDX)
      - (emQuellen[eqJX]+i)^;
  for i:=1 to aX-2 do
    (emFelder[efDAYDT,true]+i)^:=
      ( (emFelder[efAY,false]+i+1)^ - 2*(emFelder[efAY,false]+i)^ + (emFelder[efAY,false]+i-1)^ )*sqr(iDX)
      - (emQuellen[eqJY]+i)^;

  // dA/dt = dA/dt - E
  for i:=1 to aX-2 do
    (emFelder[efAX,true]+i)^:=(emFelder[efDAXDT,false]+i)^ - (emFelder[efEX,false]+i)^;
  move(emFelder[efDAYDT,false]^,emFelder[efAY,true]^,aX*sizeof(double));

  // Gradienten der Phasenraumbesetzungsdichte berechnen

  pruefeArrayEnden('Jemand hat das FFT-Array überschrieben!');

  for i:=0 to length(ffts)-1 do begin
    fftw_execute(ffts[i,0,false]);     // FFT in x
    pruefeArrayEnden('Vorwärts-FFT ('+inttostr(i)+',0) schreibt beim Ausführen über das Ende hinaus!');

    for j:=0 to aX div 2 do
      for k:=0 to aP-1 do begin  //    *-i*k*2*pi;
        tmp:=(fftTmp+j+k*(aX div 2 + 1))^.re;
        (fftTmp+j+k*(aX div 2 + 1))^.re:=  2*pi/dX/aX * j * (fftTmp+j+k*(aX div 2 + 1))^.im/aX * byte(3*j<aX);
        (fftTmp+j+k*(aX div 2 + 1))^.im:=- 2*pi/dX/aX * j * tmp/aX                             * byte(3*j<aX);
      end;

    fftw_execute(ffts[i,0,true]);
    pruefeArrayEnden('Rückwärts-FFT ('+inttostr(i)+',0) schreibt beim Ausführen über das Ende hinaus!');

    fftw_execute(ffts[i,1,false]);     // FFT in p
    pruefeArrayEnden('Vorwärts-FFT ('+inttostr(i)+',1) schreibt beim Ausführen über das Ende hinaus!');

    for j:=0 to aX-1 do
      for k:=0 to aP div 2 do begin  //    *-i*k*2*pi;
        tmp:=(fftTmp+j+k*aX)^.re;
        (fftTmp+j+k*aX)^.re:=  2*pi/dP/aP * k * (fftTmp+j+k*aX)^.im/aP * byte(4*k < aP);
        (fftTmp+j+k*aX)^.im:=- 2*pi/dP/aP * k * tmp/aP                 * byte(4*k < aP);
      end;

    fftw_execute(ffts[i,1,true]);
    pruefeArrayEnden('Rückwärts-FFT ('+inttostr(i)+',1) schreibt beim Ausführen über das Ende hinaus!');

  end;

  // Zeitableitung der Phasenraumbesetzungsdichte berechnen

  for i:=0 to length(teilchen)-1 do
    for j:=1 to aX-2 do
      for k:=1 to aP-2 do
        (impulsraum[i,true]+k+j*aP)^:=                                               // df/dt =
          - (k-(aP div 2))*dP * (iMGamma[i]+k+j*aP)^                                 //   - v
            * (impulsRaumGradient[i,0]+k+j*aP)^                                      //   * Nabla f
          - teilchen[i].eigenschaften[tsgLadung] *                                   //   - q(
            (  (emFelder[efEX,false]+j)^                                             //        E
             + (matFelder[i,mfPY]+j)^ * (iMGamma[i]+k+j*aP)^ * (emFelder[efBZ,false]+j)^) // + v/c x B)
            * (impulsRaumGradient[i,1]+k+j*aP)^;                                     //   * Nabla_p f
end;

procedure tFelder.setzeNull;
var
  i,j,k: longint;
  abl:   boolean;
  emF:   tEMFeldGroesze;
  emQ:   tEMQuellGroesze;
  maF:   tMaterieFeldGroesze;
begin
  for i:=0 to length(impulsRaumGradient)-1 do
    for j:=0 to length(impulsRaumGradient[i])-1 do
      for k:=0 to aX*aP-1 do
        (impulsRaumGradient[i,j]+k)^:=0;
  for i:=0 to length(iMGamma)-1 do
    for j:=0 to aX*aP-1 do
      (iMGamma[i]+j)^:=1;
  for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do
    for abl:=false to not(emF in [efEX,efBZ]) do
      for i:=0 to aX-1 do
        (emFelder[emF,abl]+i)^:=0;
  for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do
    for i:=0 to aX-1 do
      (emQuellen[emQ]+i)^:=0;
  for i:=0 to length(matFelder)-1 do
    for maF:=low(tMaterieFeldGroesze) to high(tMaterieFeldGroesze) do
      for j:=0 to aX-1 do
        (matFelder[i,maF]+j)^:=0;
  for i:=0 to length(impulsraum)-1 do
    for abl:=false to true do
      for j:=0 to aX*aP-1 do
        (impulsraum[i,abl]+j)^:=0;
end;

procedure tFelder.dumpErhaltungsgroeszen;
var
  i,j:  longint;
  dens: double;
begin
  dens:=0;
  for i:=0 to length(massen)-1 do
    dens:=dens+massen[i]/teilchen[i].eigenschaften[tsgMasse];
  gitter.prot.schreibe('Gesamtdefizit: '+floattostr(gesamtDefizit)+' (Anteil '+floattostr(gesamtDefizit/dens)+')',true);
  for i:=0 to length(massen)-1 do begin
    dens:=0;
    for j:=0 to aX-1 do
      dens:=dens+impulsIntegral(j,i,msN);
    dens:=dens*teilchen[i].eigenschaften[tsgMasse]*gitter.dX;
    gitter.prot.schreibe('n['+inttostr(i+1)+'] = '+floattostr(dens)+' (relative Abweichung: '+floattostr(dens/massen[i]-1)+')',true);
  end;
end;

function tFelder.impulsIntegral(ort,tlc: longint; emQ: tEMQuellGroesze): double;
begin
  if tlc<0 then begin
    result:=(emQuellen[emQ]+ort)^; // das ist leicht :-)
    exit;
  end;

  result:=0;

  case emQ of
    eqRho: result:=impulsIntegral(ort,tlc,msN);
    eqJX:  result:=impulsIntegral(ort,tlc,msVX);
    eqJY:  result:=impulsIntegral(ort,tlc,msVY);
    else begin
      gitter.prot.schreibe('unbekannte EMQuellgröße in Impulsintegral!');
      raise exception.create('Fehler in tFelder.impulsIntegral!');
    end;
  end{of case};
  result:=result * teilchen[tlc].eigenschaften[tsgLadung];
end;

function tFelder.impulsIntegral(ort,tlc: longint; maF: tMaterieSpeicherGroesze): double;
var
  i: longint;
begin
  if tlc<0 then begin
    result:=0;
    for i:=0 to length(teilchen)-1 do
      result:=result+impulsIntegral(ort,i,maF);
    exit;
  end;

  result:=0;

  case maF of
    msN:
      for i:=0 to aP-1 do
        result:=
          result +
          (impulsraum[tlc,false]+i+ort*aP)^ * dP;
    msPX:
      for i:=0 to aP-1 do
        result:=
          result +
          (impulsraum[tlc,false]+i+ort*aP)^ * dP * (i-aP/2)*dP;
    msPXSqr:
      for i:=0 to aP-1 do
        result:=
          result +
          (impulsraum[tlc,false]+i+ort*aP)^ * dP * sqr((i-aP/2)*dP);
    msPY: begin
      for i:=0 to aP-1 do
        result:=
          result +
          (impulsraum[tlc,false]+i+ort*aP)^ * dP;
      result:=result * (matFelder[tlc,mfPY]+ort)^;
    end;
    msVX: begin
      for i:=0 to aP-1 do
        result:=
          result +
          (impulsraum[tlc,false]+i+ort*aP)^ * dP * (i-aP/2)*dP
          / sqrt(1 + (sqr((i-aP/2)*dP) + sqr((matFelder[tlc,mfPY]+ort)^))*teilchen[tlc].eigenschaften[tsgIQdrMasse]);
      result:=result / teilchen[tlc].eigenschaften[tsgMasse];
    end;
    msVY: begin
      for i:=0 to aP-1 do
        result:=
          result +
          (impulsraum[tlc,false]+i+ort*aP)^ * dP
          / sqrt(1 + (sqr((i-aP/2)*dP) + sqr((matFelder[tlc,mfPY]+ort)^))*teilchen[tlc].eigenschaften[tsgIQdrMasse]);
      result:=result * (matFelder[tlc,mfPY]+ort)^ / teilchen[tlc].eigenschaften[tsgMasse];
    end;
    msPXRipple:
      for i:=1 to aP-1 do
        result:=
          result + abs((impulsraum[tlc,false]+i+ort*aP)^-(impulsraum[tlc,false]+i-1+ort*aP)^)*dP;
    else begin
      gitter.prot.schreibe('unbekannte Materiefeldgröße in Impulsintegral!');
      raise exception.create('Fehler in tFelder.impulsIntegral!');
    end;
  end{of case};
end;

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

procedure tFelder.schreibePhasenraum(tlc: longint; datNam: string);
var
  f:   textfile;
  i,j: longint;
begin
  assignfile(f,datNam);
  rewrite(f);
  for i:=0 to aP-1 do begin
    for j:=0 to aX-1 do begin
      if j>0 then write(f,#9);
      write(f,(impulsraum[tlc,false]+i+j*aP)^);
    end;
    writeln(f);
  end;
  closefile(f);
end;

procedure tFelder.nichtnegativieren; // Dichten nicht negativ machen
var
  i,j,k,l: longint;
  defizit: double;
begin
  for i:=0 to length(impulsraum)-1 do
    for j:=0 to aX*aP-1 do
      if (impulsraum[i,false]+j)^ < 0 then begin
        defizit:=-(impulsraum[i,false]+j)^;
        gesamtDefizit:=gesamtDefizit + defizit;
        (impulsraum[i,false]+j)^:=0;
        k:=j;
        l:=1;
        while (defizit>0) and (l<2*aX*aP) do begin
          k:=k+l*(1-2*byte(odd(l)));
          inc(l);
          if (k<0) or (k>=aX*aP) then continue;
          if (impulsraum[i,false]+k)^>defizit then begin
            (impulsraum[i,false]+k)^:=(impulsraum[i,false]+k)^-defizit;
            defizit:=0;
          end
          else if (impulsraum[i,false]+k)^>0 then begin
            defizit:=defizit-(impulsraum[i,false]+k)^;
            (impulsraum[i,false]+k)^:=0;
          end;
        end;
        if defizit>0 then begin
          gitter.prot.schreibe('Kann Defizit der Teilchensorte '+inttostr(i+1)+' nicht ausgleichen, '+floattostr(defizit)+' bleibt übrig!');
          gitter.abbrechen;
        end;
      end;
end;

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

constructor tGitter.create(besitzer: tSimulation; aX,aP: longint; deltaX,deltaP: double; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; name: string);
var
  i: longint;
begin
  inherited create;

  i:=aP;
  aP:=round(power(2,ceil(ln(aP-0.5)/ln(2))));
  deltaP:=deltaP*(i-1)/(aP-1);

  i:=aX;
  aX:=round(power(2,ceil(ln(aX+4-0.5)/ln(2)))); // zwei Felder links und rechts extra für Randbedingungen
  deltaX:=deltaX*(i-1)/(aX-4-1);

  abbruch:=false;
  simulation:=besitzer;
  zeitverfahren:=zv;
  kvs:=bekannteWerte;
  prot:=tProtokollant.create(simulation.prot,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(aX,aP,deltaX,deltaP,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.abbrechen;
var
  i: longint;
begin
  if not abbruch then
    for i:=0 to length(felders[aktuelleFelder].matFelder)-1 do
      felders[aktuelleFelder].schreibePhasenraum(i,extractfilepath(paramstr(0))+'phasenraum_'+inttostr(i+1)+'.dump');
  abbruch:=true;
end;

procedure tGitter.iteriereSchritt(dT: double);
    {$IFDEF negativeDichteueberwachung}
var i,j:     longint;
    pro:     tProtokollant;
    {$ENDIF}
begin
  if simulation.gotSigterm then begin
    abbrechen;
    simulation.gotSigterm:=false;
  end;

  if abbruch then begin
    t:=t+dT;
    exit;
  end;
  kvs.add('t',t);

  repeat
    felders[aktuelleFelder].berechneAbleitungen; // 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].impulsraum)-1 do
    for j:=0 to felders[aktuelleFelder].aX*felders[aktuelleFelder].aP-1 do
      if (felders[aktuelleFelder].impulsraum[i,false]+j)^<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);
        abbrechen;
      end;
  {$ENDIF}
  t:=t+dT;
end;

procedure tGitter.dumpErhaltungsgroeszen;
begin
  felders[aktuelleFelder].dumpErhaltungsgroeszen;
end;

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

constructor tSimulation.create(inName: string; protokollant: tProtokollant; name: string);
var
  ifile:                 tMyStringlist;
  zeitverfahren:         tZeitverfahren;
  s,t,aSuffix,aPrefix:   string;
  dX,breite,dP,pMax:     double;
  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);
    raise exception.create('Fehler in tSimulation.create!');
  end;

  // Standardeinstellungen Bereich 'allgemein'
  zeitverfahren:=zfRungeKuttaVier;
  dX:=1e-2;
  kvs.add('λ',1);
  kvs.add('T',1);
  kvs.add('ω',2*pi);
  kvs.add('dX',dX);
  kvs.add('qe',1);
  kvs.add('me',1);
  dT:=-1;
  sDT:=-1;
  sDX:=-1;
  tEnde:=100;
  breite:=10.0;
  fortschrittsAnzeige:=false;
  gotSighup:=false;
  gotSigterm:=false;
  dP:=-1;
  pMax:=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);
      raise exception.create('Fehler in tSimulation.create!');
    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);
          raise exception.create('Fehler in tSimulation.create!');
        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
          dX:=exprtofloat(false,s,kvs,nil);
          kvs.add('dX',dX);
          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
          dP:=exprtofloat(false,s,kvs,nil);
          continue;
        end;
        if startetMit('maximalimpuls ',s) then begin
          pMax:=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);
        raise exception.create('Fehler in tSimulation.create!');
      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);
          raise exception.create('Fehler in tSimulation.create!');
        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);
        raise exception.create('Fehler in tSimulation.create!');
      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);
      raise exception.create('Fehler in tSimulation.create!');
    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);
        raise exception.create('Fehler in tSimulation.create!');
      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);
          raise exception.create('Fehler in tSimulation.create!');
        end;
        if s='teilchen'+inttostr(length(teilchen))+'Ende' then break;
        if startetMit('ladung ',s) then begin
          teilchen[length(teilchen)-1].eigenschaften[tsgLadung]:=-exprToFloat(false,s,kvs,nil);
          kvs.add('q'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].eigenschaften[tsgLadung]);
          continue;
        end;
        if startetMit('masse ',s) then begin
          teilchen[length(teilchen)-1].eigenschaften[tsgMasse]:=exprToFloat(false,s,kvs,nil);
          teilchen[length(teilchen)-1].eigenschaften[tsgIQdrMasse]:=1/sqr(teilchen[length(teilchen)-1].eigenschaften[tsgMasse]);
          kvs.add('m'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].eigenschaften[tsgMasse]);
          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;
        if startetMit('maximales dLnN/dX',s) then begin
          teilchen[length(teilchen)-1].diffusion[0]:=exprToFloat(false,s,kvs,nil);
          continue;
        end;
        if startetMit('maximales dLnN/dP',s) then begin
          teilchen[length(teilchen)-1].diffusion[1]:=exprToFloat(false,s,kvs,nil);
          continue;
        end;
        if startetMit('impulsbreite ',s) then begin
          teilchen[length(teilchen)-1].breite:=exprToFloat(false,s,kvs,nil);
          continue;
        end;

        prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich teilchen!',true);
        raise exception.create('Fehler in tSimulation.create!');
      until false;
      continue;
    end;

    prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+'''!',true);
    raise exception.create('Fehler in tSimulation.create!');
  until false;

  ifile.free;

  if length(ausgabedateien)=0 then begin
    prot.schreibe('Du solltest irgendetwas abspeichern lassen!',true);
    raise exception.create('Fehler in tSimulation.create!');
  end;

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

  if dP<0 then
    dP:=dX;

  for i:=0 to length(teilchen)-1 do begin
    teilchen[i].diffusion[0]:=teilchen[i].diffusion[0] * dX;
    teilchen[i].diffusion[1]:=teilchen[i].diffusion[1] * dP;
  end;


  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(dX));
  pro.schreibe('dT   =   '+floattostr(dT));
  pro.schreibe('breite = '+floattostr(breite)+' (= '+inttostr(round(breite/dX)+1)+' Punkte)');
  pro.schreibe('dP   =   '+floattostr(dP));
  pro.schreibe('pMax  =  '+floattostr(pMax)+' (= '+inttostr(2*round(pMax/dP)+1)+' Punkte)');
  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(simulationen)=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(SIGUSR1, na, nil) <> 0) or
       (fPSigaction(SIGTERM, na, nil) <> 0) then begin
      pro.schreibe('Fehler: '+inttostr(fpgeterrno)+'.');
      raise exception.create('Fehler in tSimulation.create!');
    end;
    dispose(na);
  end;
  setlength(simulationen,length(simulationen)+1);
  simulationen[length(simulationen)-1]:=self;

  gitter:=tGitter.create(self,round(breite/dX)+1,2*round(pMax/dP)+1,dX,dP,kvs,teilchen,lichter,zeitverfahren,'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(simulationen)-1 downto 0 do
    if simulationen[i]=self then begin
      for j:=i+1 to length(simulationen)-1 do
        simulationen[j-1]:=simulationen[j];
      setlength(simulationen,length(simulationen)-1);
    end;
  inherited destroy;
end;

function tSimulation.iteriereSchritt(start: double; var zeitPhysik,zeitDatei: double): 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 try
    gitter.iteriereSchritt(dT);
  except
    zeitPhysik:=zeitPhysik+now;
    exit;
  end;
  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);
    gitter.dumpErhaltungsgroeszen;
  end;

  result:=gitter.t<tEnde;
end;

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

procedure SignalCapture(signal : longint); cdecl;
var
  i: longint;
begin
  case signal of
    SIGUSR1: begin
      writeln('received SIGUSR1!');
      for i:=0 to length(simulationen)-1 do
        simulationen[i].gotSighup:=true;
    end;
    SIGTERM: begin
      writeln('received SIGTERM!');
      for i:=0 to length(simulationen)-1 do
        simulationen[i].gotSigterm:=true;
    end;
  end{of case};
end;

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