diff options
Diffstat (limited to 'Physikunit.pas')
| -rw-r--r-- | Physikunit.pas | 71 |
1 files changed, 55 insertions, 16 deletions
diff --git a/Physikunit.pas b/Physikunit.pas index f95dc32..b4b1456 100644 --- a/Physikunit.pas +++ b/Physikunit.pas @@ -106,8 +106,9 @@ type _aX,_aP: longint; _dX,_dP,iDX: double; ffts: array of array[0..1,boolean] of fftw_plan_double; // fft-Algorithmen - fftTmp: Pcomplex_double; // Zwischenspeicher für ffts + fftTmp,lFftTmp: Pcomplex_double; // Zwischenspeicher für ffts procedure initialisiereDichte(ort,tlc: longint; breite,n: double); inline; + procedure pruefeArrayEnden(fehler: string); public emFelder: array[tEMFeldGroesze,boolean] of pDouble; // EM-Felder und deren Zeitableitungen // dPhiDX und B[xyz] haben keine sinnvolle Ableitung hier! @@ -585,8 +586,9 @@ 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(double)); + fftw_getmem(impulsRaumGradient[i,j],(aP*aX+1)*sizeof(double)); fillchar(impulsRaumGradient[i,j]^,aP*aX*sizeof(double),#0); + (impulsRaumGradient[i,j]+aP*aX)^:=pi; end; fillchar(iMGamma,sizeof(iMGamma),#0); setlength(iMGamma,length(teilchen)); @@ -639,8 +641,13 @@ begin lichters[true].grep('^rechts .*'); lichters[true].subst('^rechts *',''); - fftw_getmem(fftTmp,(aP+1)*(aX+1)*sizeof(complex_double)); - fillchar(fftTmp^,(aP+1)*(aX+1)*sizeof(complex_double),#0); + i:=max((aP div 2 + 1)*aX,(aX div 2 + 1)*aP); + fftw_getmem(fftTmp,(i+1)*sizeof(complex_double)); + fillchar(fftTmp^,(i+1)*sizeof(complex_double),#0); + lFftTmp:=fftTmp+i; + + lFftTmp^.re:=pi; + lFftTmp^.im:=sqr(pi); fillchar(ffts,sizeof(ffts),#0); setlength(ffts,length(teilchen)); @@ -651,11 +658,14 @@ begin 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]); + fftw_plan_many_dft_r2c(1,@_aP,aX,impulsraum[i,false],nil,1,aP,fftTmp,nil,1,aP div 2 + 1,[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]); + fftw_plan_many_dft_c2r(1,@_aP,aX,fftTmp,nil,1,aP div 2 + 1,impulsRaumGradient[i,1],nil,1,aP,[fftw_measure]); end; - (* + + pruefeArrayEnden('FFT schreibt beim Planen über das Ende hinaus!'); + + (* fftw_plan_many_dft_r2c(int rank, const int *n, int howmany, double *in, const int *inembed, int istride, int idist, @@ -729,6 +739,26 @@ begin (impulsraum[tlc,false]+i+ort*aP)^:=(impulsraum[tlc,false]+i+ort*aP)^ * ges; end; +procedure tFelder.pruefeArrayEnden(fehler: string); +var + i,j: longint; +begin + 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))); + gitter.prot.destroyall; + halt(1); + end; + for i:=0 to length(impulsRaumGradient)-1 do + for j:=0 to 1 do + if abs((impulsRaumGradient[i,j]+aP*aX)^-pi)>1e-10 then begin + gitter.prot.schreibe(fehler); + gitter.prot.schreibe(floattostr((impulsRaumGradient[i,j]+aP*aX)^)+' vs. '+floattostr(pi)); + gitter.prot.destroyall; + halt(1); + end; +end; + procedure tFelder.berechneAbleitungen; var i,j,k: longint; @@ -829,26 +859,35 @@ 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]); - for j:=0 to aX-1 do + 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 j:=0 to aX div 2 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; + 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(2*j<ax); + (fftTmp+j+k*(aX div 2 + 1))^.im:= 2*pi/dX/aX * j * tmp/aX * byte(2*j<ax); end; fftw_execute(ffts[i,0,true]); - fftw_execute(ffts[i,1,false]); + 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 j:=0 to aX-1 do - for k:=0 to aP-1 do begin // *i*k*2*pi; + for k:=0 to aP div 2 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; + (fftTmp+j+k*aX)^.re:=- 2*pi/dP/aP * k * (fftTmp+j+k*aX)^.im/aP * byte(2*k < aP); + (fftTmp+j+k*aX)^.im:= 2*pi/dP/aP * k * tmp/aP * byte(2*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; // Zeitableitung der Phasenraumbesetzungsdichte berechnen |