summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--Physikunit.pas93
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;