lauffaehig, Materie haut aber ab.

5 files changed, 304 insertions, 144 deletions

diff --git a/Physikunit.pas b/Physikunit.pas
index 142ee4b..ab99170 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -8,6 +8,7 @@ unit Physikunit;
 
 {$DEFINE Zeitschrittueberwachung}
 {$DEFINE Dichteueberwachung}
+{ $DEFINE negativeDichteueberwachung}
 
 interface
 
@@ -25,11 +26,20 @@ type
   tMaterieFeldInhalt = (
     mfN,mfDPsiDX,       // Teilchendichte, Geschwindigkeitsdichte
     mfP,mfPX,mfPY,mfPZ, // Impuls
-    mfGamma,mfIGamma   // rel. Gamma-Faktor
+    mfGamma,mfIGamma    // rel. Gamma-Faktor
   );
 
+const
+  erstesEMFMitAbleitung  =  efAX;
+  letztesEMFMitAbleitung =  efDAZDT;
+  erstesMatFMitAbleitung  = mfN;
+  letztesMatFMitAbleitung = mfDPsiDX;
+
+type
+
   tGitter = class;
   tFelder = class;
+  tSimulation = class;
 
   { tAusgabeDatei }
 
@@ -81,15 +91,16 @@ type
       // 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!
+      // A, dPhiDX, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier!
     lN,rN:    tWertePunkt;
 
     constructor create(linkerNachbar: tWertePunkt; derBesitzer: tFelder);
     destructor destroy; override;
     procedure berechneGammaUndP;
+    procedure berechneDPhiDX(dX: extended);
     procedure berechneNAbleitung(iDX,pDNMax: extended); overload;
     procedure berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean); overload;
-    procedure berechneAbleitungen(dX,iDX: extended);
+    procedure berechneAbleitungen(iDX: extended);
     procedure liKo(in1,in2: tWertePunkt; fak2: extended); overload;                            // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
     procedure liKo(in1,in2,in3: tWertePunkt; fak2,fak3: extended); overload;
     procedure liKo(in1,in2,in3,in4: tWertePunkt; fak2,fak3,fak4: extended); overload;
@@ -112,6 +123,8 @@ type
     procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31: extended); overload;
     procedure liKo(in1,in2,in3,in4,in5,in6,in7,in8,in9,in10,in11,in12,in13,in14,in15,in16,in17,in18,in19,in20,in21,in22,in23,in24,in25,in26,in27,in28,in29,in30,in31,in32: tWertePunkt; fak2,fak3,fak4,fak5,fak6,fak7,fak8,fak9,fak10,fak11,fak12,fak13,fak14,fak15,fak16,fak17,fak18,fak19,fak20,fak21,fak22,fak23,fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31,fak32: extended); overload;
     function maxDT: extended;
+    function dump(mitNachbarn: boolean): string;
+    procedure nichtnegativieren;
   end;
 
   { TFelder }
@@ -119,7 +132,8 @@ type
   tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen
   private
     matAnz:       longint;
-    spezLadungen: tExtendedArray;
+    spezLadungen,
+    massen:       tExtendedArray;
     lichters:     array[boolean] of tMyStringlist;
     procedure setzeRaender(iDX: extended);
   public
@@ -157,6 +171,7 @@ type
 
   tGitter = class(tObject)
   private
+    besitzer:                     tSimulation;
     prot:                         tProtokollant;
     zeitverfahren:                tZeitverfahren;
     pDNMaxDynamisch,abbruch:      boolean;
@@ -170,7 +185,7 @@ 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; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string);
+    constructor create(derBesitzer: tSimulation; size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string);
     destructor destroy; override;
     procedure iteriereSchritt(var dT: extended; dTMin: extended);
     function dumpErhaltungsgroessen: boolean;
@@ -429,7 +444,7 @@ var
   i: longint;
 begin
   prot.schreibe(prefix+'nMax = '+floattostr(nMax)+' * nc');
-  prot.schreibe(prefix+'q/m  = '+floattostr(spezifischeLadung)+' e/me');
+  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]);
@@ -557,6 +572,29 @@ begin
   end;
 end;
 
+procedure tWertePunkt.berechneDPhiDX(dX: extended);
+var
+  i: longint;
+begin
+  // d2Phi/dx2 = rho  ... hier muss integriert werden:
+  // dPhi/dx = dPhi/dx(x- deltax) + <rho> * deltax
+
+  if assigned(lN) then begin
+    emWerte[efDPhiDX,false]:=lN.emWerte[efDPhiDX,false];
+    for i:=0 to length(matWerte)-1 do
+      emWerte[efDPhiDX,false]:=
+        emWerte[efDPhiDX,false]
+        + besitzer.spezLadungen[i] * (lN.matWerte[i,mfN,false] + matWerte[i,mfN,false]) * dX/2;
+  end
+  else begin
+    emWerte[efDPhiDX,false]:=0;
+    for i:=0 to length(matWerte)-1 do
+      emWerte[efDPhiDX,false]:=
+        emWerte[efDPhiDX,false]
+        + besitzer.spezLadungen[i] * matWerte[i,mfN,false] * dX/2;
+  end;
+end;
+
 procedure tWertePunkt.berechneNAbleitung(iDX,pDNMax: extended); // Zeitableitung von n berechnen
 var
   i: longint;
