unit Physikunit;

{$mode objfpc}{$H+}

{$DEFINE Zeitschrittueberwachung}
{$DEFINE Dichteueberwachung}

interface

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

type
  tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaVier);
  tVerteilungsfunktion = function(x: extended): extended;
  tEMFeldInhalt = (
    efA,efAX,efAY,efAZ,
    efDAXDT,efDAYDT,efDAZDT,
    efDPhiDX
  );
  tMaterieFeldInhalt = (
    mfN,mfDPsiDX,
    mfP,mfPX,mfPY,mfPZ,
    mfGamma,mfIGamma
  );

  { tAusgabeDatei }

  tAusgabeDatei = class
    ableitung: boolean;
    emInhalt:  tEMFeldInhalt;
    matInhalt: tMaterieFeldInhalt;
    teilchen:  longint; // wenn < 0, dann EM-Feld
    name:      string;
    datei:     file;
    constructor create(feldName,prefix,suffix: string; prot: tProtokollant);
    destructor destroy; override;
  end;

  tGitter = class;

  { tWertePunkt }

  tWertePunkt = class(tObject) // repräsentiert alle Werte eines Punktes im Gitter und deren Zeitableitungen
    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, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier!
    lN,rN:    tWertePunkt;

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

  { TFelder }

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

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

  { tGitter }

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

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

    constructor create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant);
    destructor destroy; override;
    procedure iteriereSchritt(dT: extended);
    procedure macheAusgabe(aD: tAusgabedatei; sDX: extended);
    function gibErhaltungsgroessen: string;
  end;

  { tSimulation }

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

implementation

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

{ tAusgabeDatei }

constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokollant);
var
  emF:     tEMFeldInhalt;
  maF:     tMaterieFeldInhalt;
  num:     longint;
  abl,gef: boolean;
