Phasenraum Schnappschüsse sind nun korrekt orientiert

1 files changed, 99 insertions, 22 deletions

diff --git a/Physikunit.pas b/Physikunit.pas
index b5000e2..e72c48a 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -127,6 +127,8 @@ type
     _dX,_dP,iDX:        double;
     ffts:               array of array[0..1,boolean] of fftw_plan_double; // fft-Algorithmen
     fftTmp,lFftTmp:     Pcomplex_double;                         // Zwischenspeicher für ffts
+    transps:            array of fftw_plan_double;               // Algorithmen zum Transponieren des Phasenraumes
+    transpTmp:          pDouble;                                 // Zwischenspeicher dafür
     procedure initialisiereDichte(ort,tlc: longint; breite,n: double); inline;
   {$IFDEF exzessiveArrayBereichsTests}
     procedure pruefeArrayEnden(fehler: string);
@@ -157,6 +159,7 @@ type
     {$INCLUDE linearkombinationen.inc}
     {$UNDEF LiKoInterface}
     procedure schreibePhasenraum(tlc: longint; datNamPre, datNamSuf: string);
+    function transponiereFuer(tlc: longint): pDouble;
   end;
 
   { tGitter }
@@ -191,6 +194,7 @@ type
     fortschrittsAnzeige:       boolean;
     feldAusgabeDateien:        array of tFeldAusgabeDatei;
     impulsRaumAusgabeDateien:  array of tImpulsRaumAusgabeDatei;
+    transpTmp:                 array of double;
   public
     gotSigusr1,gotSigterm,
     gotSigint:                 boolean;
@@ -557,21 +561,22 @@ begin
 
   assignFile(datei,pre+'-'+inttostr(nNum-1)+'-'+inttostr(tCnt)+suf);
   rewrite(datei,1);
+  inc(tCnt);
 
   sF:=gitter.felders[gitter.aktuelleFelder];
 
-  blockWrite(datei,gitter.xl,sizeof(double));
-  tmpD:=gitter.xl + (sF.aX-1)*sF.dX;
-  blockWrite(datei,tmpD,sizeof(double));
-  tmpI:=sF.aX;
-  blockWrite(datei,tmpI,sizeof(longint));
   tmpD:=-(sF.aP-1)*sF.dP/2;
   blockWrite(datei,tmpD,sizeof(double));
   tmpD:=-tmpD;
   blockWrite(datei,tmpD,sizeof(double));
   tmpI:=sF.aP;
   blockWrite(datei,tmpI,sizeof(longint));
-  blockWrite(datei,sF.impulsraum[teilchen,false]^,sizeof(double)*sF.aX*sF.aP);
+  blockWrite(datei,gitter.xl,sizeof(double));
+  tmpD:=gitter.xl + (sF.aX-1)*sF.dX;
+  blockWrite(datei,tmpD,sizeof(double));
+  tmpI:=sF.aX;
+  blockWrite(datei,tmpI,sizeof(longint));
+  blockWrite(datei,sF.transponiereFuer(teilchen)^,sizeof(double)*sF.aX*sF.aP);
 
   closeFile(datei);
 end;
@@ -742,6 +747,7 @@ var
   emF:        tEMFeldGroesze;
   emQ:        tEMQuellGroesze;
   maF:        tMaterieFeldGroesze;
+  dims:       Pfftw_iodim;
 begin
   inherited create;
   gitter:=parent;
@@ -871,25 +877,60 @@ begin
       );
   end;
 