@@ -602,7 +640,7 @@ begin
   end
 end;
 
-procedure tWertePunkt.berechneAbleitungen(dX,iDX: extended); // Zeitableitungen berechnen
+procedure tWertePunkt.berechneAbleitungen(iDX: extended); // Zeitableitungen berechnen
 var
   i: longint;
 begin
@@ -610,16 +648,7 @@ begin
   for i:=0 to length(matWerte)-1 do
     matWerte[i,mfDPsiDX,true]:=
       besitzer.spezLadungen[i] * emWerte[efDPhiDX,false]
-        - (rN.matWerte[i,mfGamma,false] - lN.matWerte[i,mfGamma,false]) * iDX/2;
-
-  // d3Phi/dx2dt = dRho/dt  ... hier muss integriert werden:
-  // d2Phi/dxdt = d2Phi/dxdt(x- deltax) + <dRho/dt> * deltax
-  emWerte[efDPhiDX,true]:=
-    lN.emWerte[efDPhiDX,true];
-  for i:=0 to length(matWerte)-1 do
-    emWerte[efDPhiDX,true]:=
-      emWerte[efDPhiDX,true]
-      + besitzer.spezLadungen[i] * (lN.matWerte[i,mfN,true] + matWerte[i,mfN,true])*dX/2;
+;//        - (rN.matWerte[i,mfGamma,false] - lN.matWerte[i,mfGamma,false]) * iDX/2;
 
   // d2A/dt2 = Laplace(A) - Nabla(phi) ...
   emWerte[efDAXDT,true]:=
@@ -659,6 +688,73 @@ begin
       result:=min(result,-matWerte[i,mfN,false]/matWerte[i,mfN,true]);
 end;
 
+function tWertePunkt.dump(mitNachbarn: boolean): string;
+var
+  i: longint;
+begin
+  if mitNachbarn and assigned(lN) then result:=lN.dump(false)+' << '
+  else result:='';
+
+  for i:=0 to length(matWerte)-1 do
+    result:=
+      result+
+      'N['+  inttostr(i)+'] = '+floattostr(matWerte[i,mfN,false])+' '+
+      'N''['+inttostr(i)+'] = '+floattostr(matWerte[i,mfN,true])+' ';
+
+  result:=result+'maxDT: '+floattostr(maxDT);
+
+  if mitNachbarn and assigned(rN) then result:=result+' >> '+rN.dump(false);
+end;
+
+procedure tWertePunkt.nichtnegativieren; // Dichten nicht negativ machen
+var
+  i:       longint;
+  defizit: extended;
+  li,re:   tWertePunkt;
+  pro:     tProtokollant;
+begin
+  for i:=0 to length(matWerte)-1 do
+    if matWerte[i,mfN,false]<0 then begin
+      defizit:=-matWerte[i,mfN,false];
+      matWerte[i,mfN,false]:=0;
+      li:=self.lN;
+      re:=self.rN;
+      repeat
+        if assigned(li) then begin
+          if li.matWerte[i,mfN,false]>0 then begin
+            if defizit>=li.matWerte[i,mfN,false] then begin
+              defizit:=defizit-li.matWerte[i,mfN,false];
+              li.matWerte[i,mfN,false]:=0;
+            end
+            else begin
+              li.matWerte[i,mfN,false]:=li.matWerte[i,mfN,false]-defizit;
+              defizit:=0;;
+            end;
+          end;
+          li:=li.lN;
+        end;
+        if assigned(re) then begin
+          if re.matWerte[i,mfN,false]>0 then begin
+            if defizit>=re.matWerte[i,mfN,false] then begin
+              defizit:=defizit-re.matWerte[i,mfN,false];
+              re.matWerte[i,mfN,false]:=0;
+            end
+            else begin
+              re.matWerte[i,mfN,false]:=re.matWerte[i,mfN,false]-defizit;
+              defizit:=0;;
+            end;
+          end;
+          re:=re.rN;
+        end;
+      until (defizit<=0) or not (assigned(li) or assigned(re));
+      if defizit>0 then begin
+        pro:=tProtokollant.create(besitzer.par.prot,'WertePunkt.nichtnegativieren');
+        pro.schreibe('Ich konnte ein Teilchendefizit nicht auffüllen!',true);
+        besitzer.par.abbrechen;
+      end;
+    end;
+end;
+
 { linearkombinationsspezifische Methoden von tWertePunkt und tFelder }
 
 {$INCLUDE linearkombination.inc}