begin
  inherited create;
  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)));
    delete(feldName,num,length(feldName));
  end
  else
    teilchen:=0;

  emInhalt:=efA;
  matInhalt:=mfN;
  feldName:=ansiUpperCase(feldName);

  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
    prot.schreibe('tAusgabeDatei.create kennt Feld '''+feldName+''' nicht!',true);
    halt(1);
  end;

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

  if abl then
    name:='d'+name;
  name:=prefix+name+suffix;

  assignFile(datei,name);
  rewrite(datei);
end;

destructor tAusgabeDatei.destroy;
begin
  closeFile(datei);
  inherited destroy;
end;

{ tWertePunkt }

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

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

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

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

procedure tWertePunkt.berechneNAbleitung(iDX,pDNMax: extended); // Zeitableitung von n berechnen
var
  i: longint;
begin
(*  // dn/dt = -d(n/Gamma p_x)/dx
    Groessen[true,fiN]:=
      -(   rN^.rN^.Groessen[false,fiN]*rN^.rN^.Groessen[false,fiiGamma]*rN^.rN^.Groessen[false,fiPX]
         + rN^.Groessen[false,fiN]*rN^.Groessen[false,fiiGamma]*rN^.Groessen[false,fiPX]
         - lN^.Groessen[false,fiN]*lN^.Groessen[false,fiiGamma]*lN^.Groessen[false,fiPX]
         - lN^.lN^.Groessen[false,fiN]*lN^.lN^.Groessen[false,fiiGamma]*lN^.lN^.Groessen[false,fiPX]
       )*iDX/6;
    // es wird über die beiden nächsten linken bzw. rechten Nachbarn gemittelt   *)

  for i:=0 to length(matWerte)-1 do
    // dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n
    matWerte[i,mfN,true]:=
      ( rN.matWerte[i,mfN,false] * rN.matWerte[i,mfIGamma,false] *  (pDNMax - rN.matWerte[i,mfPX,false])
      - ln.matWerte[i,mfN,false] * ln.matWerte[i,mfIGamma,false] * (-pDNMax - rN.matWerte[i,mfPX,false])
      - matWerte[i,mfN,false]  *  matWerte[i,mfIGamma,false] * 2 * pDNMax) * iDX/2;
    // Der zweite Summand entspricht (in etwa) einer thermischen Verschmierung
    // und soll die numerische Stabilität bis zu p * dn/dx / n von 2*pDNMax gewährleisten.
    // Es handelt sich um einen Diffusionsterm dn/dt = alpha * Laplace n mit
    // alpha = pDNMax * dX
end;

procedure tWertePunkt.berechneAbleitungen(dT,dX,iDX: extended); // Zeitableitungen berechnen
var
  i: longint;
begin
  // d2Psi/dxdt = dPhi/dx - dGamma/dx
  for i:=0 to length(matWerte)-1 do
    matWerte[i,mfDPsiDX,true]:=
      emWerte[efDPhiDX,false]
        - (rN.matWerte[i,mfGamma,false] - lN.matWerte[i,mfGamma,false]) * iDX/2;

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

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

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

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

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

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

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

{ tFelder }

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

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

procedure tFelder.setzeRaender(iDX: extended);
var
  emF: tEMFeldInhalt;
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)

  (*  inhalt[0].emWerte[efAY,true]:=
    inhalt[0].emWerte[fiAY,true]
     + (Al(t+dTMaximum)-Al(t))/dTMaximum;
  Felders[Felder,groesse+3].Groessen[true,fiAy]:=
    Felders[Felder,groesse+3].Groessen[true,fiAy]
     + (Ar(t+dTMaximum)-Ar(t))/dTMaximum; *)
end;

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

  setzeRaender(iDX);

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

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

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

{ tGitter }

constructor tGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant);
var
  i: longint;
begin
  inherited create;
  zeitverfahren:=ZV;
  kvs:=bekannteWerte;
  Prot:=TProtokollant.create(Protokollant,'tGitter');
  dTMaximum:=deltaT;
  dX:=deltaX;
  iDX:=1/dX;
  if pDNMa < 0 then begin
    pDNMax:=0;
    pDNMaxDynamisch:=true;
  end
  else pDNMax:=pDNMa;
  case Zeitverfahren of
    zfEulerVorwaerts: Setlength(Felders,2);
    zfRungeKuttaVier: Setlength(Felders,5);
  end{of Case};
  xl:=dX/2;

  for i:=0 to length(felders)-1 do
    felders[i]:=tFelder.create(size,dichten,lichter,self);

  aktuelleFelder:=0;
  t:=0;
end;

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

procedure tGitter.iteriereSchritt(dT: extended);
var i:     longint;
    {$IFDEF Zeitschrittueberwachung}
    pro:   tProtokollant;
    j:     longint;
    {$ENDIF}
begin
  felders[aktuelleFelder].berechneAbleitungen(dT,dX,iDX,pDNMax); // 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);
    zfRungeKuttaVier:
      begin
        felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT/2); // ya = y(t) + y' dt/2
        felders[2].berechneAbleitungen(dT/2,dX,iDX,pDNMax);              // ya' = y'(t+dt/2,ya)

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

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

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

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

procedure tGitter.macheAusgabe(aD: tAusgabedatei; sDX: extended);
var i:         longint;
    sX,cX,tmp: extended;
    emF:       tEMFeldInhalt;
    maF:       tMaterieFeldInhalt;
begin
  if (aD.teilchen<0) and (aD.emInhalt=efA) then
    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;
  if (aD.teilchen>=0) and (aD.matInhalt=mfP) then
    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[aD.teilchen,maF,false]);
      felders[aktuelleFelder].inhalt[i].matWerte[aD.teilchen,maF,false]:=sqrt(tmp);
    end;

  i:=floor((length(felders[aktuelleFelder].inhalt)-1)/sDX*dX+1);
  sX:=xl+(i-1)*sDX;
  blockWrite(aD.datei,xl,sizeof(extended));
  blockWrite(aD.datei,sX,sizeof(extended));
  blockWrite(aD.datei,i,sizeof(longint));
  sX:=xl-Min(dX,sDX)/2;
  cX:=xl;
  if aD.teilchen<0 then begin
    for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do begin
      while cX>=sX do begin
        blockWrite(aD.datei,felders[aktuelleFelder].inhalt[i].emWerte[aD.emInhalt,aD.ableitung],sizeof(extended));
        sX:=sX+sDX;
      end;
      cX:=cX+dX;
    end;
  end
  else begin
    for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do begin
      while cX>=sX do begin
        blockWrite(aD.datei,felders[aktuelleFelder].inhalt[i].matWerte[aD.teilchen,aD.matInhalt,aD.ableitung],sizeof(extended));
        sX:=sX+sDX;
      end;
      cX:=cX+dX;
    end;
  end;
end;

function tGitter.gibErhaltungsgroessen: string;
var i,j: integer;
    ns:  tExtendedArray;
    Pro: TProtokollant;
begin
  pro:=tProtokollant.create(prot,'gibErhaltungsgroessen');
  setlength(ns,felders[aktuelleFelder].matAnz);
  result:='';
  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
    result:=result + ' n['+inttostr(i)+']='+floattostr(ns[i]);
    {$IFDEF Dichteueberwachung}
    if ns[i]>1000 then begin
      errorCode:=2;
      pro.schreibe(' n > 1000, es scheinen sehr viele neue Teilchen entstanden zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')');
    end;
    {$ENDIF}
  end;
  delete(result,1,1);
  pro.free;
end;

procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant);
var
  i,j: longint;
begin
  for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do
    for j:=0 to felders[aktuelleFelder].matAnz-1 do
      if abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPx,false]*
              (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
              max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10))
          > pDNMax then begin
        if pDNMaxDynamisch then begin
          pDNMax:=
            2*
            abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPX,false]*
                (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
                Max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10));
          pro.schreibe('Neues maximales px * dn/dx / n: '+floattostr(pDNMax)+' (t='+floattostr(t)+')');
        end
        else
        if pDNMax>0 then begin
          pro.schreibe('Warnung: maximaler Impuls * dn/dx / n von '+floattostr(pDNMax)+' wurde überschritten (t='+floattostr(t)+
            'T), die numerische Stabilität ist nicht mehr gewährleistet!',true);
          pro.schreibe(' Lösung: größeren Diffusionsterm wählen (-D)',true);
          pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern (-x)',true);
          pDNMax:=-1;
        end;
      end;
end;

{ tSimulation }

constructor tSimulation.create(inName: string; Protokollant: tProtokollant);
var
  ifile:                                   tMyStringlist;
  zeitverfahren:                           tZeitverfahren;
  s,t,aSuffix,aPrefix:                     string;
  deltaX,deltaT,breite,pDNMax:             extended;
  i,j:                                     longint;
  kvs:                                     tKnownValues;
  dichten,lichter:                         tMyStringlist;
begin
  inherited create;
  prot:=tProtokollant.create(Protokollant,'tSimulation');
  kvs:=tKnownValues.create;

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

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

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

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

  // Standardeinstellungen Breich 'teilchen...'
  dichten:=tMyStringlist.create(prot);

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

    if s='ausgaben' then begin
      repeat
        if not ifile.readln(s) then begin
          prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich ausgaben!',true);
          halt(1);
        end;
        if s='ausgabenEnde' then break;
        if startetMit('suffix',s) then begin
          aSuffix:=s;
          continue;
        end;
        if startetMit('prefix',s) then begin
          aPrefix:=s;
          continue;
        end;
        if startetMit('zeitschritt',s) then begin
          sDT:=exprtofloat(false,s,kvs,nil);
          continue;
        end;
        if startetMit('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);
          end;
          continue;
        end;
        prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich ausgaben!',true);
        halt(1);
      until false;
      continue;
    end;

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

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

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

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

  gitter:=tGitter.create(round(breite/deltaX),deltaT,deltaX,pDNMax,kvs,dichten,lichter,zeitverfahren,prot);
  dichten.free;
  lichter.free;
end;

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

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

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

  gitter.iteriereSchritt(dT);

  if fortschrittsAnzeige then begin
    if (floor(100*Gitter.t/tEnde) < floor(100*(Gitter.t+dT)/tEnde)) or keypressed then begin
      if keypressed then c:=readkey
      else c:=#0;
      case c of
        #27,'q': begin
          errorCode:=3;
          exit;
        end;
        ' ': begin
          writeln(' ... Pause (beliebige Taste drücken um fortzufahren) ...');
          readkey;
          writeln(' ... weiter geht''s ...');
        end;
        else begin
          Prot.schreibe(inttostr(round(100*Gitter.t/tEnde))+'% (t='+floattostr(Gitter.t)+'T)',true);
          Prot.schreibe(timetostr(now-start)+' ('+floattostr(zeitPhysik/max(1e-11,zeitPhysik+zeitDatei))+')',true);
          Prot.schreibe('ETA: '+timetostr((now-start)*(tEnde-Gitter.t)/max(Gitter.t,dT)),true);
          Prot.schreibe('aktueller Zeitschritt: '+floattostr(dT)+'T',true);
          if errorCode<2 then
            Prot.schreibe(gitter.gibErhaltungsgroessen);
        end;
      end{of case};
    end;
  end;

  result:=gitter.t>=tEnde;
end;

end.