diff options
| -rw-r--r-- | Physikunit.pas | 60 |
1 files changed, 59 insertions, 1 deletions
diff --git a/Physikunit.pas b/Physikunit.pas index 1327173..5beac01 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -9,6 +9,7 @@ unit Physikunit; {$DEFINE Zeitschrittueberwachung} {$DEFINE Dichteueberwachung} {$DEFINE negativeDichteueberwachung} +{ $DEFINE exzessiveArrayBereichsTests} interface @@ -109,7 +110,9 @@ type ffts: array of array[0..1,boolean] of fftw_plan_double; // fft-Algorithmen fftTmp,lFftTmp: Pcomplex_double; // Zwischenspeicher für ffts procedure initialisiereDichte(ort,tlc: longint; breite,n: double); inline; + {$IFDEF exzessiveArrayBereichsTests} procedure pruefeArrayEnden(fehler: string); + {$ENDIF} procedure nichtnegativieren; procedure berechnePhasenraumAbleitungen; inline; public @@ -667,7 +670,9 @@ begin fftw_plan_many_dft_c2r(1,@_aP,aX,fftTmp,nil,1,aP div 2 + 1,impulsRaumGradient[i,1],nil,1,aP,[fftw_destroy_input,fftw_exhaustive]); end; + {$IFDEF exzessiveArrayBereichsTests} pruefeArrayEnden('FFT schreibt beim Planen über das Ende hinaus!'); + {$ENDIF} (* fftw_plan_many_dft_r2c(int rank, const int *n, int howmany, @@ -762,6 +767,7 @@ begin (impulsraum[tlc,false]+i+ort*aP)^:=(impulsraum[tlc,false]+i+ort*aP)^ * ges; end; +{$IFDEF exzessiveArrayBereichsTests} procedure tFelder.pruefeArrayEnden(fehler: string); var i,j: longint; @@ -787,6 +793,7 @@ begin raise exception.create('Fehler in tFelder.pruefeArrayEnden!'); end; end; +{$ENDIF} procedure tFelder.berechneAbleitungen; var @@ -1132,11 +1139,27 @@ var i,j,k: longint; tmp: double; begin + {$IFDEF exzessiveArrayBereichsTests} pruefeArrayEnden('Jemand hat das FFT-Array überschrieben!'); + {$ENDIF} for i:=0 to length(ffts)-1 do begin + {$IFDEF exzessiveArrayBereichsTests} + for j:=0 to aP*(aX div 2 + 1) - 1 do begin + (fftTmp+j)^.re:=pi; + (fftTmp+j)^.im:=-pi; + end; + {$ENDIF} + fftw_execute(ffts[i,0,false]); // FFT in x + + {$IFDEF exzessiveArrayBereichsTests} + for j:=0 to aP*(aX div 2 + 1) - 1 do + if ((fftTmp+j)^.re=pi) or ((fftTmp+j)^.im=-pi) then + raise Exception.create('Vorwärts-FFT ('+inttostr(i)+',0) schreibt nicht auf Position '+inttostr(j)+'!'); + pruefeArrayEnden('Vorwärts-FFT ('+inttostr(i)+',0) schreibt beim Ausführen über das Ende hinaus!'); + {$ENDIF} for k:=0 to aP-1 do for j:=0 to aX div 2 do begin // *-i*k*2*pi; @@ -1145,11 +1168,35 @@ begin (fftTmp+j+k*(aX div 2 + 1))^.im:=- 2*pi/dX/aX * j * tmp/aX * byte(3*j<aX); end; + {$IFDEF exzessiveArrayBereichsTests} + for j:=0 to aP*aX - 1 do + (impulsRaumGradient[i,0]+j)^:=pi+1; + {$ENDIF} + fftw_execute(ffts[i,0,true]); + + {$IFDEF exzessiveArrayBereichsTests} + for j:=0 to aP*(aX div 2 + 1) - 1 do + if (impulsRaumGradient[i,0]+j)^=pi+1 then + raise Exception.create('Rückwärts-FFT ('+inttostr(i)+',0) schreibt nicht auf Position '+inttostr(j)+'!'); + pruefeArrayEnden('Rückwärts-FFT ('+inttostr(i)+',0) schreibt beim Ausführen über das Ende hinaus!'); + for j:=0 to aX*(aP div 2 + 1) - 1 do begin + (fftTmp+j)^.re:=pi; + (fftTmp+j)^.im:=-pi; + end; + {$ENDIF} + fftw_execute(ffts[i,1,false]); // FFT in p + + {$IFDEF exzessiveArrayBereichsTests} + for j:=0 to aX*(aP div 2 + 1) - 1 do + if ((fftTmp+j)^.re=pi) or ((fftTmp+j)^.im=-pi) then + raise Exception.create('Vorwärts-FFT ('+inttostr(i)+',0) schreibt nicht auf Position '+inttostr(j)+'!'); + pruefeArrayEnden('Vorwärts-FFT ('+inttostr(i)+',1) schreibt beim Ausführen über das Ende hinaus!'); + {$ENDIF} for k:=0 to aP div 2 do for j:=0 to aX-1 do begin // *-i*k*2*pi; @@ -1158,9 +1205,20 @@ begin (fftTmp+j+k*aX)^.im:=- 2*pi/dP/aP * k * tmp/aP * byte(4*k < aP); end; + {$IFDEF exzessiveArrayBereichsTests} + for j:=0 to aP*aX - 1 do + (impulsRaumGradient[i,1]+j)^:=pi+1; + {$ENDIF} + fftw_execute(ffts[i,1,true]); - pruefeArrayEnden('Rückwärts-FFT ('+inttostr(i)+',1) schreibt beim Ausführen über das Ende hinaus!'); + {$IFDEF exzessiveArrayBereichsTests} + for j:=0 to aP*(aX div 2 + 1) - 1 do + if (impulsRaumGradient[i,1]+j)^=pi+1 then + raise Exception.create('Rückwärts-FFT ('+inttostr(i)+',1) schreibt nicht auf Position '+inttostr(j)+'!'); + + pruefeArrayEnden('Rückwärts-FFT ('+inttostr(i)+',1) schreibt beim Ausführen über das Ende hinaus!'); + {$ENDIF} end; end; |