@@ -722,10 +818,15 @@ begin
   inhalt[0]:=tWertePunkt.create(nil,self);
   for i:=1 to length(inhalt)-1 do
     inhalt[i]:=tWertePunkt.create(inhalt[i-1],self);
+  fillchar(massen,sizeof(massen),#0);
+  setlength(massen,length(teilchen));
+  for i:=0 to length(massen)-1 do
+    massen[i]:=0;
   for i:=0 to length(inhalt)-1 do
     for j:=0 to length(teilchen)-1 do begin
       inhalt[i].matWerte[j,mfN,false]:=teilchen[j].gibDichte(par.xl+(i-2)*par.dX,parent.kvs);
       inhalt[i].matWerte[j,mfDPsiDX,false]:=0;
+      massen[j]:=massen[j]+inhalt[i].matWerte[j,mfN,false]*par.dX;
     end;
   for rechts:=false to true do begin
     lichters[rechts]:=tMyStringlist.create(nil,'');
@@ -778,6 +879,8 @@ var
 begin
   for i:=0 to length(inhalt)-1 do
     inhalt[i].berechneGammaUndP;
+  for i:=0 to length(inhalt)-1 do
+    inhalt[i].berechneDPhiDX(dX);
   for i:=0 to matAnz-1 do begin // Teilchen werden reflektiert
     inhalt[1].matWerte[i,mfPX,false]:=
        abs(inhalt[1].matWerte[i,mfPX,false]);
@@ -794,7 +897,7 @@ begin
   setzeRaender(iDX);
 
   for i:=1 to length(inhalt)-2 do
-    inhalt[i].berechneAbleitungen(dX,iDX);
+    inhalt[i].berechneAbleitungen(iDX);
 end;
 
 function tFelder.maxDT: extended;
@@ -808,12 +911,13 @@ end;
 
 { tGitter }
 
-constructor tGitter.create(size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string);
+constructor tGitter.create(derBesitzer: tSimulation; size: longint; deltaT,deltaX,pDNMa: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string);
 var
   i: longint;
 begin
   inherited create;
   abbruch:=false;
+  besitzer:=derBesitzer;
   zeitverfahren:=zv;
   kvs:=bekannteWerte;
   prot:=tProtokollant.create(protokollant,name);
@@ -859,17 +963,21 @@ end;
 
 procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant);
 var
-  i,j:    longint;
-  curMax: extended;
+  i,j,k:      longint;
+  curMax,tmp: extended;
 begin
   curMax:=0;
+  k:=0;
   for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do
-    for j:=0 to felders[aktuelleFelder].matAnz-1 do
-      curMax:=
-        max(curMax,
-            abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPx,false]*
-                (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
-                max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-100)));
+    for j:=0 to felders[aktuelleFelder].matAnz-1 do begin
+      tmp:=abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPx,false]*
+               (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
+               max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-100));
+      if tmp>curMax then begin
+        curMax:=tmp;
+        k:=i;
+      end;
+    end;
   if curMax > pDNMax then begin
     if pDNMaxDynamisch then begin
       pDNMax:=2*curMax;
@@ -878,10 +986,10 @@ begin
     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);
+        'T, x='+floattostr(xL+k*dX)+'), die numerische Stabilität ist nicht mehr gewährleistet, Simulation wird abgebrochen!',true);
       pro.schreibe(' Lösung: größeren Diffusionsterm wählen',true);
       pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern',true);
-      pDNMax:=0; // Term komplett abschalten
+      abbrechen;
     end;
   end;
 end;
@@ -892,18 +1000,17 @@ begin
     // das maximale dT, welches bei den momentanen Ableitungen nicht zu
     // unphysikalischen Effekten (z.B. negativen Dichten) führt
 
-  result:=dTMax<dT*2;
+  result:=(dTMax<dT) and (dT>1.1*dTMin);
   if result then begin // beabsichtigter Zeitschritt ist zu groß,
     dT:=dTMax/4;       // dann machen wir ihn doch kleiner!
+    if dT<dTMin then
+      dT:=dTMin
     {$IFDEF Zeitschrittueberwachung}
-    prot.schreibe('pruefeMaxDT: dT -> '+floattostr(dT)+' (t = '+floattostr(t)+')');
-    if dT<dTMin then begin
-      prot.schreibe('pruefeMaxDT: Zeitschritt geht gegen Null (ist bereits '+floattostr(dT)+' < '+floattostr(dTMin)+') - irgendwas scheint grundsätzlich kaputt zu sein!',true);
-      abbrechen;
-    end;
-    prot.spuelen;
-    {$ENDIF}
-    exit;
+    else begin
+      prot.schreibe('pruefeMaxDT: dT -> '+floattostr(dT)+' (t = '+floattostr(t)+')');
+      prot.spuelen;
+    end
+    {$ENDIF};
   end;
 end;
 
