nichtnegativieren eingebaut

3 files changed, 42 insertions, 11 deletions

diff --git a/Physikunit.pas b/Physikunit.pas
index 7dddb8c..4b73d8b 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -109,6 +109,7 @@ type
     fftTmp,lFftTmp:     Pcomplex_double;                         // Zwischenspeicher für ffts
     procedure initialisiereDichte(ort,tlc: longint; breite,n: double); inline;
     procedure pruefeArrayEnden(fehler: string);
+    procedure nichtnegativieren;
   public
     emFelder:  array[tEMFeldGroesze,boolean] of pDouble;         // EM-Felder und deren Zeitableitungen
                                                                    // dPhiDX und B[xyz] haben keine sinnvolle Ableitung hier!
@@ -547,17 +548,6 @@ begin
   end;
 end;
 
-(*
-function tRaumPunkt.nichtnegativieren: double; // Dichten nicht negativ machen
-var
-  i: longint;
-begin
-  result:=0;
-  for i:=0 to length(phasenraum)-1 do
-    result:=result+phasenraum[i].nichtnegativieren;
-end;
-
-*)
 // tFelder *********************************************************************
 
 constructor tFelder.create(anzX,anzP: longint; deltaX,deltaP: double; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
@@ -1066,6 +1056,39 @@ begin
   closefile(f);
 end;
 
+procedure tFelder.nichtnegativieren; // Dichten nicht negativ machen
+var
+  i,j,k,l: longint;
+  defizit: double;
+begin
+  for i:=0 to length(impulsraum)-1 do
+    for j:=0 to aX*aP-1 do
+      if (impulsraum[i,false]+j)^ < 0 then begin
+        defizit:=-(impulsraum[i,false]+j)^;
+        gesamtDefizit:=gesamtDefizit + defizit;
+        (impulsraum[i,false]+j)^:=0;
+        k:=j;
+        l:=1;
+        while (defizit>0) and (l<2*aX*aP) do begin
+          k:=k+l*(1-2*byte(odd(l)));
+          inc(l);
+          if (k<0) or (k>=aX*aP) then continue;
+          if (impulsraum[i,false]+k)^>defizit then begin
+            (impulsraum[i,false]+k)^:=(impulsraum[i,false]+k)^-defizit;
+            defizit:=0;
+          end
+          else if (impulsraum[i,false]+k)^>0 then begin
+            defizit:=defizit-(impulsraum[i,false]+k)^;
+            (impulsraum[i,false]+k)^:=0;
+          end;
+        end;
+        if defizit>0 then begin
+          gitter.prot.schreibe('Kann Defizit der Teilchensorte '+inttostr(i+1)+' nicht ausgleichen, '+floattostr(defizit)+' bleibt übrig!');
+          gitter.abbrechen;
+        end;
+      end;
+end;
+
 // tGitter *********************************************************************
 
 constructor tGitter.create(derBesitzer: tSimulation; aX,aP: longint; deltaX,deltaP: double; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; name: string);
diff --git a/genLikoInc b/genLikoInc
index ab0035d..e4e49d6 100755
--- a/genLikoInc
+++ b/genLikoInc
@@ -101,5 +101,9 @@ do
 done
 echo ";" >> ${file}
 echo >> ${file}
+echo "  {\$IFDEF Dichteueberwachung}" >> ${file}
+echo "  nichtnegativieren;" >> ${file}
+echo "  {\$ENDIF}" >> ${file}
+echo >> ${file}
 echo "end;" >> ${file}
 
diff --git a/linearkombination.inc b/linearkombination.inc
index 05aeeee..bca3601 100644
--- a/linearkombination.inc
+++ b/linearkombination.inc
@@ -259,4 +259,8 @@ begin
         + fak32*(in32.impulsraum[i,true]+j)^
         {$ENDIF};
 
+  {$IFDEF Dichteueberwachung}
+  nichtnegativieren;
+  {$ENDIF}
+
 end;