Phasenraumableitungen werden mit gedumpt

1 files changed, 59 insertions, 34 deletions

diff --git a/Physikunit.pas b/Physikunit.pas
index 595401a..1327173 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -111,6 +111,7 @@ type
     procedure initialisiereDichte(ort,tlc: longint; breite,n: double); inline;
     procedure pruefeArrayEnden(fehler: string);
     procedure nichtnegativieren;
+    procedure berechnePhasenraumAbleitungen; inline;
   public
     emFelder:  array[tEMFeldGroesze,boolean] of pDouble;         // EM-Felder und deren Zeitableitungen
                                                                    // dPhiDX und B[xyz] haben keine sinnvolle Ableitung hier!
@@ -134,7 +135,7 @@ type
     {$DEFINE LiKoInterface}
     {$INCLUDE linearkombinationen.inc}
     {$UNDEF LiKoInterface}
-    procedure schreibePhasenraum(tlc: longint; datNam: string);
+    procedure schreibePhasenraum(tlc: longint; datNamPre, datNamSuf: string);
   end;
 
   { tGitter }
@@ -887,36 +888,7 @@ begin
 
   // Gradienten der Phasenraumbesetzungsdichte berechnen
 
-  pruefeArrayEnden('Jemand hat das FFT-Array überschrieben!');
-
-  for i:=0 to length(ffts)-1 do begin
-    fftw_execute(ffts[i,0,false]);     // FFT in x
-    pruefeArrayEnden('Vorwärts-FFT ('+inttostr(i)+',0) schreibt beim Ausführen über das Ende hinaus!');
-
-    for k:=0 to aP-1 do
-      for j:=0 to aX div 2 do begin  //    *-i*k*2*pi;
-        tmp:=(fftTmp+j+k*(aX div 2 + 1))^.re;
-        (fftTmp+j+k*(aX div 2 + 1))^.re:=  2*pi/dX/aX * j * (fftTmp+j+k*(aX div 2 + 1))^.im/aX * byte(3*j<aX);
-        (fftTmp+j+k*(aX div 2 + 1))^.im:=- 2*pi/dX/aX * j * tmp/aX                             * byte(3*j<aX);
-      end;
-
-    fftw_execute(ffts[i,0,true]);
-    pruefeArrayEnden('Rückwärts-FFT ('+inttostr(i)+',0) schreibt beim Ausführen über das Ende hinaus!');
-
-    fftw_execute(ffts[i,1,false]);     // FFT in p
-    pruefeArrayEnden('Vorwärts-FFT ('+inttostr(i)+',1) schreibt beim Ausführen über das Ende hinaus!');
-
-    for k:=0 to aP div 2 do
-      for j:=0 to aX-1 do begin  //    *-i*k*2*pi;
-        tmp:=(fftTmp+j+k*aX)^.re;
-        (fftTmp+j+k*aX)^.re:=  2*pi/dP/aP * k * (fftTmp+j+k*aX)^.im/aP * byte(4*k < aP);
-        (fftTmp+j+k*aX)^.im:=- 2*pi/dP/aP * k * tmp/aP                 * byte(4*k < aP);
-      end;
-
-    fftw_execute(ffts[i,1,true]);
-    pruefeArrayEnden('Rückwärts-FFT ('+inttostr(i)+',1) schreibt beim Ausführen über das Ende hinaus!');
-
-  end;
+  berechnePhasenraumAbleitungen;
 
   // Zeitableitung der Phasenraumbesetzungsdichte berechnen
 
@@ -1072,12 +1044,13 @@ end;
 {$INCLUDE linearkombinationen.inc}
 {$UNDEF LiKoImplementation}
 
-procedure tFelder.schreibePhasenraum(tlc: longint; datNam: string);
+procedure tFelder.schreibePhasenraum(tlc: longint; datNamPre, datNamSuf: string);
 var
   f:   textfile;
   i,j: longint;
+  typ: byte;
 begin