@@ -928,11 +1035,8 @@ begin
 end;
 
 procedure tGitter.iteriereSchritt(var dT: extended; dTMin: extended);
-var i:       longint;
-    {$IFDEF Zeitschrittueberwachung}
-    j:       longint;
-    {$ENDIF}
-    {$IFDEF Dichteueberwachung}
+var {$IFDEF negativeDichteueberwachung}
+    i,j:     longint;
     pro:     tProtokollant;
     {$ENDIF}
     dTMax:   extended;
@@ -973,17 +1077,18 @@ begin
 
     break;
   until abbruch;
-                        prot.schreibe('dT = '+floattostr(dT));
+
   aktuelleFelder:=1-aktuelleFelder;
-  {$IFDEF Dichteueberwachung}
+  {$IFDEF negativeDichteueberwachung}
   for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do
     for j:=0 to felders[aktuelleFelder].matAnz-1 do
-        if felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false]<0 then
+        if felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false]<0 then begin
+          pro:=tProtokollant.create(prot,'iteriereSchritt');
+          pro.schreibe('n<0 bei:',true);
+          pro.schreibe(' t    = ' +floattostr(t)+',',true);
+          pro.schreibe(' i    = ' +inttostr(i)+' (x = '+floattostr(xl+i*dX)+'),',true);
+          pro.schreibe(' alte Werte:',true);
           with felders[1-aktuelleFelder] do begin
-            pro:=tProtokollant.create(prot,'iteriereSchritt');
-            pro.schreibe('n<0 bei:',true);
-            pro.schreibe(' t    = ' +floattostr(t)+',',true);
-            pro.schreibe(' i    = ' +inttostr(i)+' (x = '+floattostr(xl+i*dX)+'),',true);
             pro.schreibe(' v_l  = ' +floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]*
                                                 inhalt[i-1].matWerte[j,mfIGamma,false])+',',true);
             pro.schreibe(' n_l  = ' +floattostr(inhalt[i-1].matWerte[j,mfN,false])+',',true);
@@ -996,9 +1101,25 @@ begin
                                                 inhalt[i+1].matWerte[j,mfIGamma,false])+',',true);
             pro.schreibe(' n_r  = ' +floattostr(inhalt[i+1].matWerte[j,mfN,false])+',',true);
             pro.schreibe(' n_r'' = '+floattostr(inhalt[i+1].matWerte[j,mfN,true]),true);
-            pro.destroyall;
-            halt(1);
           end;
+          pro.schreibe(' neue Werte:',true);
+          with felders[aktuelleFelder] do begin
+            pro.schreibe(' v_l  = ' +floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]*
+                                                inhalt[i-1].matWerte[j,mfIGamma,false])+',',true);
+            pro.schreibe(' n_l  = ' +floattostr(inhalt[i-1].matWerte[j,mfN,false])+',',true);
+            pro.schreibe(' n_l'' = '+floattostr(inhalt[i-1].matWerte[j,mfN,true])+',',true);
+            pro.schreibe(' v    = ' +floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]*
+                                                inhalt[i].matWerte[j,mfIGamma,false])+',',true);
+            pro.schreibe(' n    = ' +floattostr(inhalt[i].matWerte[j,mfN,false])+',',true);
+            pro.schreibe(' n''   = '+floattostr(inhalt[i].matWerte[j,mfN,true])+',',true);
+            pro.schreibe(' v_r  = ' +floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]*
+                                                inhalt[i+1].matWerte[j,mfIGamma,false])+',',true);
+            pro.schreibe(' n_r  = ' +floattostr(inhalt[i+1].matWerte[j,mfN,false])+',',true);
+            pro.schreibe(' n_r'' = '+floattostr(inhalt[i+1].matWerte[j,mfN,true]),true);
+          end;
+          pro.free;
+          abbrechen;
+        end;
   {$ENDIF}
   t:=t+dT;
   dT:=max(dT,dTMax/100);
@@ -1025,9 +1146,9 @@ begin
     ns[i]:=ns[i]*dX;
     s:=s+ ' n['+inttostr(i)+']='+floattostr(ns[i]);
     {$IFDEF Dichteueberwachung}
-    if ns[i]>1000 then begin
+    if (not abbruch) and (abs(ns[i]-felders[aktuelleFelder].massen[i])>0.5) then begin
       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);
+      pro.schreibe(' n['+inttostr(i)+'] ('+floattostr(ns[i])+') unterscheidet sich stark vom Anfangswert ('+floattostr(felders[aktuelleFelder].massen[i])+'), es scheinen sehr viele Teilchen entstanden oder verloren gegangen zu sein. Die Simulation wird abgebrochen. (t='+floattostr(t)+')',true);
     end;
     {$ENDIF}
   end;
