extended -> double scheint viele Probleme zu lösen :-/

3 files changed, 176 insertions, 207 deletions

diff --git a/Physikunit.pas b/Physikunit.pas
index b49675b..f95dc32 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -13,11 +13,11 @@ unit Physikunit;
 interface
 
 uses
-  Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit, baseUnix, fftw_l;
+  Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit, baseUnix, fftw;
 
 type
   tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaDreiAchtel,zfRungeKuttaVier,zfRungeKuttaZehn,zfRungeKuttaZwoelf,zfRungeKuttaVierzehn);
-  tVerteilungsfunktion = function(x: extended): extended;
+  tVerteilungsfunktion = function(x: double): double;
   tEMFeldGroesze = (
     efAX,efAY,
     efDAXDT,efDAYDT,
@@ -29,7 +29,7 @@ type
   tMaterieFeldGroesze = (mfPY);
   tMaterieSpeicherGroesze = (msPX,msPXSqr,msPY,msVX,msVY,msN,msPXRipple);
   tTeilchenSpeziesGroeszen = (tsgMasse,tsgIQdrMasse,tsgLadung);
-  tEMQuellen = array[tEMQuellGroesze] of extended;
+  tEMQuellen = array[tEMQuellGroesze] of double;
 
 const
   erstesEMFMitAbleitung  =  efAX;
@@ -60,36 +60,36 @@ type
     procedure schreibeBuffer(minBufLen: longint);
   public
     zeitAnz:            longint;
-    sDT:                extended;
+    sDT:                double;
     constructor create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string);
     destructor destroy; override;
     function dump: string;
     procedure schreibeKopf;
-    procedure speichereWerte(gitter: tGitter; sT, sDX: extended);
+    procedure speichereWerte(gitter: tGitter; sT, sDX: double);
   end;
 
   { tTeilchenSpezies }
 
   tTeilchenSpezies = class
   private
-    ortsGrenzen:     array of extended;
+    ortsGrenzen:     array of double;
     dichteStuecke:   array of string;
   public
-    nMin,nMax,breite: extended;                                    // minimale und maximale Massendichte in nc; Breite der Impulsverteilung
-    eigenschaften:    array[tTeilchenSpeziesGroeszen] of extended; // q in e; m in me
-    diffusion:        array[0..1] of extended;                     // Diffusionsterm in x- und p-Richtung (= maximales d(ln n)/d[xp] * iD[XP])
+    nMin,nMax,breite: double;                                    // minimale und maximale Massendichte in nc; Breite der Impulsverteilung
+    eigenschaften:    array[tTeilchenSpeziesGroeszen] of double; // q in e; m in me
+    diffusion:        array[0..1] of double;                     // Diffusionsterm in x- und p-Richtung (= maximales d(ln n)/d[xp] * iD[XP])
     constructor create; overload;
     constructor create(original: tTeilchenSpezies); overload;
     destructor destroy; override;
     procedure dump(prot: tProtokollant; prefix: string);
-    function gibDichte(x: extended; kvs: tKnownValues): extended; // Massendichte
+    function gibDichte(x: double; kvs: tKnownValues): double; // Massendichte
     procedure liesDichteFunktion(rest: string; ifile: tMyStringList; kvs: tKnownValues; teilchen: array of tTeilchenSpezies; prot: tProtokollant);
   end;
 
   { TFelder }
 (*
-  procedure berechneAbleitungen(pX: extended); inline;
-  function nichtnegativieren: extended; inline;
+  procedure berechneAbleitungen(pX: double); inline;
+  function nichtnegativieren: double; inline;
 *)
 
   tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen
@@ -99,35 +99,35 @@ type
 
     teilchen:           array of tTeilchenSpezies;
     lichters:           array[boolean] of tMyStringlist;
-    massen:             array of extended;                         // Gesamtmassen je Teilchenspezies zur Überwachung der Erhaltungsgrößen
-    gesamtDefizit:      extended;
-    impulsRaumGradient: array of array[0..1] of pExtended;         // Ableitung des Impulsraumes nach x und px
-    iMGamma:            array of pExtended;                        // 1/m/gamma (= v/p)
+    massen:             array of double;                         // Gesamtmassen je Teilchenspezies zur Überwachung der Erhaltungsgrößen
+    gesamtDefizit:      double;
+    impulsRaumGradient: array of array[0..1] of pDouble;         // Ableitung des Impulsraumes nach x und px
+    iMGamma:            array of pDouble;                        // 1/m/gamma (= v/p)
     _aX,_aP:            longint;
-    _dX,_dP,iDX:        extended;
-    ffts:               array of array[0..1,boolean] of fftw_plan_extended; // fft-Algorithmen
-    fftTmp:             Pcomplex_extended;                         // Zwischenspeicher für ffts
-    procedure initialisiereDichte(ort,tlc: longint; breite,n: extended); inline;
+    _dX,_dP,iDX:        double;
+    ffts:               array of array[0..1,boolean] of fftw_plan_double; // fft-Algorithmen
+    fftTmp:             Pcomplex_double;                         // Zwischenspeicher für ffts
+    procedure initialisiereDichte(ort,tlc: longint; breite,n: double); inline;
   public
