lauffaehig, es entstehen aber noch "Ripple" im Phasenraum ...

1 files changed, 145 insertions, 83 deletions

diff --git a/Physikunit.pas b/Physikunit.pas
index 4861a14..5d8881a 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -27,6 +27,7 @@ type
     eqRho,eqJX,eqJY
   );
   tMaterieFeldGroesze = (mfPY);
+  tMaterieSpeicherGroesze = (msPX,msPXSqr,msPY,msVX,msVY,msN,msPXRipple);
   tTeilchenSpeziesGroeszen = (tsgMasse,tsgIQdrMasse,tsgLadung);
   tEMQuellen = array[tEMQuellGroesze] of extended;
 
@@ -46,9 +47,11 @@ type
   tAusgabeDatei = class
   private
     ableitung:          boolean;
-    emFelder:           tEMFeldGroesze;
-    emQuellen:          tEMQuellGroesze;
-    teilchen,           // -2: EM-Feld; -1: Gesamt-EM-Quelle; sonst: EM-Quelle der entsprechenden Teilchen
+    emFeld:             tEMFeldGroesze;
+    emQuelle:           tEMQuellGroesze;
+    matFeld:            tMaterieSpeicherGroesze;
+    wasSpeichern,       // 0: EM-Feld; 1: EM-Quelle; 2: Mat-Feld
+    teilchen,           // -1: Gesamt-Feld; sonst: Feld des entsprechenden Teilchen
       nNum,tCnt:        longint;
     pre,suf:            string;
     datei:              file;
@@ -132,7 +135,7 @@ type
     {$UNDEF LiKotRaumPunktHeader}
     function nichtnegativieren: extended; inline;
     function impulsIntegral(teilchen: longint; emQ: tEMQuellGroesze): extended; overload; inline;
-    function impulsIntegral(teilchen: longint): extended; overload; inline;
+    function impulsIntegral(teilchen: longint; maF: tMaterieSpeicherGroesze): extended; overload; inline;
     procedure initialisiereDichte(teilchen: longint; breite,n: extended); inline;
     procedure reflektiereTeilchen(rechts: boolean); inline;
     procedure setzeNull; inline;
@@ -211,6 +214,9 @@ const
   emQuellNamen: array[tEMQuellGroesze] of string = (
     'RHO','JX','JY'
   );
+  matSpeicherNamen: array[tMaterieSpeicherGroesze] of string = (
+    'PX','PXSQR','PY','VX','VY','N','PXRIPPLE'
+  );
   minAusgabeBuffer = 1024*1024;
 
 var
@@ -222,6 +228,7 @@ constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokol
 var
   emF:     tEMFeldGroesze;
   emQ:     tEMQuellGroesze;
+  maF:     tMaterieSpeicherGroesze;
   num:     longint;
   abl,gef: boolean;
 begin
@@ -238,8 +245,10 @@ begin
   else
     teilchen:=-1;
 
-  emFelder:=efAX;
-  emQuellen:=eqRho;
+  emFeld:=efAX;
+  emQuelle:=eqRho;
+  matFeld:=msN;
+  wasSpeichern:=0;
   feldName:=ansiUpperCase(feldName);
   nNum:=0; // Header 0 wurde also noch nicht geschrieben
   zeitAnz:=-1;
@@ -250,7 +259,8 @@ begin
     if not gef then
       for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do
         if copy('D',1,byte(abl))+emFeldNamen[emF] = feldName then begin
-          emFelder:=emF;
+          emFeld:=emF;
+          wasSpeichern:=0;
           ableitung:=abl;
           teilchen:=-2;
           gef:=true;
@@ -261,24 +271,45 @@ begin
     for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do
       if emQuellNamen[emQ] = feldName then begin
         ableitung:=false;
