Der Rand ist jetzt ein regulaerer Punkt.

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

2 files changed, 172 insertions, 39 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;
 
diff --git a/Plasmapropagation.lps b/Plasmapropagation.lps
index 9a7b886..9330883 100644
--- a/Plasmapropagation.lps
+++ b/Plasmapropagation.lps
@@ -15,13 +15,14 @@
       <Unit1>
         <Filename Value="Physikunit.pas"/>
         <IsPartOfProject Value="True"/>
+        <IsVisibleTab Value="True"/>
         <EditorIndex Value="1"/>
-        <TopLine Value="1577"/>
-        <CursorPos Y="714"/>
-        <FoldState Value=" T3mB0{F5 pjdmK0G[945jZ0N4]90jX0A2 piYj90C7]9ajL0v5]9fkW0E[j4Ejg082]93jH03013 T0%LPc071F311#"/>
-        <UsageCount Value="176"/>
+        <TopLine Value="611"/>
+        <CursorPos X="77" Y="819"/>
+        <FoldState Value=" T3mG0{F5 pjdmK0G[945jZ0L3]93jW0C31117 piXnD0F3]9bjH0J4[I4jR0Q6 pkMlL03013 T0)PPc071F]sQ"/>
+        <UsageCount Value="179"/>
         <Bookmarks Count="1">
-          <Item0 Y="1247"/>
+          <Item0 Y="1337"/>
         </Bookmarks>
         <Loaded Value="True"/>
       </Unit1>
@@ -30,16 +31,15 @@
         <IsPartOfProject Value="True"/>
         <EditorIndex Value="-1"/>
         <CursorPos X="20" Y="8"/>
-        <UsageCount Value="138"/>
+        <UsageCount Value="141"/>
       </Unit2>
       <Unit3>
         <Filename Value="input.plap"/>
         <IsPartOfProject Value="True"/>
-        <IsVisibleTab Value="True"/>
         <EditorIndex Value="2"/>
         <TopLine Value="30"/>
         <CursorPos X="23" Y="62"/>
-        <UsageCount Value="137"/>
+        <UsageCount Value="140"/>
         <Loaded Value="True"/>
         <DefaultSyntaxHighlighter Value="None"/>
       </Unit3>
@@ -49,14 +49,14 @@
         <EditorIndex Value="-1"/>
         <TopLine Value="285"/>
         <CursorPos X="20" Y="323"/>
-        <UsageCount Value="85"/>
+        <UsageCount Value="88"/>
       </Unit4>
       <Unit5>
         <Filename Value="linearkombinationen.inc"/>
         <IsPartOfProject Value="True"/>
         <EditorIndex Value="-1"/>
         <TopLine Value="64"/>
-        <UsageCount Value="48"/>
+        <UsageCount Value="51"/>
       </Unit5>
       <Unit6>
         <Filename Value="input.epost"/>
@@ -155,83 +155,126 @@
         <UsageCount Value="10"/>
       </Unit19>
     </Units>
-    <JumpHistory Count="19" HistoryIndex="18">
+    <JumpHistory Count="30" HistoryIndex="29">
       <Position1>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1683" Column="12" TopLine="1135"/>
+        <Caret Line="1229" Column="51" TopLine="1217"/>
       </Position1>
       <Position2>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1685" TopLine="1198"/>
+        <Caret Line="637" Column="74" TopLine="608"/>
       </Position2>
       <Position3>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1686" TopLine="1198"/>
+        <Caret Line="625" Column="14" TopLine="609"/>
       </Position3>
       <Position4>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1687" TopLine="1198"/>
+        <Caret Line="972" Column="8" TopLine="204"/>
       </Position4>
       <Position5>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1223" TopLine="1162"/>
+        <Caret Line="139" Column="3" TopLine="121"/>
       </Position5>
       <Position6>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1227" TopLine="1162"/>
+        <Caret Line="656" Column="16" TopLine="632"/>
       </Position6>
       <Position7>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1217" TopLine="1069"/>
+        <Caret Line="1666" TopLine="973"/>
       </Position7>
       <Position8>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1229" Column="51" TopLine="1217"/>
+        <Caret Line="959" Column="74" TopLine="926"/>
       </Position8>
       <Position9>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="637" Column="74" TopLine="608"/>
+        <Caret Line="956" Column="5" TopLine="936"/>
       </Position9>
       <Position10>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="625" Column="14" TopLine="609"/>
+        <Caret Line="734" Column="22" TopLine="718"/>
       </Position10>
       <Position11>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="972" Column="8" TopLine="204"/>
+        <Caret Line="768" Column="27" TopLine="729"/>
       </Position11>
       <Position12>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="139" Column="3" TopLine="121"/>
+        <Caret Line="15"/>
       </Position12>
       <Position13>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="656" Column="16" TopLine="632"/>
+        <Caret Line="731" Column="24" TopLine="703"/>
       </Position13>
       <Position14>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="1666" TopLine="973"/>
+        <Caret Line="730" Column="28" TopLine="703"/>
       </Position14>
       <Position15>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="622" TopLine="713"/>
+        <Caret Line="765" Column="66" TopLine="729"/>
       </Position15>
       <Position16>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="30" Column="61" TopLine="46"/>
+        <Caret Line="632" Column="42" TopLine="603"/>
       </Position16>
       <Position17>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="959" Column="74" TopLine="926"/>
+        <Caret Line="644" Column="26" TopLine="612"/>
       </Position17>
       <Position18>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="956" Column="5" TopLine="936"/>
+        <Caret Line="112" Column="60" TopLine="82"/>
       </Position18>
       <Position19>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="30" Column="64" TopLine="75"/>
+        <Caret Line="630" Column="30" TopLine="609"/>
       </Position19>
+      <Position20>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="610" Column="26" TopLine="587"/>
+      </Position20>
+      <Position21>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="101" Column="71" TopLine="81"/>
+      </Position21>
+      <Position22>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="839" Column="60" TopLine="802"/>
+      </Position22>
+      <Position23>
+        <Filename Value="Physikunit.pas"/>
+      </Position23>
+      <Position24>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="610" Column="49" TopLine="587"/>
+      </Position24>
+      <Position25>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="101" Column="61" TopLine="75"/>
+      </Position25>
+      <Position26>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="628" Column="21" TopLine="523"/>
+      </Position26>
+      <Position27>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="114" TopLine="88"/>
+      </Position27>
+      <Position28>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="643" TopLine="588"/>
+      </Position28>
+      <Position29>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="638" Column="33" TopLine="619"/>
+      </Position29>
+      <Position30>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="640" Column="27" TopLine="620"/>
+      </Position30>
     </JumpHistory>
   </ProjectSession>
 </CONFIG>