1 files changed, 53 insertions, 9 deletions

diff --git a/Physikunit.pas b/Physikunit.pas
index 0ba1418..2a87c63 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -6,9 +6,15 @@ unit Physikunit;
   {$ERROR This program can be compiled only on/for Unix/Linux based systems.}
 {$ENDIF}
 
-{$DEFINE Dichteueberwachung}
-{$DEFINE negativeDichteueberwachung}
-{ $DEFINE exzessiveArrayBereichsTests}
+{$MACRO on}
+
+// nach einer LiKo negative Dichten entfernen:
+{ $DEFINE DichteNichtnegativieren:=perSkalieren}  // insgesamt runter skalieren um negative Dichten aufzufüllen
+{$DEFINE DichteNichtnegativieren:=perEinzelklau} // negative Dichten aus zufälligen Phasenraumpunkten auffüllen
+
+{ $DEFINE negativeDichteueberwachung} // prüfen, ob die Dichten negativ werden
+
+{ $DEFINE exzessiveArrayBereichsTests} // zeitaufwändige(!) Array-Bereichstest für fftw
 
 interface
 
@@ -132,7 +138,9 @@ type
   {$IFDEF exzessiveArrayBereichsTests}
     procedure pruefeArrayEnden(fehler: string);
   {$ENDIF}
+  {$IFDEF dichteNichtnegativieren}
     procedure nichtnegativieren; inline;
+  {$ENDIF}
     procedure berechnePhasenraumAbleitungen; inline;
   public
     emFelder:  array[tEMFeldGroesze,boolean] of pDouble;         // EM-Felder und deren Zeitableitungen
@@ -171,7 +179,9 @@ type
     abbruch:                      boolean;
     dX,iDX,xl,t:                  double;
     kvs:                          tKnownValues;
+  {$IF dichteNichtnegativieren = perEinzelklau}
     impulsRaumPermutationen:      array of tLongintArray;
+  {$ENDIF}
     procedure abbrechen;
   public
 
@@ -182,7 +192,9 @@ type
     destructor destroy; override;
     procedure iteriereSchritt(dT: double);
     procedure dumpErhaltungsgroeszen;
+  {$IF dichteNichtnegativieren = perEinzelklau}
     function gibImpulsRaumPermutation: pTLongintArray;
+  {$ENDIF}
   end;
 
   { tSimulation }
@@ -1349,22 +1361,35 @@ begin
   end;
 end;
 
+{$IFDEF dichteNichtnegativieren}
 procedure tFelder.nichtnegativieren; // Dichten nicht negativ machen
 var
-  i,j:     longint;
-  defizit: double;
-  perm:    pTLongintArray;
+  i,j:             longint;
+  defizit
+  {$IF dichteNichtnegativieren = perSkalieren}
+  ,enfizit
+  {$ENDIF}:        double;
+  {$IF dichteNichtnegativieren = perEinzelklau}
+  perm:            pTLongintArray;
+  {$ENDIF}
 begin
   for i:=0 to length(teilchen)-1 do begin
     defizit:=0;
+    {$IF dichteNichtnegativieren = perSkalieren}
+    enfizit:=0;
+    {$ENDIF}
     for j:=0 to aP*aX-1 do
       if (impulsraum[i,false]+j)^ < 0 then begin
         defizit:=defizit-(impulsraum[i,false]+j)^;
         (impulsraum[i,false]+j)^:=0;
-      end;
+      end
+      {$IF dichteNichtnegativieren = perSkalieren}
+      else enfizit:=enfizit+(impulsraum[i,false]+j)^
+      {$ENDIF};
 
     gesamtDefizit:=gesamtDefizit+defizit;
 
+    {$IF dichteNichtnegativieren = perEinzelklau}
     if defizit>0 then begin
       perm:=gitter.gibImpulsRaumPermutation;
       for j:=0 to aP*aX-1 do begin
@@ -1376,16 +1401,31 @@ begin
         end
         else if (impulsraum[i,false]+perm^[j])^>0 then begin
           defizit:=defizit-(impulsraum[i,false]+perm^[j])^;
-          (impulsraum[i,false]+perm^[j])^:=0;
+         (impulsraum[i,false]+perm^[j])^:=0;
         end;
       end;
     end;
     if defizit>0 then begin
+    {$ENDIF}
+    {$IF dichteNichtnegativieren = perSkalieren}
+    if enfizit=0 then begin
+    {$ENDIF}
       gitter.prot.schreibe('Kann Defizit der Teilchensorte '+inttostr(i+1)+' nicht ausgleichen, '+floattostr(defizit)+' bleibt übrig!',true);
+      gitter.prot.schreibe('Kann Defizit der Teilchensorte '+inttostr(i+1)+' nicht ausgleichen, es ist nichts positives mehr übrig!',true);
       gitter.abbrechen;
+      continue;
     end;
+
+    {$IF dichteNichtnegativieren = perSkalieren}
+    enfizit:=1-defizit/enfizit;
+
+    if defizit>0 then
+      for j:=0 to aP*aX-1 do
+        (impulsraum[i,false]+j)^:=(impulsraum[i,false]+j)^ * enfizit;
+    {$ENDIF}
   end;
 end;
+{$ENDIF}
 
 procedure tFelder.berechnePhasenraumAbleitungen;
 var
@@ -1535,9 +1575,11 @@ begin
   end{of Case};
   xl:=dX/2;
 
-  setlength(impulsRaumPermutationen,10);
+{$IF dichteNichtnegativieren = perEinzelklau}
+  setlength(impulsRaumPermutationen,100);
   for i:=0 to length(impulsRaumPermutationen)-1 do
     impulsRaumPermutationen[i]:=permutation(aX*aP);
+{$ENDIF}
 
   for i:=0 to length(felders)-1 do
     felders[i]:=tFelder.create(aX,aP,deltaX,deltaP,teilchen,lichter,self);
@@ -1633,10 +1675,12 @@ begin
   felders[aktuelleFelder].dumpErhaltungsgroeszen;
 end;
 
+{$IF dichteNichtnegativieren = perEinzelklau}
 function tGitter.gibImpulsRaumPermutation: pTLongintArray;
 begin
   result:=@(impulsRaumPermutationen[random(length(impulsRaumPermutationen))]);
 end;
+{$ENDIF}
 
 // tSimulation *****************************************************************