+    (*
+       fftw_plan_many_dft_r2c(int rank, const int *n, int howmany,
+                                        double *in, const int *inembed,
+                                        int istride, int idist,
+                                        fftw_complex *out, const int *onembed,
+                                        int ostride, int odist,
+                                        unsigned flags);
+       fftw_plan fftw_plan_many_dft_c2r(int rank, const int *n, int howmany,
+                                        fftw_complex *in, const int *inembed,
+                                        int istride, int idist,
+                                        double *out, const int *onembed,
+                                        int ostride, int odist,
+                                        unsigned flags);
+    *)
+
+  fillchar(transps,sizeof(transps),#0);
+  setlength(transps,length(teilchen));
+  fftw_getmem(transpTmp,(aX*aP+1)*sizeof(double));
+  (transpTmp+aX*aP)^:=pi;
+  getmem(dims,2*sizeof(fftw_iodim));
+  dims^.n:=aX;       // aX Spalten
+  dims^.ins:=aP;     // Werte im Input um jeweils aP versetzt
+  dims^.ous:=1;      // und im Output um jeweils 1
+  (dims+1)^.n:=aP;   // aP Zeilen
+  (dims+1)^.ins:=1;  // Werte im Input um jeweils 1 versetzt
+  (dims+1)^.ous:=aX; // und im Output um jeweils aX
+  for i:=0 to length(transps)-1 do begin
+    transps[i]:=
+      fftw_plan_guru_r2r(
+        0,                   // null Dimensionen für die FFT, also nur kopieren!
+        nil,                 // Bockwurst
+        2,                   // 2 Dimensionen
+        dims,                // @(aX;aP)
+        impulsraum[i,false], // von hier
+        transpTmp,           // nach dort
+        nil,                 // Bockwurst
+        [fftw_estimate]      // einfach irgendwas nehmen
+      );
+  end;
+  freemem(dims);
+
+    (*
+       fftw_plan fftw_plan_guru_r2r(int rank, const fftw_iodim *dims,
+                                    int howmany_rank,
+                                    const fftw_iodim *howmany_dims,
+                                    double *in, double *out,
+                                    const fftw_r2r_kind *kind,
+                                    unsigned flags);
+    *)
+
   {$IFDEF exzessiveArrayBereichsTests}
-  pruefeArrayEnden('FFT schreibt beim Planen über das Ende hinaus!');
+  pruefeArrayEnden('FFT/Transponieren schreibt beim Planen über das Ende hinaus!');
   {$ENDIF}
 
-  (*
-     fftw_plan_many_dft_r2c(int rank, const int *n, int howmany,
-                                      double *in, const int *inembed,
-                                      int istride, int idist,
-                                      fftw_complex *out, const int *onembed,
-                                      int ostride, int odist,
-                                      unsigned flags);
-     fftw_plan fftw_plan_many_dft_c2r(int rank, const int *n, int howmany,
-                                      fftw_complex *in, const int *inembed,
-                                      int istride, int idist,
-                                      double *out, const int *onembed,
-                                      int ostride, int odist,
-                                      unsigned flags);
-  *)
-
   fillchar(massen,sizeof(massen),#0);
   setlength(massen,length(teilchen));
   for i:=0 to length(massen)-1 do
@@ -950,6 +991,10 @@ begin
         fftw_destroy_plan(ffts[i,j,abl]);
   setlength(ffts,0);
   fftw_freemem(fftTmp);
+  for i:=0 to length(transps)-1 do
+    fftw_destroy_plan(transps[i]);
+  setlength(transps,0);
+  fftw_freemem(transpTmp);
   inherited destroy;
 end;
 
@@ -974,6 +1019,11 @@ var
   i,j: longint;
   abl: boolean;
 begin
+  if abs((transpTmp+aP*aX)^-pi)>1e-10 then begin
+    gitter.prot.schreibe(fehler);
+    gitter.prot.schreibe(floattostr((transpTmp+aP*aX)^)+' vs. '+floattostr(pi));
+    raise exception.create('Fehler in tFelder.pruefeArrayEnden!');
+  end;
   if (abs(lFftTmp^.re-pi)>1e-10) or (abs(lFftTmp^.im-sqr(pi))>1e-10) then begin
     gitter.prot.schreibe(fehler);
     gitter.prot.schreibe(floattostr(lFftTmp^.re)+' vs. '+floattostr(pi)+' und '+floattostr(lFftTmp^.im)+' vs. '+floattostr(sqr(pi)));
@@ -1423,6 +1473,29 @@ begin
   end;
 end;
 
+function tFelder.transponiereFuer(tlc: longint): pDouble;
+{$IFDEF exzessiveArrayBereichsTests}
+var
+  i: longint;
+{$ENDIF}
+begin
+  {$IFDEF exzessiveArrayBereichsTests}
+  for i:=0 to aP*aX-1 do
+    (transpTmp+i)^:=sqrt(2);
+  {$ENDIF}
+  fftw_execute(transps[tlc]);
+  {$IFDEF exzessiveArrayBereichsTests}
+  pruefeArrayEnden('TransporniereFuer('+inttostr(tlc)+') schreibt über das Ende hinaus!');
+  for i:=0 to aP*aX-1 do
+    if (transpTmp+i)^=sqrt(2) then begin
+      gitter.prot.schreibe('TransporniereFuer('+inttostr(tlc)+') schreibt nicht alle Werte!');
+      gitter.prot.schreibe('Position '+inttostr(i)+' wurde ausgelassen!');
+      raise exception.create('Fehler in tFelder.transponiereFuer!');
+    end;
+  {$ENDIF}
+  result:=transpTmp;
+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);
@@ -1608,6 +1681,9 @@ begin
   setlength(feldAusgabeDateien,0);
   setlength(impulsRaumAusgabeDateien,0);
 
+  fillchar(transpTmp,sizeof(transpTmp),#0); // Transpositions-Puffer initialisieren
+  setlength(transpTmp,0);
+
   // Standardeinstellungen Bereich 'lichtVon...'
   lichter:=tMyStringlist.create(prot,'lichter');
 
@@ -1941,6 +2017,7 @@ begin
         simulationen[j-1]:=simulationen[j];
       setlength(simulationen,length(simulationen)-1);
     end;
+  setlength(transpTmp,0);
   inherited destroy;
 end;