@@ -1276,7 +1397,7 @@ begin
         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);
+          teilchen[length(teilchen)-1].spezifischeLadung:=-exprToFloat(false,s,kvs,nil);
           kvs.add('qm'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].spezifischeLadung);
           continue;
         end;
@@ -1385,7 +1506,7 @@ begin
   setlength(sighupSimulationen,length(sighupSimulationen)+1);
   sighupSimulationen[length(sighupSimulationen)-1]:=self;
 
-  gitter:=tGitter.create(round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,teilchen,lichter,zeitverfahren,prot,'gitter');
+  gitter:=tGitter.create(self,round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,teilchen,lichter,zeitverfahren,prot,'gitter');
   for i:=0 to length(teilchen)-1 do
     teilchen[i].free;
   setlength(teilchen,0);
diff --git a/Plasmapropagation.lps b/Plasmapropagation.lps
index de2ba33..ba80993 100644
--- a/Plasmapropagation.lps
+++ b/Plasmapropagation.lps
@@ -9,19 +9,20 @@
         <IsPartOfProject Value="True"/>
         <TopLine Value="54"/>
         <CursorPos X="37" Y="34"/>
-        <UsageCount Value="188"/>
+        <UsageCount Value="199"/>
         <Loaded Value="True"/>
       </Unit0>
       <Unit1>
         <Filename Value="Physikunit.pas"/>
         <IsPartOfProject Value="True"/>
+        <IsVisibleTab Value="True"/>
         <EditorIndex Value="1"/>
-        <TopLine Value="943"/>
-        <CursorPos Y="1413"/>
-        <FoldState Value=" T3lf0,5 pjSlM0A[94BjZ084]9Ija09]9ifkK07[I45jO0I013 ppasC0G2 T0=HQ/0F17121b"/>
-        <UsageCount Value="129"/>
+        <TopLine Value="619"/>
+        <CursorPos X="13" Y="648"/>
+        <FoldState Value=" T3mD0,5 pjSlM0A[94BjZ0812114]DLnd0B[944jR08013]Bdl20U2 pp6q10G2 T0\9Q/0F171A"/>
+        <UsageCount Value="140"/>
         <Bookmarks Count="1">
-          <Item0 Y="974"/>
+          <Item0 Y="1078"/>
         </Bookmarks>
         <Loaded Value="True"/>
       </Unit1>
@@ -30,16 +31,15 @@
         <IsPartOfProject Value="True"/>
         <EditorIndex Value="-1"/>
         <CursorPos X="15" Y="46"/>
-        <UsageCount Value="91"/>
+        <UsageCount Value="102"/>
       </Unit2>
       <Unit3>
         <Filename Value="input.plap"/>
         <IsPartOfProject Value="True"/>
-        <IsVisibleTab Value="True"/>
         <EditorIndex Value="3"/>
-        <TopLine Value="6"/>
-        <CursorPos X="19" Y="14"/>
-        <UsageCount Value="90"/>
+        <TopLine Value="31"/>
+        <CursorPos X="26" Y="36"/>
+        <UsageCount Value="101"/>
         <Loaded Value="True"/>
         <DefaultSyntaxHighlighter Value="None"/>
       </Unit3>
@@ -47,32 +47,34 @@
         <Filename Value="linearkombination.inc"/>
         <IsPartOfProject Value="True"/>
         <EditorIndex Value="2"/>
-        <CursorPos X="48" Y="220"/>
-        <UsageCount Value="38"/>
+        <TopLine Value="88"/>
+        <CursorPos X="50" Y="100"/>
+        <UsageCount Value="49"/>
         <Loaded Value="True"/>
       </Unit4>
       <Unit5>
         <Filename Value="input.epost"/>
         <EditorIndex Value="4"/>
-        <CursorPos X="37" Y="22"/>
-        <UsageCount Value="78"/>
+        <TopLine Value="13"/>
+        <CursorPos X="5" Y="22"/>
+        <UsageCount Value="84"/>
         <Loaded Value="True"/>
         <DefaultSyntaxHighlighter Value="None"/>
       </Unit5>
       <Unit6>
         <Filename Value="../units/matheunit.pas"/>
         <EditorIndex Value="-1"/>
-        <TopLine Value="53"/>
-        <CursorPos Y="53"/>
-        <FoldState Value=" T3i905B pj0jV034 piaj60U2-"/>
-        <UsageCount Value="17"/>
+        <TopLine Value="366"/>
+        <CursorPos X="27" Y="390"/>
+        <FoldState Value=" T3i905B pj0jV034 piaj60U511+"/>
+        <UsageCount Value="16"/>
       </Unit6>
       <Unit7>
         <Filename Value="../units/lowlevelunit.pas"/>
         <EditorIndex Value="-1"/>
