Tagesendstand, alte (= inhaltliche) Bugs sind noch zu beseitigen

5 files changed, 362 insertions, 129 deletions

diff --git a/Physikunit.pas b/Physikunit.pas
index 1d50c5a..8c216c9 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -47,6 +47,22 @@ type
     procedure speichereWerte(gitter: tGitter; sDX: extended);
   end;
 
+  { tTeilchenSpezies }
+
+  tTeilchenSpezies = class
+  private
+    ortsGrenzen:   array of extended;
+    dichteStuecke: array of string;
+  public
+    nMax,                        // maximale Massendichte in nc
+    spezifischeLadung: extended; // q/m in e/me
+    constructor create;
+    destructor destroy; override;
+    procedure dump(prot: tProtokollant; prefix: string);
+    function gibDichte(x: extended; kvs: tKnownValues): extended; // Massendichte
+    procedure liesDichteFunktion(rest: string; ifile: tMyStringList; kvs: tKnownValues; teilchen: array of tTeilchenSpezies; prot: tProtokollant);
+  end;
+
   { tWertePunkt }
 
   tWertePunkt = class(tObject) // repräsentiert alle Werte eines Punktes im Gitter und deren Zeitableitungen
@@ -64,6 +80,7 @@ type
     procedure berechneAbleitungen(dX,iDX: extended);
     procedure liKo(in1,in2: tWertePunkt; fak2: extended); overload;                            // Werte werden auf (in1 + fak2*in2') gesetzt
     procedure liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+    function maxDT: extended;
   end;
 
   { TFelder }
@@ -77,11 +94,12 @@ type
     inhalt: array of tWertePunkt;
     par:    tGitter;
 
-    constructor create(groesse: longint; dichten,lichter: tMyStringList; parent: tGitter);
+    constructor create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
     destructor destroy; override;
     procedure berechneAbleitungen(dT,dX,iDT,iDX,pDNMax: extended);
     procedure liKo(in1,in2: tFelder; fak2: extended); overload;                            // Werte werden auf (in1 + fak2*in2') gesetzt
     procedure liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+    function maxDT: extended;
   end;
 
   { tGitter }
@@ -99,10 +117,10 @@ type
     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; name: string);
+    constructor create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string);
     destructor destroy; override;
-    procedure iteriereSchritt(dT,iDT: extended);
-    function gibErhaltungsgroessen: string;
+    procedure iteriereSchritt(var dT: extended);
+    function dumpErhaltungsgroessen: boolean;
     procedure berechne(was: char; teilchen: longint);
   end;
 
@@ -147,7 +165,7 @@ begin
     dec(num);
   inc(num);
   if num<=length(feldName) then begin
-    teilchen:=strtoint(copy(feldName,num,length(feldName)));
+    teilchen:=strtoint(copy(feldName,num,length(feldName)))-1;
     delete(feldName,num,length(feldName));
   end
   else
@@ -186,7 +204,7 @@ begin
   end;
 
   if teilchen>=0 then
-    pre:=matFeldNamen[maF]+inttostr(teilchen)
+    pre:=matFeldNamen[maF]+inttostr(teilchen+1)
   else
     pre:=emFeldNamen[emF];
 
@@ -201,6 +219,7 @@ destructor tAusgabeDatei.destroy;
 begin
   if (tCnt<>zeitAnz) and ((nNum<>0) or (tCnt<>-1)) then begin
     pro.schreibe('Falsche Anzahl an Zeitschritten in '+pre+'-'+inttostr(nNum-1)+suf+' geschrieben ('+inttostr(tCnt)+' statt '+inttostr(zeitAnz)+')!',true);
+    pro.destroyall;
     halt(1);
   end;
   inherited destroy;
@@ -221,6 +240,7 @@ procedure tAusgabeDatei.schreibeKopf;
 begin
   if (tCnt<>zeitAnz) and ((nNum<>0) or (tCnt<>-1)) then begin
     pro.schreibe('Falsche Anzahl an Zeitschritten in '+pre+'-'+inttostr(nNum-1)+suf+' geschrieben ('+inttostr(tCnt)+' statt '+inttostr(zeitAnz)+')!',true);
+    pro.destroyall;
     halt(1);
   end;
   tCnt:=0;
@@ -236,6 +256,11 @@ procedure tAusgabeDatei.speichereWerte(gitter: tGitter; sDX: extended);
 var i,cnt: longint;
     sX,cX: extended;
 begin
