Der Rand ist jetzt ein regulaerer Punkt.

Gradient in Funktion gekapselt, ist jetzt O(6).

1 files changed, 99 insertions, 9 deletions

diff --git a/Physikunit.pas b/Physikunit.pas
index 5d8881a..27cca80 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -91,12 +91,14 @@ type
 
   tImpulsPunkt = class(tObject) // Besetzungsdichten verschiedener Teilchenspezies für ein Paar (r;p)
   private
-    raumpunkt: tRaumPunkt;
-    nachbarn:  array[0..1,boolean] of tImpulsPunkt; // x- und px-Nachbarn (jeweil in Richtung - und +)
+    raumpunkt:  tRaumPunkt;
+    nachbarn:   array[0..1,boolean] of tImpulsPunkt; // x- und px-Nachbarn (jeweil in Richtung - und +)
+    istNullPkt: boolean;
   public
     werte:   array of array[boolean] of extended; // Besetzungsdichten und deren Ableitungen für verschiedene Teilchenspezies
     iMGamma: array of extended;                   // 1/m/gamma (= v/p)
-    constructor create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint);
+    constructor create(anzTeilchen: longint); overload; // NullPunkt erstellen
+    constructor create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint); overload;
     destructor destroy; override;
     {$DEFINE LiKotImpulsPunktHeader}
     {$INCLUDE linearkombinationen.inc}
@@ -106,6 +108,9 @@ type
     function gibMirWas(var defizit: extended; entfernung, teilchen,richtung: longint; positiv: boolean): boolean; inline; // ich habe jemanden erreicht
     procedure setzeNull; inline;
     procedure berechneAbleitungen(pX,iDX,iDPX: extended); inline;
+    function gradient(rtg,tlc: longint; iD,fak: extended): extended; inline;
+    property istNullPunkt: boolean read istNullPkt;
+    function gibNullPunkt: tImpulsPunkt;
   end;
 
   { tRaumPunkt }
@@ -580,12 +585,37 @@ end;
 
 // tImpulsPunkt ****************************************************************
 
-constructor tImpulsPunkt.create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint);
+constructor tImpulsPunkt.create(anzTeilchen: longint);
 var
   i: longint;
   b: boolean;
 begin
   inherited create;