-        <TopLine Value="4"/>
-        <CursorPos X="86" Y="23"/>
-        <UsageCount Value="7"/>
+        <TopLine Value="46"/>
+        <CursorPos X="3" Y="12"/>
+        <UsageCount Value="9"/>
       </Unit7>
       <Unit8>
         <Filename Value="../units/mystringlistunit.pas"/>
@@ -80,186 +82,186 @@
         <TopLine Value="367"/>
         <CursorPos X="17" Y="390"/>
         <FoldState Value=" T3i3075 piZjD0WQ"/>
-        <UsageCount Value="8"/>
+        <UsageCount Value="7"/>
       </Unit8>
       <Unit9>
         <Filename Value="../epost/werteunit.pas"/>
         <EditorIndex Value="-1"/>
         <TopLine Value="950"/>
         <CursorPos X="30" Y="1054"/>
-        <UsageCount Value="5"/>
+        <UsageCount Value="4"/>
       </Unit9>
       <Unit10>
         <Filename Value="../epost/typenunit.pas"/>
         <EditorIndex Value="-1"/>
         <TopLine Value="347"/>
         <CursorPos X="62" Y="358"/>
-        <UsageCount Value="5"/>
+        <UsageCount Value="4"/>
       </Unit10>
       <Unit11>
         <Filename Value="../units/systemunit.pas"/>
         <EditorIndex Value="-1"/>
         <CursorPos X="3" Y="79"/>
-        <UsageCount Value="5"/>
+        <UsageCount Value="4"/>
       </Unit11>
       <Unit12>
         <Filename Value="/usr/lib/fpc/src/rtl/inc/objpash.inc"/>
         <EditorIndex Value="-1"/>
         <TopLine Value="232"/>
         <CursorPos X="23" Y="192"/>
-        <UsageCount Value="5"/>
+        <UsageCount Value="4"/>
       </Unit12>
       <Unit13>
         <Filename Value="rk14.inc"/>
         <EditorIndex Value="-1"/>
-        <UsageCount Value="8"/>
+        <UsageCount Value="7"/>
       </Unit13>
       <Unit14>
         <Filename Value="rk1210.inc"/>
         <EditorIndex Value="-1"/>
         <TopLine Value="492"/>
         <CursorPos X="28" Y="565"/>
-        <UsageCount Value="9"/>
+        <UsageCount Value="8"/>
       </Unit14>
       <Unit15>
         <Filename Value="rk1412.inc"/>
         <EditorIndex Value="-1"/>
         <TopLine Value="945"/>
         <CursorPos X="5" Y="985"/>
-        <UsageCount Value="9"/>
+        <UsageCount Value="8"/>
       </Unit15>
       <Unit16>
         <Filename Value="rk108.inc"/>
         <EditorIndex Value="-1"/>
         <CursorPos X="28" Y="8"/>
-        <UsageCount Value="9"/>
+        <UsageCount Value="8"/>
       </Unit16>
       <Unit17>
         <Filename Value="rk3_8.inc"/>
         <EditorIndex Value="-1"/>
         <CursorPos X="75" Y="23"/>
-        <UsageCount Value="10"/>
+        <UsageCount Value="9"/>
       </Unit17>
     </Units>
     <JumpHistory Count="30" HistoryIndex="29">
       <Position1>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="341" Column="8" TopLine="320"/>
       </Position1>
       <Position2>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1430" Column="42" TopLine="1419"/>
+        <Caret Line="126" Column="17" TopLine="93"/>
       </Position2>
       <Position3>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="648" Column="104" TopLine="619"/>
+        <Caret Line="521" Column="50" TopLine="433"/>
       </Position3>
       <Position4>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="608" Column="52" TopLine="213"/>
+        <Caret Line="550" Column="51" TopLine="509"/>
       </Position4>
       <Position5>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="977" Column="63" TopLine="953"/>
+        <Caret Line="551" Column="51" TopLine="510"/>
       </Position5>
       <Position6>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="92" Column="38" TopLine="72"/>
+        <Caret Line="665" Column="28" TopLine="326"/>
       </Position6>
       <Position7>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="947" Column="12" TopLine="910"/>
+        <Caret Line="798" Column="44" TopLine="765"/>
       </Position7>
       <Position8>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="930" Column="64" TopLine="707"/>
+        <Caret Line="799" Column="25" TopLine="766"/>
       </Position8>
       <Position9>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="166" Column="72" TopLine="145"/>
+        <Caret Line="800" Column="34" TopLine="767"/>
       </Position9>
       <Position10>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="175" Column="60" TopLine="165"/>
+        <Caret Line="801" Column="51" TopLine="774"/>
       </Position10>
       <Position11>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="946" Column="84" TopLine="929"/>
+        <Caret Line="67" Column="22" TopLine="47"/>
       </Position11>
       <Position12>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="945" Column="20" TopLine="909"/>
+        <Caret Line="425" Column="4" TopLine="210"/>
       </Position12>
       <Position13>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1423" Column="40" TopLine="1394"/>
