nichtnegativieren in Physikunit.pas beschleunigt

1 files changed, 28 insertions, 40 deletions

diff --git a/Physikunit.pas b/Physikunit.pas
index 935dcad..4dbc2da 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -171,6 +171,7 @@ type
     abbruch:                      boolean;
     dX,iDX,xl,t:                  double;
     kvs:                          tKnownValues;
+    impulsRaumPermutation:        tLongintArray;
     procedure abbrechen;
   public
 
@@ -1348,52 +1349,37 @@ end;
 
 procedure tFelder.nichtnegativieren; // Dichten nicht negativ machen
 var
-  i,j,k,l,m,koo: longint;
-  defizit:       double;
-  schr,zSchr:    array[0..3] of longint;
+  i,j:     longint;
+  defizit: double;
 begin
-  schr[0]:=-aP;
-  schr[1]:=-1;
-  schr[2]:=aP;
-  schr[3]:=1;
-  for i:=0 to 3 do
-    zSchr[i]:=schr[(i+1) mod 4] + schr[(i+2) mod 4];
-  for i:=0 to length(impulsraum)-1 do
-    for j:=0 to aX*aP-1 do
+  for i:=0 to length(teilchen)-1 do begin
+    defizit:=0;
+    for j:=0 to aP*aX-1 do
       if (impulsraum[i,false]+j)^ < 0 then begin
-        defizit:=-(impulsraum[i,false]+j)^;
-        gesamtDefizit:=gesamtDefizit + defizit;
+        defizit:=defizit-(impulsraum[i,false]+j)^;
         (impulsraum[i,false]+j)^:=0;
+      end;
+
+    gesamtDefizit:=gesamtDefizit+defizit;
 
-        for k:=1 to aX+aP do
-          if defizit>0 then
-            for l:=0 to 3 do
-              if defizit>0 then
-                for m:=0 to k-1 do begin
-                  koo:=j + schr[l]*k + zSchr[l]*m;
-                  if (koo<0) or (koo>=aX*aP) then continue;
-
-                  if (impulsraum[i,false]+koo)^>defizit then begin
-                    (impulsraum[i,false]+koo)^:=(impulsraum[i,false]+koo)^-defizit;
-                    defizit:=0;
-                    break;
-                  end
-                  else if (impulsraum[i,false]+koo)^>=0 then begin
-                    defizit:=defizit-(impulsraum[i,false]+koo)^;
-                    (impulsraum[i,false]+koo)^:=0;
-                  end
-                  else begin
-                    defizit:=defizit-(impulsraum[i,false]+koo)^;
-                    gesamtDefizit:=gesamtDefizit-(impulsraum[i,false]+koo)^;
-                    (impulsraum[i,false]+koo)^:=0;
-                  end;
-                end;
-
-        if defizit>0 then begin
-          gitter.prot.schreibe('Kann Defizit der Teilchensorte '+inttostr(i+1)+' nicht ausgleichen, '+floattostr(defizit)+' bleibt übrig!',true);
-          gitter.abbrechen;
+    if defizit>0 then
+      for j:=0 to aP*aX-1 do begin
+        if (impulsraum[i,false]+gitter.impulsRaumPermutation[j])^>=defizit then begin
+          (impulsraum[i,false]+gitter.impulsRaumPermutation[j])^:=
+            (impulsraum[i,false]+gitter.impulsRaumPermutation[j])^-defizit;
+          defizit:=0;
+          break;
+        end
+        else if (impulsraum[i,false]+gitter.impulsRaumPermutation[j])^>0 then begin
+          defizit:=defizit-(impulsraum[i,false]+gitter.impulsRaumPermutation[j])^;
+          (impulsraum[i,false]+gitter.impulsRaumPermutation[j])^:=0;
         end;
       end;
+    if defizit>0 then begin
+      gitter.prot.schreibe('Kann Defizit der Teilchensorte '+inttostr(i+1)+' nicht ausgleichen, '+floattostr(defizit)+' bleibt übrig!',true);
+      gitter.abbrechen;
+    end;
+  end;
 end;
 
 procedure tFelder.berechnePhasenraumAbleitungen;
@@ -1544,6 +1530,8 @@ begin
   end{of Case};
   xl:=dX/2;
 
+  impulsRaumPermutation:=permutation(aX*aP);
+
   for i:=0 to length(felders)-1 do
     felders[i]:=tFelder.create(aX,aP,deltaX,deltaP,teilchen,lichter,self);