+  if (teilchen>=gitter.felders[gitter.aktuelleFelder].matAnz) then begin
+    pro.schreibe('Teilchen '+inttostr(teilchen+1)+' gibt es nicht, da kann ich auch nichts speichern!',true);
+    pro.destroyall;
+    halt(1);
+  end;
   if (teilchen<0) and (emInhalt=efA) then
     gitter.berechne('A',teilchen);
   if (teilchen>=0) and (matInhalt=mfP) then
@@ -278,10 +303,107 @@ begin
   closeFile(datei);
   if cnt<>0 then begin
     pro.schreibe('Falsche Anzahl an Ortsschritten geschrieben ('+inttostr(cnt)+')!',true);
+    pro.destroyall;
     halt(1);
   end;
 end;
 
+{ tTeilchenSpezies }
+
+constructor tTeilchenSpezies.create;
+begin
+  inherited create;
+  fillchar(ortsGrenzen,sizeof(ortsGrenzen),#0);
+  setlength(ortsGrenzen,0);
+  fillchar(dichteStuecke,sizeof(dichteStuecke),#0);
+  setlength(dichteStuecke,0);
+  nMax:=0;
+  spezifischeLadung:=0;
+end;
+
+destructor tTeilchenSpezies.destroy;
+begin
+  setlength(ortsGrenzen,0);
+  setlength(dichteStuecke,0);
+  inherited destroy;
+end;
+
+procedure tTeilchenSpezies.dump(prot: tProtokollant; prefix: string);
+var
+  i: longint;
+begin
+  prot.schreibe(prefix+'nMax = '+floattostr(nMax)+' * nc');
+  prot.schreibe(prefix+'q/m  = '+floattostr(spezifischeLadung)+' e/me');
+  prot.schreibe(prefix+'n(x):');
+  for i:=0 to length(dichteStuecke)-1 do begin
+    prot.schreibe(prefix+' '+dichteStuecke[i]);
+    if i<length(ortsGrenzen) then
+      prot.schreibe(prefix+'  '+floattostr(ortsGrenzen[i]));
+  end;
+end;
+
+function tTeilchenSpezies.gibDichte(x: extended; kvs: tKnownValues): extended;
+var
+  i: longint;
+begin
+  i:=0;
+  while (i<length(ortsGrenzen)) and (ortsGrenzen[i]<x) do
+    inc(i);
+  kvs.add('x',x);
+  result:=nMax*exprToFloat(false,dichteStuecke[i],kvs,nil);
+  kvs.rem('x');
+end;
+
+procedure tTeilchenSpezies.liesDichteFunktion(rest: string; ifile: tMyStringList; kvs: tKnownValues; teilchen: array of tTeilchenSpezies; prot: tProtokollant);
+var
+  s:     string;
+  ori,i: longint;
+begin
+  if rest='stückweise' then begin
+    setlength(ortsGrenzen,0);
+    setlength(dichteStuecke,0);
+    repeat
+      if not ifile.readln(s) then begin
+        prot.schreibe('Unerwartetes Dateiende in Parameterdatei im Bereich teilchen -> verteilung stückweise!',true);
+        prot.destroyall;
+        halt(1);
+      end;
+
+      if length(dichteStuecke)=length(ortsGrenzen) then begin // neues Funktionsstück wird definiert
+        setlength(dichteStuecke,length(dichteStuecke)+1);
+        dichteStuecke[length(dichteStuecke)-1]:=s;
+        continue;
+      end;
+
+      // kein neues Funktionsstück =>
+      if s='stückweiseEnde' then break; // Ende oder
+      // neue Ortsgrenze
+
+      setlength(ortsGrenzen,length(ortsGrenzen)+1);
+      ortsGrenzen[length(ortsGrenzen)-1]:=exprToFloat(false,s,kvs,nil);
+    until false;
+  end
+  else if startetMit('wie teilchen',rest) then begin
+    ori:=strtoint(rest)-1;
+    if (ori<0) or (ori>=length(teilchen)-1) then begin
+      prot.schreibe('Kann mich nicht auf die Dichte von Teilchen '+inttostr(ori+1)+' beziehen, weil es noch nicht definiert wurde!',true);
+      prot.destroyall;
+      halt(1);
+    end;
+    setlength(ortsGrenzen,length(teilchen[ori].ortsGrenzen));
+    for i:=0 to length(ortsGrenzen)-1 do
+      ortsGrenzen[i]:=teilchen[ori].ortsGrenzen[i];
+    setlength(dichteStuecke,length(teilchen[ori].dichteStuecke));
+    for i:=0 to length(dichteStuecke)-1 do
+      dichteStuecke[i]:=teilchen[ori].dichteStuecke[i];
+  end
+  else begin
+    setlength(ortsGrenzen,0);
+    setlength(dichteStuecke,1);
+    dichteStuecke[0]:=rest;
+  end;
+end;
+
 { tWertePunkt }
 
 constructor tWertePunkt.create(linkerNachbar: tWertePunkt; materien: longint);
@@ -470,27 +592,37 @@ begin
         + fak5 * in5.matWerte[i,maF,true];
 end;
 
+function tWertePunkt.maxDT: extended;
+var
+  i: longint;
+begin
+  result:=1;
+  for i:=0 to length(matWerte)-1 do
+    if (matWerte[i,mfN,true]<0) and
+       (matWerte[i,mfN,false]>0) then
+      result:=min(result,matWerte[i,mfN,false]/matWerte[i,mfN,true]);
+end;
+
 { tFelder }
 
-constructor tFelder.create(groesse: longint; dichten,lichter: tMyStringList; parent: tGitter);
+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:=dichten.count;
+  matAnz:=length(teilchen);
   fillchar(inhalt,sizeof(inhalt),#0);
   setlength(inhalt,groesse+4); // zwei Felder links und rechts extra für Randbedingungen
   inhalt[0]:=tWertePunkt.create(nil,matAnz);
   for i:=1 to length(inhalt)-1 do
     inhalt[i]:=tWertePunkt.create(inhalt[i-1],matAnz);
-  for i:=0 to length(inhalt)-1 do 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 i:=0 to length(inhalt)-1 do
+    for j:=0 to length(teilchen)-1 do begin
+      inhalt[i].matWerte[j,mfN,false]:=teilchen[j].gibDichte(par.xl+(i-2)*par.dX,parent.kvs);
+      inhalt[i].matWerte[j,mfDPsiDX,false]:=0;
+    end;
   for rechts:=false to true do begin
     lichters[rechts]:=tMyStringlist.create(nil,'');
     lichters[rechts].text:=lichter.text;
@@ -572,9 +704,18 @@ begin
     inhalt[i].liKo(in1.inhalt[i],in2.inhalt[i],in3.inhalt[i],in4.inhalt[i],in5.inhalt[i],fak2,fak3,fak4,fak5);
 end;
 
+function tFelder.maxDT: extended;
+var
+  i: longint;
+begin
+  result:=inhalt[0].maxDT;
+  for i:=1 to length(inhalt)-1 do
+    result:=min(result,inhalt[i].maxDT);
+end;
+
 { tGitter }
 
-constructor tGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; dichten, lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string);
+constructor tGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string);
 var
   i: longint;
 begin
@@ -597,7 +738,7 @@ begin
   xl:=dX/2;
 
   for i:=0 to length(felders)-1 do
-    felders[i]:=tFelder.create(size,dichten,lichter,self);
+    felders[i]:=tFelder.create(size,teilchen,lichter,self);
 
   aktuelleFelder:=0;
   t:=0;
@@ -614,16 +755,45 @@ begin
   inherited destroy;
 end;
 
-procedure tGitter.iteriereSchritt(dT,iDT: extended);
-var i:     longint;
+procedure tGitter.iteriereSchritt(var dT: extended);
+var i:       longint;
     {$IFDEF Zeitschrittueberwachung}
-    pro:   tProtokollant;
-    j:     longint;
+    pro:     tProtokollant;
+    j:       longint;
     {$ENDIF}
+    iDT,mDT: extended;
 begin
-  t:=t+dT;
-  kvs.add('t',t);
-  felders[aktuelleFelder].berechneAbleitungen(dT,dX,iDT,iDX,pDNMax); // y' = y'(t,y(t))
+  kvs.add('t',t+dT);
+  diffusionsTermAnpassen(prot);
+
+  repeat
+    iDT:=1/dT;
+    felders[aktuelleFelder].berechneAbleitungen(dT,dX,iDT,iDX,pDNMax); // y' = y'(t,y(t))
+
+    mDT:=felders[aktuelleFelder].maxDT;
+      // das maximale dT, welches bei den Momentanen Ableitungen nicht zu
+      // unphysikalischen Effekten (z.B. negativen Dichten) führt
+
+    if mDT<dT*2 then begin // beabsichtigter Zeitschritt ist zu groß
+      dT:=mDT/4;           // machen wir ihn kleiner!
+      {$IFDEF Zeitschrittueberwachung}
+      pro:=tProtokollant.create(prot,'iteriereSchritt');
+      pro.schreibe('dT -> '+floattostr(dT));
+      if dT<1E-30 then begin
+        pro.schreibe('Zeitschritt geht gegen Null (ist bereits '+floattostr(dT)+' < 10^-30) - irgendwas scheint grundsätzlich kaputt zu sein!',true);
+        pro.destroyall;
+        halt(1);
+      end;
+      pro.free;
+      {$ENDIF}
+      continue;
+    end;
+
+    break;
+  until false;
+
+  if dT<mDT/1000 then
+    dT:=min(1,mDT)/1000;
 
   case zeitverfahren of
     zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt
@@ -680,31 +850,38 @@ begin
           halt(1);
         end;
   {$ENDIF}
+  t:=t+dT;
 end;
 
-function tGitter.gibErhaltungsgroessen: string;
+function tGitter.dumpErhaltungsgroessen: boolean;
 var i,j: integer;
     ns:  tExtendedArray;
     pro: tProtokollant;
+    s:   string;
 begin
-  pro:=tProtokollant.create(prot,'gibErhaltungsgroessen');
+  pro:=tProtokollant.create(prot,'dumpErhaltungsgroessen');
   setlength(ns,felders[aktuelleFelder].matAnz);
-  result:='';
+
+  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
-    result:=result + ' n['+inttostr(i)+']='+floattostr(ns[i]);
+    ns[i]:=ns[i]*dX;
+    s:=s+ ' 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)+')');
+      result:=false;
+      pro.schreibe(' n['+inttostr(i)+'] > 1000, es scheinen sehr viele neue Teilchen entstanden zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')',true);
     end;
     {$ENDIF}
   end;