+        <Caret Line="438" Column="63" TopLine="216"/>
       </Position13>
       <Position14>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="92" Column="38" TopLine="72"/>
+        <Caret Line="801" Column="51" TopLine="768"/>
       </Position14>
       <Position15>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="605" Column="46" TopLine="585"/>
+        <Caret Line="1385" Column="60" TopLine="1232"/>
       </Position15>
       <Position16>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="775" Column="39" TopLine="704"/>
+        <Caret Line="789" Column="63" TopLine="773"/>
       </Position16>
       <Position17>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="797" Column="38" TopLine="778"/>
+        <Caret Line="10" Column="28"/>
       </Position17>
       <Position18>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="131" Column="35" TopLine="111"/>
+        <Caret Line="1024" Column="31" TopLine="784"/>
       </Position18>
       <Position19>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="775" Column="39" TopLine="753"/>
+        <Caret Line="1067" Column="29" TopLine="1056"/>
       </Position19>
       <Position20>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="946" Column="38" TopLine="910"/>
+        <Caret Line="1019" TopLine="861"/>
       </Position20>
       <Position21>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="166" Column="25" TopLine="158"/>
+        <Caret Line="11" Column="3"/>
       </Position21>
       <Position22>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="901" Column="17" TopLine="708"/>
+        <Caret Line="1023" Column="5" TopLine="1000"/>
       </Position22>
       <Position23>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="889" Column="18" TopLine="860"/>
+        <Caret Line="9" Column="33"/>
       </Position23>
       <Position24>
         <Filename Value="Physikunit.pas"/>
+        <Caret Line="991" Column="36" TopLine="787"/>
       </Position24>
       <Position25>
-        <Filename Value="Physikunit.pas"/>
-        <Caret Line="1423" Column="37" TopLine="1394"/>
+        <Filename Value="input.epost"/>
+        <Caret Line="27" TopLine="5"/>
       </Position25>
       <Position26>
-        <Filename Value="Physikunit.pas"/>
-        <Caret TopLine="145"/>
+        <Filename Value="input.epost"/>
+        <Caret Line="44" Column="15" TopLine="5"/>
       </Position26>
       <Position27>
-        <Filename Value="Physikunit.pas"/>
-        <Caret Line="1412" TopLine="774"/>
+        <Filename Value="input.epost"/>
+        <Caret Line="56" Column="27" TopLine="19"/>
       </Position27>
       <Position28>
-        <Filename Value="Physikunit.pas"/>
-        <Caret Line="175" Column="54" TopLine="161"/>
+        <Filename Value="input.epost"/>
+        <Caret Line="47" TopLine="25"/>
       </Position28>
       <Position29>
-        <Filename Value="Physikunit.pas"/>
-        <Caret Line="946" Column="12" TopLine="774"/>
+        <Filename Value="input.epost"/>
+        <Caret Line="57" TopLine="25"/>
       </Position29>
       <Position30>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1424" Column="36" TopLine="988"/>
+        <Caret Line="989" Column="35" TopLine="903"/>
       </Position30>
     </JumpHistory>
   </ProjectSession>
diff --git a/input.epost b/input.epost
index 5a448af..39b7444 100644
--- a/input.epost
+++ b/input.epost
@@ -15,25 +15,47 @@ Palette
   000000
 Ende
 
-Threadanzahl:     10
+Palette
+  Name:            rotblau
+  ff0000
+  ffffff
+  0000ff
+Ende
 
-# !Schleife: $Feld: dDPHIDX dDAYDT dAY dAX
-!Schleife: $Feld: AY DAYDT N1 N2
+Threadanzahl:     7
+
+!Schleife: $symmetrie: symmetrisch asymmetrisch
+
+?$symmetrie  = symmetrisch:  !Schleife: $Feld: DPHIDX DPSIDX1 # AY DAYDT
+?$symmetrie = asymmetrisch: !Schleife: $Feld: N1 # N2
 
 Daten einlesen
   Genauigkeit:     extended
+  xmax:            5
+  tmax:            40
   SpaceTime-Datei: /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/$Feld-*_test.dat
-#  tmin:            0
-#  tmax:            5
-#  xmin:            0.005
-#  xmax:            5
-#  Gamma:           1
 Ende
 
 lineares Bild
-  Palette:        erweiterter Regenbogen
-  maximale und minimale Dichten bestimmen
+?$symmetrie = symmetrisch:   Palette:         rotblau
+?$symmetrie = symmetrisch:   maximale und minimale Dichten bestimmen (symmetrisch)
+?$symmetrie = asymmetrisch:  Palette:         erweiterter Regenbogen
+?$symmetrie = asymmetrisch:  maximale und minimale Dichten bestimmen
   Datei:           /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/$Feld_test.bmp