+  istNullPkt:=true; // als NullPunkt erstellen
+  raumPunkt:=nil;
+  fillchar(werte,sizeof(werte),#0);
+  setlength(werte,anzTeilchen);
+  for i:=0 to length(werte)-1 do
+    for b:=false to true do
+      werte[i,b]:=0;
+  fillchar(iMGamma,sizeof(iMGamma),#0);
+  setlength(iMGamma,length(werte));
+  for i:=0 to length(iMGamma)-1 do
+    iMGamma[i]:=1;
+  nachbarn[0,false]:=self;
+  nachbarn[0,true]:=self;
+  nachbarn[1,false]:=self;
+  nachbarn[1,true]:=self;
+end;
+
+constructor tImpulsPunkt.create(rp: tRaumPunkt; linksRaum,linksImpuls: tImpulsPunkt; anzTeilchen: longint);
+var
+  i:  longint;
+  b:  boolean;
+  np: tImpulsPunkt;
+begin
+  inherited create;
+  istNullPkt:=false;
   raumPunkt:=rp;
   fillchar(werte,sizeof(werte),#0);
   setlength(werte,anzTeilchen);
@@ -601,8 +631,21 @@ begin
   nachbarn[1,false]:=linksImpuls;
   nachbarn[1,true]:=nil;
   for i:=0 to 1 do
-    if assigned(nachbarn[i,false]) then
+    if assigned(nachbarn[i,false]) and
+       not nachbarn[i,false].istNullPunkt then
       nachbarn[i,false].nachbarn[i,true]:=self;
+  np:=nil;
+  for i:=0 to 1 do
+    if assigned(nachbarn[i,false]) and
+       not assigned(np) then
+      np:=nachbarn[i,false].gibNullPunkt;
+  if not assigned(np) then
+    np:=tImpulsPunkt.create(anzTeilchen);
+  for i:=0 to 1 do
+    for b:=false to true do
+      if not assigned(nachbarn[i,b]) then begin
+        nachbarn[i,b]:=np;
+      end;
 end;
 
 destructor tImpulsPunkt.destroy;
@@ -612,7 +655,8 @@ var
 begin
   for i:=0 to 1 do
     for b:=false to true do
-      if assigned(nachbarn[i,b]) then
+      if assigned(nachbarn[i,b]) and
+         not nachbarn[i,b].istNullPunkt then
         nachbarn[i,b].nachbarn[i,not b]:=nil;
   setlength(werte,0);
   setlength(iMGamma,0);
@@ -726,7 +770,9 @@ begin
         (  raumpunkt.emFelder[efEX,false]
          + raumpunkt.matFelder[i,mfPY,false] * iMGamma[i] * raumpunkt.emFelder[efBZ,false]);
     werte[i,true]:=
-      tX * (
+      gradient(0,i,iDX,tX) +
+      gradient(1,i,iDPX,tP);
+(*      tX * (
           (nachbarn[0,true].werte[i,false] - nachbarn[0,false].werte[i,false])
          + sign(tX)*raumpunkt.felder.teilchen[i].diffusion[0] * (
              nachbarn[0,true].werte[i,false] +
@@ -741,7 +787,7 @@ begin
              nachbarn[1,false].werte[i,false] -
              2*werte[i,false]
            )
-        )* iDPX/2;
+        )* iDPX/2; *)
   end;
 
     // Der jeweils zweite Summand entspricht (in etwa) einer thermischen
@@ -749,6 +795,50 @@ begin
     // d(ln n)/d[xp] * d[XP] = "diffusion" gewährleisten.
 end;
 
+function tImpulsPunkt.gradient(rtg,tlc: longint; iD,fak: extended): extended;
+begin
+
+  {
+  // Variante mit Diffusionsterm
+  result:=
+    fak * ((nachbarn[rtg,true].werte[tlc,false] - nachbarn[rtg,false].werte[tlc,false])
+           + sign(fak)*raumpunkt.felder.teilchen[tlc].diffusion[rtg] * (
+               nachbarn[rtg,true].werte[tlc,false] +
+               nachbarn[rtg,false].werte[tlc,false] -
+               2*werte[tlc,false]
+             )
+          ) * iD/2;
+  }
+
+  result:=
+    fak * ( - 1/60 * nachbarn[rtg,false].nachbarn[rtg,false].nachbarn[rtg,false].werte[tlc,false]
+            + 3/20 * nachbarn[rtg,false].nachbarn[rtg,false].werte[tlc,false]
+            - 3/4  * nachbarn[rtg,false].werte[tlc,false]
+            + 3/4  * nachbarn[rtg,true].werte[tlc,false]
+            - 3/20 * nachbarn[rtg,true].nachbarn[rtg,true].werte[tlc,false]
+            + 1/60 * nachbarn[rtg,true].nachbarn[rtg,true].nachbarn[rtg,true].werte[tlc,false]
+          ) * iD;
+
+
+//  probier mal den Gradienten anders zu berechnen und zwar so:
+//    https://en.wikipedia.org/wiki/Finite_difference_coefficient
+end;
+
+function tImpulsPunkt.gibNullPunkt: tImpulsPunkt;
+var
+  i: longint;
+begin
+  if istNullPunkt then begin
+    result:=self;
+    exit;
+  end;
+  result:=nil;
+  for i:=0 to 1 do
+    if assigned(nachbarn[i,false]) and
+       not assigned(result) then
+      result:=nachbarn[i,false].gibNullPunkt;
+end;
+
 // tRaumPunkt ******************************************************************
 
 constructor tRaumPunkt.create(linkerNachbar: tRaumPunkt; derBesitzer: tFelder; teilchen,_aP: longint; _dP: extended);
@@ -871,7 +961,7 @@ var
   i: longint;
 begin
   // df/dt = - v Nabla f - q(E + v/c x B) Nabla_p f
-  for i:=1 to aP-2 do // Rand wird ausgelassen!
+  for i:=0 to aP-1 do
     phasenraum[i].berechneAbleitungen((i-(aP div 2))*dP,iDX,1/dP);
 end;