-        emQuellen:=emQ;
+        emQuelle:=emQ;
+        wasSpeichern:=1;
+        gef:=true;
+        break;
+      end;
+  if not gef then
+    for maF:=low(tMaterieSpeicherGroesze) to high(tMaterieSpeicherGroesze) do
+      if matSpeicherNamen[maF] = feldName then begin
+        ableitung:=false;
+        matFeld:=maF;
+        wasSpeichern:=2;
         gef:=true;
         break;
       end;
 
   if not gef then begin
     pro.schreibe('tAusgabeDatei.create kennt Größe '''+feldName+''' nicht!',true);
+    nam:='';
+    for abl:=false to true do
+      for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do
+        nam:=nam+', '+copy('D',1,byte(abl))+emFeldNamen[emF];
+    for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do
+      nam:=nam+', '+emQuellNamen[emQ];
+    for maF:=low(tMaterieSpeicherGroesze) to high(tMaterieSpeicherGroesze) do
+      nam:=nam+', '+matSpeicherNamen[maF];
+    delete(nam,1,2);
+    pro.schreibe('Ich kenne nur: '+nam,true);
+
     pro.destroyAll;
     halt(1);
   end;
 
-  if teilchen=-2 then
-    pre:=emFeldNamen[emF]
-  else begin
-    pre:=emQuellNamen[emQ];
-    if teilchen>=0 then
-      pre:=pre+inttostr(teilchen+1);
+  case wasSpeichern of
+    0: pre:=emFeldNamen[emFeld];
+    1: pre:=emQuellNamen[emQuelle];
+    2: pre:=matSpeicherNamen[matFeld];
   end;
+  if teilchen>=0 then
+    pre:=pre+inttostr(teilchen+1);
 
   if ableitung then
     pre:='d'+pre;
@@ -320,13 +351,13 @@ end;
 
 function tAusgabeDatei.dump: string;
 begin
-  if teilchen=-2 then
-    result:=emFeldNamen[emFelder]
-  else begin
-    result:=emQuellNamen[emQuellen];
-    if teilchen>=0 then
-      result:=result+inttostr(teilchen);
+  case wasSpeichern of
+    0: result:=emFeldNamen[emFeld];
+    1: result:=emQuellNamen[emQuelle];
+    2: result:=matSpeicherNamen[matFeld];
   end;
+  if teilchen>=0 then
+    result:=result+'['+inttostr(teilchen+1)+']';
   if ableitung then
     result:=result+'''';
   result:=result+' -> '+pre+'-$i'+suf;
@@ -385,40 +416,40 @@ begin
   inc(bufPos,sizeof(integer));
 
   sX:=cX-Min(gitter.dX,sDX)/2;
-  if teilchen=-2 then begin
-    for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
-      while cX>=sX do begin
-        dec(cnt);
-        move(gitter.felders[gitter.aktuelleFelder].inhalt[i].emFelder[emFelder,ableitung],(buf+bufPos)^,sizeof(extended));
-        inc(bufPos,sizeof(extended));
-        sX:=sX+sDX;
+  case wasSpeichern of
+    0: // em-Feld speichern
+      for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
+        while cX>=sX do begin
+          dec(cnt);
+          move(gitter.felders[gitter.aktuelleFelder].inhalt[i].emFelder[emFeld,ableitung],(buf+bufPos)^,sizeof(extended));
+          inc(bufPos,sizeof(extended));
+          sX:=sX+sDX;
+        end;
+        cX:=cX+gitter.dX;
       end;
-      cX:=cX+gitter.dX;
-    end;
-  end
-  else if teilchen=-1 then begin
-    for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
-      while cX>=sX do begin
-        dec(cnt);
-        move(gitter.felders[gitter.aktuelleFelder].inhalt[i].emQuellen[emQuellen],(buf+bufPos)^,sizeof(extended));
-        inc(bufPos,sizeof(extended));
-        sX:=sX+sDX;
+    1:  // em-Quelle
+      for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
+        while cX>=sX do begin
+          dec(cnt);
+          val:=gitter.felders[gitter.aktuelleFelder].inhalt[i].impulsIntegral(teilchen,emQuelle);
+          move(val,(buf+bufPos)^,sizeof(extended));
+          inc(bufPos,sizeof(extended));
+          sX:=sX+sDX;
+        end;
+        cX:=cX+gitter.dX;
       end;
-      cX:=cX+gitter.dX;
-    end;
-  end
-  else begin
-    for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
-      while cX>=sX do begin
-        dec(cnt);
-        val:=gitter.felders[gitter.aktuelleFelder].inhalt[i].impulsIntegral(teilchen,emQuellen);
-        move(val,(buf+bufPos)^,sizeof(extended));
-        inc(bufPos,sizeof(extended));
-        sX:=sX+sDX;
+    2:  // Materiefeld
+      for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
+        while cX>=sX do begin
+          dec(cnt);
+          val:=gitter.felders[gitter.aktuelleFelder].inhalt[i].impulsIntegral(teilchen,matFeld);
+          move(val,(buf+bufPos)^,sizeof(extended));
+          inc(bufPos,sizeof(extended));
+          sX:=sX+sDX;
+        end;
+        cX:=cX+gitter.dX;
       end;
-      cX:=cX+gitter.dX;
-    end;
-  end;
+  end{of case};
   if cnt<>0 then begin
     pro.schreibe('Falsche Anzahl an Ortsschritten geschrieben ('+inttostr(cnt)+')!',true);
     pro.destroyall;
@@ -590,17 +621,18 @@ end;
 
 function tImpulsPunkt.nichtnegativieren: extended; // Dichten nicht negativ machen
 var
-  i,richtung,entfernung:     longint;
   defizit:                   extended;
+  i:                         longint;
   pro:                       tProtokollant;
-  jemandErreicht,vorZurueck: boolean;
+(*  i,richtung,entfernung:     longint;
+  jemandErreicht,vorZurueck: boolean; *)
 begin
   result:=0;
-  for i:=0 to length(werte)-1 do
+  for i:=0 to length(werte)-1 do begin
     if werte[i,false]<0 then begin
       defizit:=-werte[i,false];
       result:=result+defizit;
-      werte[i,false]:=0;
+      werte[i,false]:=0; (*
       entfernung:=1;
       jemandErreicht:=true;
       while (defizit>0) and jemandErreicht do begin
@@ -616,8 +648,15 @@ begin
         pro:=tProtokollant.create(raumpunkt.felder.gitter.prot,'WertePunkt.nichtnegativieren');
         pro.schreibe('Ich konnte ein Teilchendefizit der Sorte '+inttostr(i+1)+' nicht auffüllen!',true);
         raumpunkt.felder.gitter.abbrechen;
-      end;
+      end;  *)
     end;
+    if werte[i,false]>1E10 then begin
+      pro:=tProtokollant.create(raumpunkt.felder.gitter.prot,'impulsPunkt.nichtnegativieren');
+      pro.schreibe('An einer Stelle im Phasenraum sind bereits mehr als 10^100 Teilchen, ich breche ab!',true);
+      pro.free;
+      raumpunkt.felder.gitter.abbrechen;
+    end;
+  end;
 end;
 
 procedure tImpulsPunkt.akkumuliereEMQuellen(var emQuellen: tEMQuellen; pX,dVP: extended);
@@ -688,16 +727,16 @@ begin
          + raumpunkt.matFelder[i,mfPY,false] * iMGamma[i] * raumpunkt.emFelder[efBZ,false]);
     werte[i,true]:=
       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] - nachbarn[0,false].werte[i,false])