-  delete(result,1,1);
+  delete(s,1,1);
+  prot.schreibe(s);
   pro.free;
 end;
 
@@ -763,6 +940,7 @@ begin
     exit;
   end;
   prot.schreibe('Kann '''+was+''' nicht berechnen, weil ich es nicht verstehe!',true);
+  prot.destroyall;
   halt;
 end;
 
@@ -776,7 +954,8 @@ var
   deltaX,breite,pDNMax:                    extended;
   i:                                       longint;
   kvs:                                     tKnownValues;
-  dichten,lichter:                         tMyStringlist;
+  teilchen:                                array of tTeilchenSpezies;
+  lichter:                                 tMyStringlist;
   pro:                                     tProtokollant;
 begin
   inherited create;
@@ -787,6 +966,7 @@ begin
   ifile.loadfromfile(inName);
   if not ifile.unfoldMacros then begin
     prot.schreibe('Fehlerhafte Macros in Parameterdatei '''+inName+'''!',true);
+    prot.destroyall;
     halt(1);
   end;
 
@@ -797,6 +977,8 @@ begin
   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;
@@ -814,11 +996,12 @@ begin
   lichter:=tMyStringlist.create(prot,'lichter');
 
   // Standardeinstellungen Breich 'teilchen...'
-  dichten:=tMyStringlist.create(prot,'teilchen');
+  setlength(teilchen,0);
 
   repeat
     if not ifile.readln(s) then begin
       prot.schreibe('Unerwartetes Dateiende in '''+inName+'''!',true);
+      prot.destroyall;
       halt(1);
     end;
 
@@ -829,6 +1012,7 @@ 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;
@@ -840,25 +1024,25 @@ begin
           Zeitverfahren:=zfEulerVorwaerts;
           continue;
         end;
-        if startetMit('ortsschritt',s) then begin
+        if startetMit('ortsschritt ',s) then begin
           deltaX:=exprtofloat(false,s,kvs,nil);
           kvs.add('dX',deltaX);
           continue;
         end;
-        if startetMit('zeitschritt',s) then begin
+        if startetMit('zeitschritt ',s) then begin
           dT:=exprtofloat(false,s,kvs,nil);
           kvs.add('dT',dT);
           continue;
         end;
-        if startetMit('diffusionsterm',s) then begin
+        if startetMit('diffusionsterm ',s) then begin
           pDNMax:=exprtofloat(false,s,kvs,nil);
           continue;
         end;
-        if startetMit('zeit',s) then begin
+        if startetMit('zeit ',s) then begin
           tEnde:=exprtofloat(false,s,kvs,nil);
           continue;
         end;
-        if startetMit('breite',s) then begin
+        if startetMit('breite ',s) then begin
           breite:=exprtofloat(false,s,kvs,nil);
           continue;
         end;
@@ -871,6 +1055,7 @@ begin
           continue;
         end;
         prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich allgemein!',true);
+        prot.destroyall;
         halt(1);
       until false;
       continue;
@@ -880,26 +1065,27 @@ 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
+        if startetMit('suffix ',s) then begin
           aSuffix:=s;
           continue;
         end;
-        if startetMit('prefix',s) then begin
+        if startetMit('prefix ',s) then begin
           aPrefix:=s;
           continue;
         end;
-        if startetMit('zeitschritt',s) then begin
+        if startetMit('zeitschritt ',s) then begin
           sDT:=exprtofloat(false,s,kvs,nil);
           continue;
         end;
-        if startetMit('ortsschritt',s) then begin
+        if startetMit('ortsschritt ',s) then begin
           sDX:=exprtofloat(false,s,kvs,nil);
           continue;
         end;
-        if startetMit('felder',s) then begin
+        if startetMit('felder ',s) then begin
           s:=s+',';
           while s<>'' do begin
             t:=erstesArgument(s,',');
@@ -909,6 +1095,7 @@ begin
           continue;
         end;
         prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich ausgaben!',true);
+        prot.destroyall;
         halt(1);
       until false;
       continue;
@@ -924,10 +1111,49 @@ begin
         continue;
       end;
       prot.schreibe('Licht kann nicht von '''+erstesArgument(s)+''' kommen!',true);
+      prot.destroyall;
       halt(1);
     end;
 
+    if startetMit('teilchen',s) then begin
+      if (s<>'') and (strtoint(s)<>length(teilchen)+1) then begin
+        prot.schreibe('Ich erwarte die Teilchen beginnend bei 1 der Reihe nach in Parameterdatei '''+inName+'''!',true);
+        prot.destroyall;
+        halt(1);
+      end;
+      setlength(teilchen,length(teilchen)+1);
+      teilchen[length(teilchen)-1]:=tTeilchenSpezies.create;
+      repeat
+        if not ifile.readln(s) then begin
+          prot.schreibe('Unerwartetes Dateiende in Parameterdatei '''+inName+''' im Bereich teilchen!',true);
+          prot.destroyall;
+          halt(1);
+        end;
+        if s='teilchen'+inttostr(length(teilchen))+'Ende' then break;
+        if startetMit('spezifische Ladung ',s) then begin
+          teilchen[length(teilchen)-1].spezifischeLadung:=exprToFloat(false,s,kvs,nil);
+          kvs.add('qm'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].spezifischeLadung);
+          continue;
+        end;
+        if startetMit('maximaldichte ',s) then begin
+          teilchen[length(teilchen)-1].nMax:=exprToFloat(false,s,kvs,nil);
+          kvs.add('nmax'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].nMax);
+          continue;
+        end;
+        if startetMit('verteilung ',s) then begin
+          teilchen[length(teilchen)-1].liesDichteFunktion(s,ifile,kvs,teilchen,prot);
+          continue;
+        end;
+
+        prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+''' im Bereich teilchen!',true);
+        prot.destroyall;
+        halt(1);
+      until false;
+      continue;
+    end;
+
     prot.schreibe('Unbekannter Befehl '''+s+''' in Parameterdatei '''+inName+'''!',true);
+    prot.destroyall;
     halt(1);
   until false;
 
@@ -935,6 +1161,7 @@ begin
 
   if length(ausgabedateien)=0 then begin
     prot.schreibe('Du solltest irgendetwas abspeichern lassen!',true);
+    prot.destroyall;
     halt(1);
   end;
 
@@ -962,9 +1189,10 @@ begin
   pro.schreibe('pDNMax = '+floattostr(pDNMax));
   pro.schreibe('bekannte Werte:');
   kvs.dump(pro,'  ');
-  if dichten.count>0 then begin
-    pro.schreibe('dichten:');
-    dichten.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:');
@@ -980,8 +1208,10 @@ begin
   end;
   pro.free;
 
-  gitter:=tGitter.create(round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,dichten,lichter,zeitverfahren,prot,'gitter');
-  dichten.free;
+  gitter:=tGitter.create(round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,teilchen,lichter,zeitverfahren,prot,'gitter');
+  for i:=0 to length(teilchen)-1 do
+    teilchen[i].free;
+  setlength(teilchen,0);
   lichter.free;
 end;
 
@@ -1000,7 +1230,7 @@ begin
 
   zeitPhysik:=zeitPhysik-now;
   if errorCode<2 then
-    gitter.iteriereSchritt(dT,1/dT);
+    gitter.iteriereSchritt(dT);
   zeitPhysik:=zeitPhysik+now;
   zeitDatei:=zeitDatei-now;
   while (gitter.t>=sT) and (sT<tEnde) do begin
@@ -1012,12 +1242,15 @@ begin
 
   if fortschrittsAnzeige then begin
     if (floor(100*Gitter.t/tEnde) < floor(100*(Gitter.t+dT)/tEnde)) then begin
-      Prot.schreibe(inttostr(round(100*Gitter.t/tEnde))+'% (t='+floattostr(Gitter.t)+'T)',true);
-      Prot.schreibe(timetostr(now-start)+' ('+floattostr(zeitPhysik/max(1e-11,zeitPhysik+zeitDatei))+')',true);
-      Prot.schreibe('ETA: '+timetostr((now-start)*(tEnde-Gitter.t)/max(Gitter.t,dT)),true);
-      Prot.schreibe('aktueller Zeitschritt: '+floattostr(dT)+'T',true);
+      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);
+        if not gitter.dumpErhaltungsgroessen then begin
+          errorCode:=2;
+          exit;
+        end;
     end;
   end;
 
diff --git a/Plasmapropagation b/Plasmapropagation
index 3501611..561167f 100755
--- a/Plasmapropagation
+++ b/Plasmapropagation
diff --git a/Plasmapropagation.lps b/Plasmapropagation.lps
index cbe3abd..b5afdb3 100644
--- a/Plasmapropagation.lps
+++ b/Plasmapropagation.lps
@@ -7,20 +7,22 @@
       <Unit0>
         <Filename Value="Plasmapropagation.lpr"/>
         <IsPartOfProject Value="True"/>
+        <TopLine Value="25"/>
         <CursorPos X="34" Y="11"/>
-        <UsageCount Value="149"/>
+        <UsageCount Value="153"/>
         <Loaded Value="True"/>
       </Unit0>
       <Unit1>
         <Filename Value="Physikunit.pas"/>
         <IsPartOfProject Value="True"/>
+        <IsVisibleTab Value="True"/>
         <EditorIndex Value="1"/>
-        <TopLine Value="438"/>
-        <CursorPos Y="538"/>
-        <FoldState Value=" T3k60-4 piYkM0O]9j8jg0M[943jb0J5]9SkI0O3 pnBo2051g"/>
-        <UsageCount Value="90"/>
+        <TopLine Value="760"/>
+        <CursorPos X="45" Y="781"/>
+        <FoldState Value=" T3kO0-4 pigkU095 pjDja083]90jG0M[I4Bk40J513[c4nL0~E2%"/>
+        <UsageCount Value="94"/>
         <Bookmarks Count="1">
-          <Item0 Y="627"/>
+          <Item0 Y="790"/>
         </Bookmarks>
         <Loaded Value="True"/>
       </Unit1>
@@ -30,33 +32,34 @@
         <EditorIndex Value="2"/>
         <TopLine Value="20"/>
         <CursorPos X="22" Y="31"/>
-        <UsageCount Value="52"/>
+        <UsageCount Value="56"/>
         <Loaded Value="True"/>
       </Unit2>
       <Unit3>
         <Filename Value="input.plap"/>
         <IsPartOfProject Value="True"/>
         <EditorIndex Value="4"/>
-        <CursorPos Y="22"/>
-        <UsageCount Value="51"/>
+        <TopLine Value="9"/>
+        <CursorPos X="30" Y="43"/>
+        <UsageCount Value="55"/>
         <Loaded Value="True"/>
         <DefaultSyntaxHighlighter Value="None"/>
       </Unit3>
       <Unit4>
         <Filename Value="input.epost"/>
-        <IsVisibleTab Value="True"/>
         <EditorIndex Value="5"/>
-        <CursorPos X="27" Y="21"/>
-        <UsageCount Value="58"/>
+        <CursorPos X="32" Y="21"/>
+        <UsageCount Value="60"/>
         <Loaded Value="True"/>
         <DefaultSyntaxHighlighter Value="None"/>
       </Unit4>
       <Unit5>
         <Filename Value="../units/matheunit.pas"/>
         <EditorIndex Value="3"/>
-        <TopLine Value="316"/>
-        <CursorPos Y="356"/>
-        <UsageCount Value="19"/>
+        <TopLine Value="83"/>
+        <CursorPos Y="53"/>
+        <FoldState Value=" T3i905B pj0jV034 piaj60U2-"/>
+        <UsageCount Value="21"/>
         <Loaded Value="True"/>
       </Unit5>
       <Unit6>
@@ -105,122 +108,122 @@
     <JumpHistory Count="30" HistoryIndex="29">
       <Position1>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="533" Column="45" TopLine="502"/>
+        <Caret Line="112" Column="37" TopLine="92"/>
       </Position1>
       <Position2>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="679" Column="54" TopLine="522"/>
+        <Caret Line="736" TopLine="686"/>
       </Position2>
       <Position3>
-        <Filename Value="../units/protokollunit.pas"/>
-        <Caret Line="28" Column="23"/>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="1195" Column="41" TopLine="891"/>
       </Position3>
       <Position4>
-        <Filename Value="../units/protokollunit.pas"/>
-        <Caret Line="29" Column="23"/>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="121" Column="35" TopLine="101"/>
       </Position4>
       <Position5>
-        <Filename Value="../units/protokollunit.pas"/>
-        <Caret Line="76" Column="15" TopLine="58"/>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="806" Column="58" TopLine="678"/>
       </Position5>
       <Position6>
-        <Filename Value="../units/protokollunit.pas"/>
-        <Caret Line="114" TopLine="77"/>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="807" Column="39" TopLine="678"/>
       </Position6>
       <Position7>
-        <Filename Value="../units/protokollunit.pas"/>
-        <Caret Line="29" TopLine="29"/>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="832" Column="36" TopLine="678"/>
       </Position7>
       <Position8>
-        <Filename Value="../units/protokollunit.pas"/>
-        <Caret Line="72" Column="23" TopLine="52"/>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="1195" Column="31" TopLine="862"/>
       </Position8>
       <Position9>
-        <Filename Value="../units/protokollunit.pas"/>
-        <Caret Line="31" Column="24"/>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="121" Column="35" TopLine="101"/>
       </Position9>
       <Position10>
-        <Filename Value="../units/protokollunit.pas"/>
-        <Caret Line="114" Column="25" TopLine="80"/>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="813" Column="76" TopLine="678"/>
       </Position10>
       <Position11>
-        <Filename Value="Plasmapropagation.lpr"/>
-        <Caret Line="35" Column="63"/>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="121" Column="37" TopLine="98"/>
       </Position11>
       <Position12>
-        <Filename Value="Plasmapropagation.lpr"/>
-        <Caret Line="26" Column="11" TopLine="15"/>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="811" Column="17" TopLine="736"/>
       </Position12>
       <Position13>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="509" TopLine="369"/>
+        <Caret Line="827" Column="144" TopLine="736"/>
       </Position13>
       <Position14>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="11" Column="84"/>
+        <Caret Line="1198" Column="20" TopLine="894"/>
       </Position14>
       <Position15>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1010" Column="69" TopLine="990"/>
       </Position15>
       <Position16>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1017" TopLine="764"/>
+        <Caret Line="827" Column="40" TopLine="735"/>
       </Position16>
       <Position17>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="937" Column="40" TopLine="915"/>
+        <Caret Line="1200" Column="5" TopLine="894"/>
       </Position17>
       <Position18>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="262" Column="50" TopLine="242"/>
+        <Caret Line="820" Column="52" TopLine="686"/>
       </Position18>
       <Position19>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="55" Column="53" TopLine="23"/>
+        <Caret Line="121" Column="45" TopLine="103"/>
       </Position19>
       <Position20>
         <Filename Value="Physikunit.pas"/>
+        <Caret Line="828" Column="13" TopLine="736"/>
       </Position20>
       <Position21>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="15" Column="57"/>
+        <Caret Line="807" TopLine="736"/>
       </Position21>
       <Position22>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="410" Column="48" TopLine="382"/>
+        <Caret Line="1175" TopLine="782"/>
       </Position22>
       <Position23>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="528" Column="149" TopLine="415"/>
+        <Caret Line="120" Column="33" TopLine="113"/>
       </Position23>
       <Position24>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="511"/>
+        <Caret Line="736" Column="35" TopLine="695"/>
       </Position24>
       <Position25>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="529" Column="35" TopLine="473"/>
+        <Caret Line="120" Column="35" TopLine="87"/>
       </Position25>
       <Position26>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="511" Column="19" TopLine="475"/>
+        <Caret Line="802" Column="51" TopLine="790"/>
       </Position26>
       <Position27>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="75" Column="31" TopLine="61"/>
+        <Caret Line="6"/>
       </Position27>
       <Position28>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="534" Column="61" TopLine="472"/>
+        <Caret Line="760" Column="121" TopLine="726"/>
       </Position28>
       <Position29>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="64" Column="35" TopLine="47"/>
+        <Caret Line="736" TopLine="558"/>
       </Position29>
       <Position30>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="369" Column="43" TopLine="344"/>
+        <Caret Line="701" Column="41" TopLine="465"/>
       </Position30>
     </JumpHistory>
   </ProjectSession>
diff --git a/input.epost b/input.epost
index b85c634..5a448af 100644
--- a/input.epost
+++ b/input.epost
@@ -18,7 +18,7 @@ Ende
 Threadanzahl:     10
 
 # !Schleife: $Feld: dDPHIDX dDAYDT dAY dAX
-!Schleife: $Feld: AY DAYDT
+!Schleife: $Feld: AY DAYDT N1 N2
 
 Daten einlesen
   Genauigkeit:     extended
diff --git a/input.plap b/input.plap
index 34fcf17..5e98ccc 100644
--- a/input.plap
+++ b/input.plap
@@ -5,48 +5,45 @@ allgemein
   ortsschritt    10^-2 * λ
   zeitschritt    10^-2 * T
   zeit           20 * T
-  breite         5 * λ
+!setze $breite: (5 * λ)
+  breite         $breite
   mit Fortschrittsanzeige
 allgemeinEnde
 
 ausgaben
   prefix         /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/
   suffix         _test.dat
-  felder         AX,AY,dAYDT,dPhiDX
-#  felder         AX,AY,dAYDT,N,dPhiDX,dPsiDX,Gamma
+  felder         AX,AY,dAYDT,dPhiDX,N1,dPsiDX1,Gamma1,N2,dPsiDX2,Gamma2
 ausgabenEnde
 
 !setze $tFwhm:   (2.5 * T)
 !setze $tMitte:  (1 * T)
 
-licht von links    0.5 * 2^(-2*((t-$tMitte)/$tFwhm)^2) * sin(ω*t)
-
-Dateiende
+licht von links      0.5 * 2^(-2*((t-$tMitte)/$tFwhm)^2) * sin(ω*t)
 
 teilchen1
-  ladung         -q
-  masse          me
-  maximaldichte  0
-  beweglich
-  !setze $profilbreite: 5 λ
-  !setze $randbreite: .1 λ
+  spezifische Ladung -q/me
+  maximaldichte      10
+!setze $profilbreite: (4 * λ)
+!setze $randbreite:   (0.1 * λ)
   verteilung stückweise
     0
-      (breite-$profilbreite)/2 - $randbreite
-    sin((x - (breite-$profilbreite)/2 - $randbreite)*π/2/$randbreite)^2
-      (breite-$profilbreite)/2
+      ($breite-$profilbreite)/2 - $randbreite
+    sin((x - ($breite-$profilbreite)/2 - $randbreite)*π/2/$randbreite)^2
+      ($breite-$profilbreite)/2
     1
-      (breite+$profilbreite)/2
-    sin((x - (breite+$profilbreite)/2 + $randbreite)*π/2/$randbreite)^2
-      (breite+$profilbreite)/2 + $randbreite
+      ($breite+$profilbreite)/2
+    sin((x - ($breite+$profilbreite)/2 + $randbreite)*π/2/$randbreite)^2
+      ($breite+$profilbreite)/2 + $randbreite
     0
   stückweiseEnde
 teilchen1Ende
 
 teilchen2
-  ladung         q
-  maximaldichte  nmax1
-  unbeweglich
+  spezifische Ladung (q/me)/2
+  maximaldichte      nmax1*2
   verteilung wie teilchen1
 teilchen2Ende
 
+Dateiende
+