-  assignfile(f,datNam);
+  assignfile(f,datNamPre+datNamSuf);
   rewrite(f);
   for i:=0 to aP-1 do begin
     for j:=0 to aX-1 do begin
@@ -1087,6 +1060,21 @@ begin
     writeln(f);
   end;
   closefile(f);
+
+  berechnePhasenraumAbleitungen;
+
+  for typ:=0 to 1 do begin
+    assignfile(f,datNamPre+'_d'+char(ord('x')+(ord('p')-ord('x'))*typ)+datNamSuf);
+    rewrite(f);
+    for i:=0 to aP-1 do begin
+      for j:=0 to aX-1 do begin
+        if j>0 then write(f,#9);
+        write(f,(impulsRaumGradient[tlc,typ]+i+j*aP)^);
+      end;
+      writeln(f);
+    end;
+    closefile(f);
+  end;
 end;
 
 procedure tFelder.nichtnegativieren; // Dichten nicht negativ machen
@@ -1139,6 +1127,43 @@ begin
       end;
 end;
 
+procedure tFelder.berechnePhasenraumAbleitungen;
+var
+  i,j,k: longint;
+  tmp:   double;
+begin
+  pruefeArrayEnden('Jemand hat das FFT-Array überschrieben!');
+
+  for i:=0 to length(ffts)-1 do begin
+    fftw_execute(ffts[i,0,false]);     // FFT in x
+    pruefeArrayEnden('Vorwärts-FFT ('+inttostr(i)+',0) schreibt beim Ausführen über das Ende hinaus!');
+
+    for k:=0 to aP-1 do
+      for j:=0 to aX div 2 do begin  //    *-i*k*2*pi;
+        tmp:=(fftTmp+j+k*(aX div 2 + 1))^.re;
+        (fftTmp+j+k*(aX div 2 + 1))^.re:=  2*pi/dX/aX * j * (fftTmp+j+k*(aX div 2 + 1))^.im/aX * byte(3*j<aX);
+        (fftTmp+j+k*(aX div 2 + 1))^.im:=- 2*pi/dX/aX * j * tmp/aX                             * byte(3*j<aX);
+      end;
+
+    fftw_execute(ffts[i,0,true]);
+    pruefeArrayEnden('Rückwärts-FFT ('+inttostr(i)+',0) schreibt beim Ausführen über das Ende hinaus!');
+
+    fftw_execute(ffts[i,1,false]);     // FFT in p
+    pruefeArrayEnden('Vorwärts-FFT ('+inttostr(i)+',1) schreibt beim Ausführen über das Ende hinaus!');
+
+    for k:=0 to aP div 2 do
+      for j:=0 to aX-1 do begin  //    *-i*k*2*pi;
+        tmp:=(fftTmp+j+k*aX)^.re;
+        (fftTmp+j+k*aX)^.re:=  2*pi/dP/aP * k * (fftTmp+j+k*aX)^.im/aP * byte(4*k < aP);
+        (fftTmp+j+k*aX)^.im:=- 2*pi/dP/aP * k * tmp/aP                 * byte(4*k < aP);
+      end;
+
+    fftw_execute(ffts[i,1,true]);
+    pruefeArrayEnden('Rückwärts-FFT ('+inttostr(i)+',1) schreibt beim Ausführen über das Ende hinaus!');
+
+  end;
+end;
+
 // tGitter *********************************************************************
 
 constructor tGitter.create(besitzer: tSimulation; aX,aP: longint; deltaX,deltaP: double; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; name: string);
@@ -1200,7 +1225,7 @@ var
 begin
   if not abbruch then
     for i:=0 to length(felders[aktuelleFelder].matFelder)-1 do
-      felders[aktuelleFelder].schreibePhasenraum(i,extractfilepath(paramstr(0))+'phasenraum_'+inttostr(i+1)+'.dump');
+      felders[aktuelleFelder].schreibePhasenraum(i,extractfilepath(paramstr(0))+'phasenraum_'+inttostr(i+1),'.dump');
   abbruch:=true;
 end;