diff options
| -rw-r--r-- | Physikunit.pas | 93 |
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; |