nichtnegativieren eingebaut

1 files changed, 34 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);