+         + sign(tX)*raumpunkt.felder.teilchen[i].diffusion[0] * (
              nachbarn[0,true].werte[i,false] +
              nachbarn[0,false].werte[i,false] -
              2*werte[i,false]
            )
         ) * iDX/2
     + tP * (
-        (  nachbarn[1,true].werte[i,false] - nachbarn[1,false].werte[i,false])
-         - sign(tP)*raumpunkt.felder.teilchen[i].diffusion[1] * (
+          (nachbarn[1,true].werte[i,false] - nachbarn[1,false].werte[i,false])
+         + sign(tP)*raumpunkt.felder.teilchen[i].diffusion[1] * (
              nachbarn[1,true].werte[i,false] +
              nachbarn[1,false].werte[i,false] -
              2*werte[i,false]
@@ -846,8 +885,6 @@ begin
 end;
 
 function tRaumPunkt.impulsIntegral(teilchen: longint; emQ: tEMQuellGroesze): extended;
-var
-  i: longint;
 begin
   if teilchen<0 then begin
     result:=emQuellen[emQ]; // das ist leicht :-)
@@ -857,12 +894,50 @@ begin
   result:=0;
 
   case emQ of
-    eqRho:
+    eqRho: result:=impulsIntegral(teilchen,msN);
+    eqJX:  result:=impulsIntegral(teilchen,msVX);
+    eqJY:  result:=impulsIntegral(teilchen,msVY);
+  end{of case};
+  result:=result * felder.teilchen[teilchen].eigenschaften[tsgLadung];
+end;
+
+function tRaumPunkt.impulsIntegral(teilchen: longint; maF: tMaterieSpeicherGroesze): extended;
+var
+  i: longint;
+begin
+  if teilchen<0 then begin
+    result:=0;
+    for i:=0 to length(felder.teilchen)-1 do
+      result:=result+impulsIntegral(i,maF);
+    exit;
+  end;
+
+  result:=0;
+
+  case maF of
+    msN:
       for i:=0 to aP-1 do
         result:=
           result +
           phasenraum[i].werte[teilchen,false]*dP;
-    eqJX: begin
+    msPX:
+      for i:=0 to aP-1 do
+        result:=
+          result +
+          phasenraum[i].werte[teilchen,false]*dP * (i-aP/2)*dP;
+    msPXSqr:
+      for i:=0 to aP-1 do
+        result:=
+          result +
+          phasenraum[i].werte[teilchen,false]*dP * sqr((i-aP/2)*dP);
+    msPY: begin
+      for i:=0 to aP-1 do
+        result:=
+          result +
+          phasenraum[i].werte[teilchen,false]*dP;
+      result:=result * matFelder[teilchen,mfPY,false];
+    end;
+    msVX: begin
       for i:=0 to aP-1 do
         result:=
           result +
@@ -870,7 +945,7 @@ begin
           / sqrt(1 + (sqr((i-aP/2)*dP) + sqr(matFelder[teilchen,mfPY,false]))*felder.teilchen[teilchen].eigenschaften[tsgIQdrMasse]);
       result:=result / felder.teilchen[teilchen].eigenschaften[tsgMasse];
     end;
-    eqJY: begin
+    msVY: begin
       for i:=0 to aP-1 do
         result:=
           result +
@@ -878,24 +953,11 @@ begin
           / sqrt(1 + (sqr((i-aP/2)*dP) + sqr(matFelder[teilchen,mfPY,false]))*felder.teilchen[teilchen].eigenschaften[tsgIQdrMasse]);
       result:=result * matFelder[teilchen,mfPY,false] / felder.teilchen[teilchen].eigenschaften[tsgMasse];
     end;
+    msPXRipple:
+      for i:=1 to aP-1 do
+        result:=
+          result + abs(phasenraum[i].werte[teilchen,false]-phasenraum[i-1].werte[teilchen,false])*dP;
   end{of case};
-  result:=result * felder.teilchen[teilchen].eigenschaften[tsgLadung];
-end;
-
-function tRaumPunkt.impulsIntegral(teilchen: longint): extended;
-var
-  i,j: longint;
-begin
-  result:=0;
-
-  if teilchen<0 then begin
-    for j:=0 to length(matFelder)-1 do
-      for i:=0 to aP-1 do
-        result:=result + phasenraum[i].werte[j,false]*dP;
-  end
-  else
-    for i:=0 to aP-1 do
-      result:=result + phasenraum[i].werte[teilchen,false]*dP;
 end;
 
 procedure tRaumPunkt.initialisiereDichte(teilchen: longint; breite,n: extended);
@@ -1086,7 +1148,7 @@ begin
   for i:=0 to length(massen)-1 do begin
     dens:=0;
     for j:=0 to length(inhalt)-1 do
-      dens:=dens+inhalt[j].impulsIntegral(i);
+      dens:=dens+inhalt[j].impulsIntegral(i,msN);
     dens:=dens*teilchen[i].eigenschaften[tsgMasse]*gitter.dX;
     prot.schreibe('n['+inttostr(i+1)+'] = '+floattostr(dens)+' (relative Abweichung: '+floattostr(dens/massen[i]-1)+')',true);
   end;