+  Schriftgröße: 60
+  Achse: oben   10+
+  Achse: links  20+
+  Achse: unten  10+
+  Achse: rechts 20+
+  Rahmen
 Ende
 
+externer Befehl:   convert /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/$Feld_test.bmp /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/$Feld_test.png&
+
+?$symmetrie = symmetrisch:  !Schleifenende
+?$symmetrie = asymmetrisch: !Schleifenende
 !Schleifenende
+
+warte auf externe Befehle
+
+externer Befehl:   rm /home_raid/erich/Dokumente/Prograemmchen/Plasmapropagation/Daten/*_test.bmp
diff --git a/input.plap b/input.plap
index cad3cd7..94a7c54 100644
--- a/input.plap
+++ b/input.plap
@@ -7,10 +7,10 @@ allgemein
   runge-Kutta-3/8
 #  runge-Kutta-4
 #  euler-Vorwärts
-  ortsschritt    10^-3 * λ
-  zeitschritt    10^-3 * T
-  minZeitschritt 10^-10 * T
-  zeit           20 * T
+  ortsschritt    2*10^-3 * λ
+  zeitschritt    2*10^-3 * T
+  minZeitschritt 2*10^-3 * T
+  zeit           40 * T
   diffusionsterm 2                # max. p/Lp
 !setze $breite: (5 * λ)
   breite         $breite
@@ -26,12 +26,14 @@ ausgabenEnde
 !setze $tFwhm:   (2.5 * T)
 !setze $tMitte:  (1 * T)
 
-licht von links      0.1 * ω * 2^(-2*((t-$tMitte)/$tFwhm)^2) * cos(ω*t)    # Zeitableitung des A-Feldes
+# licht von links      0.1 * 2^(-2*((t-$tMitte)/$tFwhm)^2) * ω * cos(ω*t)    # Zeitableitung des A-Feldes
+
+!setze $IonenMassenFaktor: (1836.15267245 + 1838.68366158)
 
 teilchen1
   spezifische Ladung -q/me
   maximaldichte      10
-#  minimaldichte      1
+  minimaldichte      0 # 10^-2
 !setze $profilbreite: (4 * λ)
 !setze $randbreite:   (0.1 * λ)
   verteilung stückweise
@@ -48,10 +50,21 @@ teilchen1
 teilchen1Ende
 
 teilchen2
-  spezifische Ladung (q/me)/2
-  maximaldichte      nmax1*2
-#  minimaldichte      nmin1*2
-  verteilung wie teilchen1
+  spezifische Ladung q/$IonenMassenFaktor/me
+  maximaldichte      nmax1*$IonenMassenFaktor
+  minimaldichte      nmin1*$IonenMassenFaktor
+#  verteilung wie teilchen1
+  verteilung stückweise
+    0
+      (0.1*λ + $breite-$profilbreite)/2 - $randbreite
+    sin((x - (0.1*λ + $breite-$profilbreite)/2 - $randbreite)*π/2/$randbreite)^2
+      (0.1*λ + $breite-$profilbreite)/2
+    1
+      (0.1*λ + $breite+$profilbreite)/2
+    sin((x - (0.1*λ + $breite+$profilbreite)/2 + $randbreite)*π/2/$randbreite)^2
+      (0.1*λ + $breite+$profilbreite)/2 + $randbreite
+    0
+  stückweiseEnde
 teilchen2Ende
 
 Dateiende
diff --git a/linearkombination.inc b/linearkombination.inc
index 0f6ceca..6250006 100644
--- a/linearkombination.inc
+++ b/linearkombination.inc
@@ -55,7 +55,7 @@ begin
       efDAXDT,efDAYDT,efDAZDT,
       efDPhiDX
     ); *)
-  for emF:=efAX to efDPhiDX do // alles außer efA, welchen Ableitung ja nicht berechnet wurde
+  for emF:=erstesEMFMitAbleitung to letztesEMFMitAbleitung do
     emWerte[emF,false]:=
       in1.emWerte[emF,false]
       + fak2  * in2.emWerte[emF,true]  {$IFDEF lkA3}
@@ -97,7 +97,7 @@ begin
       mfGamma,mfIGamma
     ); *)
   for i:=0 to length(matWerte)-1 do // siehe oben
-    for maF:=mfN to mfDPsiDX do
+    for maF:=erstesMatFMitAbleitung to letztesMatFMitAbleitung do
       matWerte[i,maF,false]:=
         in1.matWerte[i,maF,false]
         + fak2  * in2.matWerte[i,maF,true]  {$IFDEF lkA3}
@@ -229,6 +229,8 @@ begin
       {$IFDEF lkA31},fak24,fak25,fak26,fak27,fak28,fak29,fak30,fak31
       {$IFDEF lkA32},fak32
       {$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF}{$ENDIF});
+  for i:=0 to length(inhalt)-1 do
+    inhalt[i].nichtnegativieren;
 end;