-    emFelder:  array[tEMFeldGroesze,boolean] of pExtended;         // EM-Felder und deren Zeitableitungen
+    emFelder:  array[tEMFeldGroesze,boolean] of pDouble;         // EM-Felder und deren Zeitableitungen
                                                                    // dPhiDX und B[xyz] haben keine sinnvolle Ableitung hier!
-    emQuellen: array[tEMQuellGroesze] of pExtended;                // Quellen für die EM-Felder
-    matFelder: array of array[tMaterieFeldGroesze] of pExtended;   // px-unabhängige Felder, nach Teilchenspezies sortiert
-    impulsraum: array of array[boolean] of pExtended;              // px-abhängige Felder, nach Teilchenspezies sortiert, sowie deren Ableitung
+    emQuellen: array[tEMQuellGroesze] of pDouble;                // Quellen für die EM-Felder
+    matFelder: array of array[tMaterieFeldGroesze] of pDouble;   // px-unabhängige Felder, nach Teilchenspezies sortiert
+    impulsraum: array of array[boolean] of pDouble;              // px-abhängige Felder, nach Teilchenspezies sortiert, sowie deren Ableitung
 
     gitter: tGitter;
 
-    constructor create(anzX,anzP: longint; deltaX,deltaP: extended; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
+    constructor create(anzX,anzP: longint; deltaX,deltaP: double; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
     destructor destroy; override;
     procedure berechneAbleitungen;
     procedure setzeNull; inline;
     procedure dumpErhaltungsgroeszen(prot: tProtokollant);
     property aX: longint read _aX;
     property aP: longint read _aP;
-    property dX: extended read _dX;
-    property dP: extended read _dP;
-    function impulsIntegral(ort,tlc: longint; emQ: tEMQuellGroesze): extended; overload; inline;
-    function impulsIntegral(ort,tlc: longint; maF: tMaterieSpeicherGroesze): extended; overload; inline;
+    property dX: double read _dX;
+    property dP: double read _dP;
+    function impulsIntegral(ort,tlc: longint; emQ: tEMQuellGroesze): double; overload; inline;
+    function impulsIntegral(ort,tlc: longint; maF: tMaterieSpeicherGroesze): double; overload; inline;
     {$DEFINE LiKoInterface}
     {$INCLUDE linearkombinationen.inc}
     {$UNDEF LiKoInterface}
@@ -141,7 +141,7 @@ type
     prot:                         tProtokollant;
     zeitverfahren:                tZeitverfahren;
     abbruch:                      boolean;
-    dX,iDX,xl,t:                  extended;
+    dX,iDX,xl,t:                  double;
     kvs:                          tKnownValues;
     procedure abbrechen;
   public
@@ -149,9 +149,9 @@ type
     aktuelleFelder: longint;
     felders: array of tFelder; // mehrere komplette Simulationsboxen von Feldern, benötigt um Zwischenschritte für die Zeitentwicklung zu speichern
 
-    constructor create(derBesitzer: tSimulation; aX,aP: longint; deltaX,deltaP: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string);
+    constructor create(derBesitzer: tSimulation; aX,aP: longint; deltaX,deltaP: double; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string);
     destructor destroy; override;
-    procedure iteriereSchritt(dT: extended);
+    procedure iteriereSchritt(dT: double);
     procedure dumpErhaltungsgroeszen;
   end;
 
@@ -161,14 +161,14 @@ type
   private
     prot:                      tProtokollant;
     gitter:                    tGitter;
-    dT,tEnde,sT,sDT,sDX:       extended;
+    dT,tEnde,sT,sDT,sDX:       double;
     fortschrittsAnzeige:       boolean;
     ausgabeDateien:            array of tAusgabeDatei;
   public
     gotSighup:                 boolean;
     constructor create(inName: string; protokollant: tProtokollant; name: string);
     destructor destroy; override;
-    function iteriereSchritt(start: extended; var zeitPhysik,zeitDatei: extended): boolean; // noch nicht zu Ende?
+    function iteriereSchritt(start: double; var zeitPhysik,zeitDatei: double): boolean; // noch nicht zu Ende?
   end;
 
 procedure SignalCapture(signal : longint); cdecl;
@@ -349,9 +349,9 @@ begin
   closefile(datei);
 end;
 
-procedure tAusgabeDatei.speichereWerte(gitter: tGitter; sT, sDX: extended);
+procedure tAusgabeDatei.speichereWerte(gitter: tGitter; sT, sDX: double);
 var i,cnt:     longint;
-    sX,cX,val: extended;
+    sX,cX,val: double;
 begin
   if (teilchen>=length(gitter.felders[gitter.aktuelleFelder].teilchen)) then begin
     pro.schreibe('Teilchen '+inttostr(teilchen+1)+' gibt es nicht, da kann ich auch nichts speichern!',true);
@@ -367,19 +367,19 @@ begin
   sX:=cX+(cnt-1)*sDX;
 
   inc(tCnt);
-  if bufLen-bufPos < sizeof(integer) + (2+cnt)*sizeof(extended) then
-    schreibeBuffer(sizeof(integer) + (2+cnt)*sizeof(extended));
+  if bufLen-bufPos < sizeof(integer) + (2+cnt)*sizeof(double) then
+    schreibeBuffer(sizeof(integer) + (2+cnt)*sizeof(double));
 
-  if bufLen-bufPos < sizeof(integer) + (2+cnt)*sizeof(extended) then begin
-    pro.schreibe('Schreibbuffer ist zu klein ('+inttostr(bufLen)+' statt mindestens '+inttostr(sizeof(integer) + (2+cnt)*sizeof(extended))+' Bytes)!');
+  if bufLen-bufPos < sizeof(integer) + (2+cnt)*sizeof(double) then begin
+    pro.schreibe('Schreibbuffer ist zu klein ('+inttostr(bufLen)+' statt mindestens '+inttostr(sizeof(integer) + (2+cnt)*sizeof(double))+' Bytes)!');
     pro.destroyall;
     halt(1);
   end;
 
-  move(cX,(buf+bufPos)^,sizeof(extended));
-  inc(bufPos,sizeof(extended));
-  move(sX,(buf+bufPos)^,sizeof(extended));
-  inc(bufPos,sizeof(extended));
+  move(cX,(buf+bufPos)^,sizeof(double));
+  inc(bufPos,sizeof(double));
+  move(sX,(buf+bufPos)^,sizeof(double));
+  inc(bufPos,sizeof(double));
   move(cnt,(buf+bufPos)^,sizeof(integer));
   inc(bufPos,sizeof(integer));
 
@@ -389,8 +389,8 @@ begin
       for i:=0 to gitter.felders[gitter.aktuelleFelder].aX-1 do begin
         while cX>=sX do begin
           dec(cnt);
-          move((gitter.felders[gitter.aktuelleFelder].emFelder[emFeld,ableitung]+i)^,(buf+bufPos)^,sizeof(extended));
-          inc(bufPos,sizeof(extended));
+          move((gitter.felders[gitter.aktuelleFelder].emFelder[emFeld,ableitung]+i)^,(buf+bufPos)^,sizeof(double));
+          inc(bufPos,sizeof(double));
           sX:=sX+sDX;
         end;
         cX:=cX+gitter.dX;
@@ -400,8 +400,8 @@ begin
         while cX>=sX do begin
           dec(cnt);
           val:=gitter.felders[gitter.aktuelleFelder].impulsIntegral(i,teilchen,emQuelle);
-          move(val,(buf+bufPos)^,sizeof(extended));
-          inc(bufPos,sizeof(extended));
+          move(val,(buf+bufPos)^,sizeof(double));
+          inc(bufPos,sizeof(double));
           sX:=sX+sDX;
         end;
         cX:=cX+gitter.dX;
@@ -411,8 +411,8 @@ begin
         while cX>=sX do begin
           dec(cnt);
           val:=gitter.felders[gitter.aktuelleFelder].impulsIntegral(i,teilchen,matFeld);
-          move(val,(buf+bufPos)^,sizeof(extended));
-          inc(bufPos,sizeof(extended));
+          move(val,(buf+bufPos)^,sizeof(double));
+          inc(bufPos,sizeof(double));
           sX:=sX+sDX;
         end;
         cX:=cX+gitter.dX;
@@ -484,7 +484,7 @@ begin
   end;
 end;
 
-function tTeilchenSpezies.gibDichte(x: extended; kvs: tKnownValues): extended;
+function tTeilchenSpezies.gibDichte(x: double; kvs: tKnownValues): double;
 var
   i: longint;
 begin
@@ -547,7 +547,7 @@ begin
 end;
 
 (*
-function tRaumPunkt.nichtnegativieren: extended; // Dichten nicht negativ machen
+function tRaumPunkt.nichtnegativieren: double; // Dichten nicht negativ machen
 var
   i: longint;
 begin
@@ -559,11 +559,11 @@ end;
 *)
 // tFelder *********************************************************************
 
-constructor tFelder.create(anzX,anzP: longint; deltaX,deltaP: extended; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
+constructor tFelder.create(anzX,anzP: longint; deltaX,deltaP: double; _teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
 var
   i,j:        longint;
   rechts,abl: boolean;
-  dens:       extended;
+  dens:       double;
   emF:        tEMFeldGroesze;
   emQ:        tEMQuellGroesze;
   maF:        tMaterieFeldGroesze;
@@ -585,39 +585,39 @@ begin
   setlength(impulsRaumGradient,length(teilchen));
   for i:=0 to length(impulsRaumGradient)-1 do
     for j:=0 to length(impulsRaumGradient[i])-1 do begin
-      fftw_getmem(impulsRaumGradient[i,j],aP*aX*sizeof(extended));
-      fillchar(impulsRaumGradient[i,j]^,aP*aX*sizeof(extended),#0);
+      fftw_getmem(impulsRaumGradient[i,j],aP*aX*sizeof(double));
+      fillchar(impulsRaumGradient[i,j]^,aP*aX*sizeof(double),#0);
     end;
   fillchar(iMGamma,sizeof(iMGamma),#0);
   setlength(iMGamma,length(teilchen));
   for i:=0 to length(iMGamma)-1 do begin
-    fftw_getmem(iMGamma[i],aP*aX*sizeof(extended));
-    fillchar(iMGamma[i]^,aP*aX*sizeof(extended),#0);
+    fftw_getmem(iMGamma[i],aP*aX*sizeof(double));
+    fillchar(iMGamma[i]^,aP*aX*sizeof(double),#0);
   end;
   for emF:=low(tEMFeldGroesze) to high(tEMFeldGroesze) do begin
     emFelder[emF,true]:=nil;
     for abl:=false to not(emF in [efEX,efBZ]) do begin
-      fftw_getmem(emFelder[emF,abl],aX*sizeof(extended));
-      fillchar(emFelder[emF,abl]^,aX*sizeof(extended),#0);
+      fftw_getmem(emFelder[emF,abl],aX*sizeof(double));
+      fillchar(emFelder[emF,abl]^,aX*sizeof(double),#0);
     end;
   end;
   for emQ:=low(tEMQuellGroesze) to high(tEMQuellGroesze) do begin
-    fftw_getmem(emQuellen[emQ],aX*sizeof(extended));
-    fillchar(emQuellen[emQ]^,aX*sizeof(extended),#0);
+    fftw_getmem(emQuellen[emQ],aX*sizeof(double));
+    fillchar(emQuellen[emQ]^,aX*sizeof(double),#0);
   end;
   fillchar(matFelder,sizeof(matFelder),#0);
   setlength(matFelder,length(teilchen));
   for i:=0 to length(matFelder)-1 do
     for maF:=low(tMaterieFeldGroesze) to high(tMaterieFeldGroesze) do begin
-      fftw_getmem(matFelder[i,maF],aX*sizeof(extended));
-      fillchar(matFelder[i,maF]^,aX*sizeof(extended),#0);
+      fftw_getmem(matFelder[i,maF],aX*sizeof(double));
+      fillchar(matFelder[i,maF]^,aX*sizeof(double),#0);
     end;
   fillchar(impulsraum,sizeof(impulsraum),#0);
   setlength(impulsraum,length(teilchen));
   for i:=0 to length(impulsraum)-1 do
     for abl:=false to true do begin
-      fftw_getmem(impulsraum[i,abl],aX*aP*sizeof(extended));
-      fillchar(impulsraum[i,abl]^,aX*aP*sizeof(extended),#0);
+      fftw_getmem(impulsraum[i,abl],aX*aP*sizeof(double));
+      fillchar(impulsraum[i,abl]^,aX*aP*sizeof(double),#0);
     end;
 
   fillchar(massen,sizeof(massen),#0);
@@ -639,21 +639,22 @@ begin
   lichters[true].grep('^rechts .*');
   lichters[true].subst('^rechts *','');
 
-  fftw_getmem(fftTmp,(aP+1)*(aX+1)*sizeof(complex_extended));
+  fftw_getmem(fftTmp,(aP+1)*(aX+1)*sizeof(complex_double));
+  fillchar(fftTmp^,(aP+1)*(aX+1)*sizeof(complex_double),#0);
 
   fillchar(ffts,sizeof(ffts),#0);
-  setlength(ffts,length(teilchen));  write('.');
+  setlength(ffts,length(teilchen));
   for i:=0 to length(ffts)-1 do begin
-//    ffts[i,0,false]:=  // Planung der Hintransformationen über x
-//      fftw_plan_many_dft_r2c(1,@_aX,aP,impulsraum[i,false],nil,aP,1,fftTmp,nil,aP,1,[fftw_preserve_input,fftw_measure]); write('1');
+    ffts[i,0,false]:=  // Planung der Hintransformationen über x
+      fftw_plan_many_dft_r2c(1,@_aX,aP,impulsraum[i,false],nil,aP,1,fftTmp,nil,aP,1,[fftw_preserve_input,fftw_measure]);
     ffts[i,0,true]:=  // Planung der Rücktransformationen über x
-      fftw_plan_many_dft_c2r(1,@_aX,aP,fftTmp,nil,aP,1,impulsRaumGradient[i,0],nil,aP,1,[fftw_measure]);                 write('2');
+      fftw_plan_many_dft_c2r(1,@_aX,aP,fftTmp,nil,aP,1,impulsRaumGradient[i,0],nil,aP,1,[fftw_measure]);
 
-//    ffts[i,1,false]:=  // Planung der Hintransformationen über p
-//      fftw_plan_many_dft_r2c(1,@_aP,aX,impulsraum[i,false],nil,1,aP,fftTmp,nil,1,aP,[fftw_preserve_input,fftw_measure]); write('3');
-//    ffts[i,1,true]:=  // Planung der Rücktransformationen über p
-//      fftw_plan_many_dft_c2r(1,@_aP,aX,fftTmp,nil,1,aP,impulsRaumGradient[i,1],nil,1,aP,[fftw_measure]);                 write('4');
-  end;                               write('>');
+    ffts[i,1,false]:=  // Planung der Hintransformationen über p
+      fftw_plan_many_dft_r2c(1,@_aP,aX,impulsraum[i,false],nil,1,aP,fftTmp,nil,1,aP,[fftw_preserve_input,fftw_measure]);
+    ffts[i,1,true]:=  // Planung der Rücktransformationen über p
+      fftw_plan_many_dft_c2r(1,@_aP,aX,fftTmp,nil,1,aP,impulsRaumGradient[i,1],nil,1,aP,[fftw_measure]);
+  end;
    (*
      fftw_plan_many_dft_r2c(int rank, const int *n, int howmany,
                                       double *in, const int *inembed,
@@ -713,10 +714,10 @@ begin
   inherited destroy;
 end;
 
-procedure tFelder.initialisiereDichte(ort,tlc: longint; breite,n: extended);
+procedure tFelder.initialisiereDichte(ort,tlc: longint; breite,n: double);
 var
   i:   longint;
-  ges: extended;
+  ges: double;
 begin
   ges:=0;
   for i:=0 to aP-1 do begin
@@ -734,7 +735,7 @@ var
   emQ:    tEMQuellGroesze;
   emF:    tEMFeldGroesze;
   rechts: boolean;
-  tmp,err:extended;
+  tmp:    double;
 begin
 
   // PX ggf. korrigieren
@@ -824,35 +825,32 @@ begin
   // dA/dt = dA/dt - E
   for i:=1 to aX-2 do
     (emFelder[efAX,true]+i)^:=(emFelder[efDAXDT,false]+i)^ - (emFelder[efEX,false]+i)^;
-  move(emFelder[efDAYDT,false]^,emFelder[efAY,true]^,aX*sizeof(extended));
+  move(emFelder[efDAYDT,false]^,emFelder[efAY,true]^,aX*sizeof(double));
 
   // Gradienten der Phasenraumbesetzungsdichte berechnen
 
   for i:=0 to length(ffts)-1 do begin
-//    fftw_execute(ffts[i,0,false]);
-  //    *i*k*2*pi;
-//    fftw_execute(ffts[i,0,true]);
-    write('A');
+    fftw_execute(ffts[i,0,false]);
+    for j:=0 to aX-1 do
+      for k:=0 to aP-1 do begin  //    *i*k*2*pi;
+        tmp:=(fftTmp+j+k*aX)^.re;
+        (fftTmp+j+k*aX)^.re:=- 2*pi/dX/aX * j * (fftTmp+j+k*aX)^.im/aX;
+        (fftTmp+j+k*aX)^.im:=  2*pi/dX/aX * j * tmp/aX;
+      end;
+
+    fftw_execute(ffts[i,0,true]);
     fftw_execute(ffts[i,1,false]);
-    write('B');
-  //    *i*k*2*pi;
+
+    for j:=0 to aX-1 do
+      for k:=0 to aP-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;
+        (fftTmp+j+k*aX)^.im:=  2*pi/dP/aP * k * tmp/aP;
+      end;
+
     fftw_execute(ffts[i,1,true]);
-    write('C');
   end;
 
-  tmp:=0;
-  err:=0;
-  for i:=0 to length(impulsraum)-1 do
-    for j:=0 to aP*aX-1 do begin
-      tmp:=tmp+abs((impulsraum[i,false]+j)^)+abs((impulsRaumGradient[i,1]+j)^);
-      err:=err+abs((impulsraum[i,false]+j)^-(impulsRaumGradient[i,1]+j)^);
-//      tmp:=tmp+abs((impulsraum[i,false]+j)^)+abs((impulsRaumGradient[i,0]+j)^)+abs((impulsRaumGradient[i,1]+j)^);
-//      err:=err+abs((impulsraum[i,false]+j)^-(impulsRaumGradient[i,0]+j)^)
-//              +abs((impulsraum[i,false]+j)^-(impulsRaumGradient[i,1]+j)^);
-    end;
-  writeln(tmp,' ',err);
-  halt;
-
   // Zeitableitung der Phasenraumbesetzungsdichte berechnen
 
   for i:=0 to length(teilchen)-1 do
@@ -902,7 +900,7 @@ end;
 procedure tFelder.dumpErhaltungsgroeszen(prot: tProtokollant);
 var
   i,j:  longint;
-  dens: extended;
+  dens: double;
 begin
   dens:=0;
   for i:=0 to length(massen)-1 do
@@ -917,7 +915,7 @@ begin
   end;
 end;
 
-function tFelder.impulsIntegral(ort,tlc: longint; emQ: tEMQuellGroesze): extended;
+function tFelder.impulsIntegral(ort,tlc: longint; emQ: tEMQuellGroesze): double;
 begin
   if tlc<0 then begin
     result:=(emQuellen[emQ]+ort)^; // das ist leicht :-)
@@ -934,7 +932,7 @@ begin
   result:=result * teilchen[tlc].eigenschaften[tsgLadung];
 end;
 
-function tFelder.impulsIntegral(ort,tlc: longint; maF: tMaterieSpeicherGroesze): extended;
+function tFelder.impulsIntegral(ort,tlc: longint; maF: tMaterieSpeicherGroesze): double;
 var
   i: longint;
 begin
@@ -999,7 +997,7 @@ end;
 
 // tGitter *********************************************************************
 
-constructor tGitter.create(derBesitzer: tSimulation; aX,aP: longint; deltaX,deltaP: extended; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string);
+constructor tGitter.create(derBesitzer: tSimulation; aX,aP: longint; deltaX,deltaP: double; bekannteWerte: tKnownValues; teilchen: array of tTeilchenSpezies; lichter: tMyStringlist; zv: tZeitverfahren; protokollant: tProtokollant; name: string);
 var
   i: longint;
 begin
@@ -1061,7 +1059,7 @@ begin
     felders[i].setzeNull;
 end;
 
-procedure tGitter.iteriereSchritt(dT: extended);
+procedure tGitter.iteriereSchritt(dT: double);
     {$IFDEF negativeDichteueberwachung}
 var i,j:     longint;
     pro:     tProtokollant;
@@ -1157,7 +1155,7 @@ var
   ifile:                 tMyStringlist;
   zeitverfahren:         tZeitverfahren;
   s,t,aSuffix,aPrefix:   string;
-  dX,breite,dP,pMax:     extended;
+  dX,breite,dP,pMax:     double;
   i:                     longint;
   kvs:                   tKnownValues;
   teilchen:              array of tTeilchenSpezies;
@@ -1520,7 +1518,7 @@ begin
   inherited destroy;
 end;
 
-function tSimulation.iteriereSchritt(start: extended; var zeitPhysik,zeitDatei: extended): boolean; // noch nicht zu Ende?
+function tSimulation.iteriereSchritt(start: double; var zeitPhysik,zeitDatei: double): boolean; // noch nicht zu Ende?
 var
   i: longint;
 begin
diff --git a/Plasmapropagation.lpr b/Plasmapropagation.lpr
index d3b8f0a..e1346c5 100644
--- a/Plasmapropagation.lpr
+++ b/Plasmapropagation.lpr
@@ -8,7 +8,7 @@ uses
   {$ENDIF}{$ENDIF}
   Classes, SysUtils, CustApp,
   { you can add units after this }
-  math, Physikunit, protokollunit, fftw_l;
+  math, Physikunit, protokollunit, lowlevelunit;
 
 type
 
@@ -27,49 +27,10 @@ type
 procedure TPlasmapropagation.DoRun;
 var
   simulation:                 tSimulation;
-  start,zeitPhysik,zeitDatei: extended;
+  start,zeitPhysik,zeitDatei: double;
   prot:                       tProtokollant;
   s,t,u:                      string;
-
-  xlen,ylen,i,j:              longint;
-  ein:                        Pcomplex_extended;
-  aus:                        Pextended;
-  plan:                       fftw_plan_extended;
 begin
-
-  xlen:=4;
-  ylen:=4;
-
-  fftw_getmem(ein,xlen*ylen*sizeof(complex_extended));
-  fftw_getmem(aus,xlen*ylen*sizeof(extended));
-
-  plan:=fftw_plan_many_dft_c2r(1,@xlen,ylen,ein,nil,ylen,1,aus,nil,ylen,1,[fftw_measure]);
-
-  for i:=0 to xlen-1 do
-    for j:=0 to ylen-1 do begin
-      (ein+j+i*ylen)^.re:=j;
-      (ein+j+i*ylen)^.im:=i;
-    end;
-  for i:=0 to xlen-1 do begin
-    for j:=0 to ylen-1 do
-      write((ein+j+i*ylen)^.re,' + ',(ein+j+i*ylen)^.im,' *I ; ');
-    writeln;
-  end;
-
-  fftw_execute(plan);
-
-  for i:=0 to xlen-1 do begin
-    for j:=0 to ylen-1 do
-      write((aus+j+i*ylen)^,' ');
-    writeln;
-  end;
-
-  fftw_destroy_plan(plan);
-  fftw_freemem(ein);
-  fftw_freemem(aus);
-
-  halt;
-
   prot:=tProtokollant.create('error');
 
   if paramcount<>1 then begin
diff --git a/Plasmapropagation.lps b/Plasmapropagation.lps
index eee3207..d9cdf52 100644
--- a/Plasmapropagation.lps
+++ b/Plasmapropagation.lps
@@ -7,19 +7,19 @@
       <Unit0>
         <Filename Value="Plasmapropagation.lpr"/>
         <IsPartOfProject Value="True"/>
-        <IsVisibleTab Value="True"/>
-        <TopLine Value="17"/>
-        <CursorPos X="49" Y="37"/>
+        <TopLine Value="19"/>
+        <CursorPos Y="46"/>
         <UsageCount Value="200"/>
         <Loaded Value="True"/>
       </Unit0>
       <Unit1>
         <Filename Value="Physikunit.pas"/>
         <IsPartOfProject Value="True"/>
+        <IsVisibleTab Value="True"/>
         <EditorIndex Value="1"/>
-        <TopLine Value="626"/>
-        <CursorPos Y="645"/>
-        <FoldState Value=" T3lM0{F5 pjdmK0G51313]D2uD0f2111 p0sU0G2 T0!2Pc071F311'"/>
+        <TopLine Value="618"/>
+        <CursorPos X="59" Y="643"/>
+        <FoldState Value=" T3lM0{F5 pjdmK0G51214]DWuB0f3]9BkH031 pp4pD0G2 T0!0Pc071F311)"/>
         <UsageCount Value="200"/>
         <Loaded Value="True"/>
       </Unit1>
@@ -28,14 +28,14 @@
         <IsPartOfProject Value="True"/>
         <EditorIndex Value="-1"/>
         <CursorPos X="20" Y="8"/>
-        <UsageCount Value="165"/>
+        <UsageCount Value="168"/>
       </Unit2>
       <Unit3>
         <Filename Value="input.plap"/>
         <IsPartOfProject Value="True"/>
-        <EditorIndex Value="4"/>
+        <EditorIndex Value="3"/>
         <CursorPos X="33" Y="23"/>
-        <UsageCount Value="164"/>
+        <UsageCount Value="167"/>
         <Loaded Value="True"/>
         <DefaultSyntaxHighlighter Value="None"/>
       </Unit3>
@@ -44,18 +44,18 @@
         <IsPartOfProject Value="True"/>
         <EditorIndex Value="-1"/>
         <CursorPos X="54" Y="16"/>
-        <UsageCount Value="112"/>
+        <UsageCount Value="115"/>
       </Unit4>
       <Unit5>
         <Filename Value="linearkombinationen.inc"/>
         <IsPartOfProject Value="True"/>
         <EditorIndex Value="-1"/>
         <CursorPos X="29" Y="16"/>
-        <UsageCount Value="75"/>
+        <UsageCount Value="78"/>
       </Unit5>
       <Unit6>
         <Filename Value="input.epost"/>
-        <EditorIndex Value="5"/>
+        <EditorIndex Value="4"/>
         <TopLine Value="22"/>
         <CursorPos X="26" Y="29"/>
         <UsageCount Value="100"/>
@@ -64,19 +64,19 @@
       </Unit6>
       <Unit7>
         <Filename Value="../units/matheunit.pas"/>
-        <EditorIndex Value="3"/>
+        <EditorIndex Value="2"/>
         <TopLine Value="19"/>
         <CursorPos X="66" Y="48"/>
         <FoldState Value=" T3iA05B pj0jV034 piaj60U6 plLmM041p"/>
-        <UsageCount Value="22"/>
+        <UsageCount Value="24"/>
         <Loaded Value="True"/>
       </Unit7>
       <Unit8>
         <Filename Value="../units/lowlevelunit.pas"/>
         <EditorIndex Value="-1"/>
-        <TopLine Value="558"/>
-        <CursorPos X="41" Y="567"/>
-        <UsageCount Value="8"/>
+        <TopLine Value="626"/>
+        <CursorPos X="27" Y="657"/>
+        <UsageCount Value="11"/>
       </Unit8>
       <Unit9>
         <Filename Value="../units/mystringlistunit.pas"/>
@@ -126,146 +126,156 @@
       </Unit15>
       <Unit16>
         <Filename Value="../units/fftw_l.pas"/>
-        <EditorIndex Value="2"/>
-        <TopLine Value="68"/>
-        <CursorPos X="29" Y="96"/>
-        <UsageCount Value="14"/>
-        <Loaded Value="True"/>
+        <EditorIndex Value="-1"/>
+        <TopLine Value="70"/>
+        <CursorPos Y="85"/>
+        <UsageCount Value="15"/>
       </Unit16>
       <Unit17>
         <Filename Value="../units/fftw_s.pas"/>
         <EditorIndex Value="-1"/>
-        <CursorPos X="38" Y="25"/>
-        <UsageCount Value="9"/>
+        <TopLine Value="79"/>
+        <CursorPos X="29" Y="96"/>
+        <UsageCount Value="10"/>
       </Unit17>
       <Unit18>
         <Filename Value="../units/fftw.pas"/>
         <EditorIndex Value="-1"/>
-        <TopLine Value="76"/>
-        <UsageCount Value="9"/>
+        <TopLine Value="67"/>
+        <CursorPos X="60" Y="95"/>
+        <UsageCount Value="10"/>
       </Unit18>
     </Units>
     <JumpHistory Count="30" HistoryIndex="29">
       <Position1>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="643" Column="48" TopLine="623"/>
+        <Caret Line="595" TopLine="576"/>
       </Position1>
       <Position2>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="108" Column="77" TopLine="94"/>
+        <Caret Line="597" TopLine="576"/>
       </Position2>
       <Position3>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="645" TopLine="622"/>
+        <Caret Line="604" TopLine="585"/>
       </Position3>
       <Position4>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="646" Column="59" TopLine="626"/>
+        <Caret Line="605" TopLine="585"/>
       </Position4>
       <Position5>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="824" Column="25" TopLine="798"/>
+        <Caret Line="606" TopLine="585"/>
       </Position5>
       <Position6>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="843" Column="45" TopLine="818"/>
+        <Caret Line="608" TopLine="585"/>
       </Position6>
       <Position7>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="845" Column="13" TopLine="819"/>
+        <Caret Line="609" TopLine="585"/>
       </Position7>
       <Position8>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="650" TopLine="632"/>
+        <Caret Line="615" TopLine="600"/>
       </Position8>
       <Position9>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="640" TopLine="632"/>
+        <Caret Line="616" TopLine="600"/>
       </Position9>
       <Position10>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="642" TopLine="632"/>
+        <Caret Line="617" TopLine="600"/>
       </Position10>
       <Position11>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="644" TopLine="632"/>
+        <Caret Line="618" TopLine="600"/>
       </Position11>
       <Position12>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="645" TopLine="632"/>
+        <Caret Line="565" Column="39" TopLine="548"/>
       </Position12>
       <Position13>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="646" TopLine="632"/>
+        <Caret Line="16" Column="100"/>
       </Position13>
       <Position14>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="648" TopLine="632"/>
+        <Caret Line="108" Column="73" TopLine="88"/>
       </Position14>
       <Position15>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="647" TopLine="632"/>
+        <Caret Line="109" Column="40" TopLine="89"/>
       </Position15>
       <Position16>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="648" TopLine="632"/>
+        <Caret Line="642" Column="57" TopLine="622"/>
       </Position16>
       <Position17>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="650" TopLine="632"/>
+        <Caret Line="650" Column="57" TopLine="630"/>
       </Position17>
       <Position18>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="640" TopLine="632"/>
+        <Caret Line="112" Column="49" TopLine="83"/>
       </Position18>
       <Position19>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="642" TopLine="632"/>
       </Position19>
       <Position20>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="644" TopLine="632"/>
+        <Caret Line="88" Column="62" TopLine="83"/>
       </Position20>
       <Position21>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="645" TopLine="632"/>
+        <Caret Line="936" Column="69" TopLine="486"/>
       </Position21>
       <Position22>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="646" TopLine="632"/>
       </Position22>
       <Position23>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="648" TopLine="632"/>
+        <Caret Line="108" Column="9" TopLine="93"/>
       </Position23>
       <Position24>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="647" TopLine="632"/>
+        <Caret Line="109" Column="95" TopLine="93"/>
       </Position24>
       <Position25>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="648" TopLine="632"/>
+        <Caret Line="562" TopLine="551"/>
       </Position25>
       <Position26>
-        <Filename Value="Plasmapropagation.lpr"/>
-        <Caret Line="47" TopLine="14"/>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="644" Column="28" TopLine="627"/>
       </Position26>
       <Position27>
-        <Filename Value="Plasmapropagation.lpr"/>
-        <Caret Line="63" Column="32" TopLine="38"/>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="672" Column="3" TopLine="548"/>
       </Position27>
       <Position28>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="650" Column="47" TopLine="611"/>
+        <Caret Line="715" Column="2" TopLine="463"/>
       </Position28>
       <Position29>
         <Filename Value="Physikunit.pas"/>
-        <Caret Line="81" TopLine="81"/>
+        <Caret Line="843" Column="46" TopLine="811"/>
       </Position29>
       <Position30>
-        <Filename Value="Plasmapropagation.lpr"/>
-        <Caret Line="37" Column="40" TopLine="16"/>
+        <Filename Value="Physikunit.pas"/>
+        <Caret Line="834" Column="43" TopLine="819"/>
       </Position30>
     </JumpHistory>
   </ProjectSession>
+  <Debugging>
+    <BreakPoints Count="1">
+      <Item1>
+        <Kind Value="bpkSource"/>
+        <WatchScope Value="wpsLocal"/>
+        <WatchKind Value="wpkWrite"/>
+        <Source Value="Plasmapropagation.lpr"/>
+        <Line Value="46"/>
+      </Item1>
+    </BreakPoints>
+  </Debugging>
 </CONFIG>