summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorErich Eckner <git@eckner.net>2015-07-30 15:15:24 +0200
committerErich Eckner <git@eckner.net>2015-07-30 15:15:24 +0200
commitefb7499d2a7e6cb0e384431ae773371606851a77 (patch)
tree76af6833dbf77a44e43597ed497bedfa4af126ad
parenta2fe917e7cf1a6e9eacf92dccabe4114fd145332 (diff)
downloadPlasmapropagation-efb7499d2a7e6cb0e384431ae773371606851a77.tar.xz
Daten und Programm aus git entfernt, RK3/8 eingefuegt
-rw-r--r--.gitignore1
-rw-r--r--AXtestbin15079072 -> 0 bytes
-rw-r--r--AYtestbin15079072 -> 0 bytes
-rw-r--r--Gammatestbin15079072 -> 0 bytes
-rw-r--r--Ntestbin15079072 -> 0 bytes
-rw-r--r--Physikunit.pas520
-rwxr-xr-xPlasmapropagationbin3316560 -> 0 bytes
-rw-r--r--Plasmapropagation.lps113
-rw-r--r--dAYDTtestbin15079072 -> 0 bytes
-rw-r--r--dPhiDXtestbin15079072 -> 0 bytes
-rw-r--r--dPsiDXtestbin15079072 -> 0 bytes
-rw-r--r--input.plap6
12 files changed, 437 insertions, 203 deletions
diff --git a/.gitignore b/.gitignore
index 17147b4..69a5ac3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,6 +7,7 @@ error
*.o
*.zip
*.tar.gz
+Plasmapropagation
Ergebnisse
lib
epost
diff --git a/AXtest b/AXtest
deleted file mode 100644
index a67e5fc..0000000
--- a/AXtest
+++ /dev/null
Binary files differ
diff --git a/AYtest b/AYtest
deleted file mode 100644
index 5c8f0b7..0000000
--- a/AYtest
+++ /dev/null
Binary files differ
diff --git a/Gammatest b/Gammatest
deleted file mode 100644
index eb6cdf5..0000000
--- a/Gammatest
+++ /dev/null
Binary files differ
diff --git a/Ntest b/Ntest
deleted file mode 100644
index 534a0a7..0000000
--- a/Ntest
+++ /dev/null
Binary files differ
diff --git a/Physikunit.pas b/Physikunit.pas
index 8c216c9..da9d56e 100644
--- a/Physikunit.pas
+++ b/Physikunit.pas
@@ -2,16 +2,20 @@ unit Physikunit;
{$mode objfpc}{$H+}
+{$IFNDEF UNIX}
+ {$ERROR This program can be compiled only on/for Unix/Linux based systems.}
+{$ENDIF}
+
{$DEFINE Zeitschrittueberwachung}
{$DEFINE Dichteueberwachung}
interface
uses
- Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit;
+ Classes, SysUtils, Math, protokollunit, matheunit, mystringlistunit, lowlevelunit, baseUnix;
type
- tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaVier);
+ tZeitverfahren = (zfEulerVorwaerts,zfRungeKuttaDreiAchtel,zfRungeKuttaVier);
tVerteilungsfunktion = function(x: extended): extended;
tEMFeldInhalt = (
efA,efAX,efAY,efAZ,
@@ -54,7 +58,7 @@ type
ortsGrenzen: array of extended;
dichteStuecke: array of string;
public
- nMax, // maximale Massendichte in nc
+ nMin,nMax, // minimale und maximale Massendichte in nc
spezifischeLadung: extended; // q/m in e/me
constructor create;
destructor destroy; override;
@@ -66,9 +70,9 @@ type
{ tWertePunkt }
tWertePunkt = class(tObject) // repräsentiert alle Werte eines Punktes im Gitter und deren Zeitableitungen
- matWerte: array of array[tMaterieFeldInhalt,boolean] of double;
+ matWerte: array of array[tMaterieFeldInhalt,boolean] of extended;
// Materiefelder (getrennt für verschiedene Teilchenspezies) und deren Zeitableitungen
- emWerte: array[tEMFeldInhalt,boolean] of double;
+ emWerte: array[tEMFeldInhalt,boolean] of extended;
// EM-Felder und deren Zeitableitungen
// A, p[xyz]? und i?Gamma haben keine sinnvolle Ableitung hier!
lN,rN: tWertePunkt;
@@ -76,10 +80,13 @@ type
constructor create(linkerNachbar: tWertePunkt; materien: longint);
destructor destroy; override;
procedure berechneGammaUndP;
- procedure berechneNAbleitung(iDX,pDNMax: extended);
+ procedure berechneNAbleitung(iDX,pDNMax: extended); overload;
+ procedure berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean); overload;
procedure berechneAbleitungen(dX,iDX: extended);
- procedure liKo(in1,in2: tWertePunkt; fak2: extended); overload; // Werte werden auf (in1 + fak2*in2') gesetzt
- procedure liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+ procedure liKo(in1,in2: tWertePunkt; fak2: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+ procedure liKo(in1,in2,in3: tWertePunkt; fak2,fak3: extended); overload;
+ procedure liKo(in1,in2,in3,in4: tWertePunkt; fak2,fak3,fak4: extended); overload;
+ procedure liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); overload;
function maxDT: extended;
end;
@@ -87,8 +94,9 @@ type
tFelder = class(tObject) // repräsentiert eine Simulationsbox von Feldern inklusive deren Ableitungen
private
- matAnz: longint;
- lichters: array[boolean] of tMyStringlist;
+ matAnz: longint;
+ spezLadungen: tExtendedArray;
+ lichters: array[boolean] of tMyStringlist;
procedure setzeRaender(dT,iDX,iDT: extended);
public
inhalt: array of tWertePunkt;
@@ -97,8 +105,10 @@ type
constructor create(groesse: longint; teilchen: array of tTeilchenSpezies; lichter: tMyStringList; parent: tGitter);
destructor destroy; override;
procedure berechneAbleitungen(dT,dX,iDT,iDX,pDNMax: extended);
- procedure liKo(in1,in2: tFelder; fak2: extended); overload; // Werte werden auf (in1 + fak2*in2') gesetzt
- procedure liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+ procedure liKo(in1,in2: tFelder; fak2: extended); overload; // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+ procedure liKo(in1,in2,in3: tFelder; fak2,fak3: extended); overload;
+ procedure liKo(in1,in2,in3,in4: tFelder; fak2,fak3,fak4: extended); overload;
+ procedure liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); overload;
function maxDT: extended;
end;
@@ -112,6 +122,7 @@ type
dX,iDX,pDNMax,dTMaximum,xl,t: extended;
kvs: tKnownValues;
procedure diffusionsTermAnpassen(pro: tProtokollant);
+ function pruefeMaxDT(aF: longint; var mDT,dT: extended): boolean; // wurde verkleinert?
public
aktuelleFelder: longint;
@@ -134,11 +145,14 @@ type
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?
end;
+procedure SignalCapture(signal : longint); cdecl;
+
implementation
const
@@ -149,6 +163,9 @@ const
'N','DPSIDX','P','PX','PY','PZ','GAMMA','IGAMMA'
);
+var
+ sighupSimulationen: array of tSimulation;
+
{ tAusgabeDatei }
constructor tAusgabeDatei.create(feldName,prefix,suffix: string; prot: tProtokollant; nam: string);
@@ -284,7 +301,7 @@ begin
for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
while cX>=sX do begin
dec(cnt);
- blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].emWerte[emInhalt,ableitung],sizeof(double));
+ blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].emWerte[emInhalt,ableitung],sizeof(extended));
sX:=sX+sDX;
end;
cX:=cX+gitter.dX;
@@ -294,7 +311,7 @@ begin
for i:=0 to length(gitter.felders[gitter.aktuelleFelder].inhalt)-1 do begin
while cX>=sX do begin
dec(cnt);
- blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].matWerte[teilchen,matInhalt,ableitung],sizeof(double));
+ blockWrite(datei,gitter.felders[gitter.aktuelleFelder].inhalt[i].matWerte[teilchen,matInhalt,ableitung],sizeof(extended));
sX:=sX+sDX;
end;
cX:=cX+gitter.dX;
@@ -317,7 +334,8 @@ begin
setlength(ortsGrenzen,0);
fillchar(dichteStuecke,sizeof(dichteStuecke),#0);
setlength(dichteStuecke,0);
- nMax:=0;
+ nMin:=1E-3;
+ nMax:=nMin;
spezifischeLadung:=0;
end;
@@ -350,7 +368,7 @@ begin
while (i<length(ortsGrenzen)) and (ortsGrenzen[i]<x) do
inc(i);
kvs.add('x',x);
- result:=nMax*exprToFloat(false,dichteStuecke[i],kvs,nil);
+ result:=(nMax-nMin)*exprToFloat(false,dichteStuecke[i],kvs,nil)+nMin;
kvs.rem('x');
end;
@@ -479,15 +497,35 @@ begin
for i:=0 to length(matWerte)-1 do
// dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n
matWerte[i,mfN,true]:=
- ( rN.matWerte[i,mfN,false] * rN.matWerte[i,mfIGamma,false] * (pDNMax - rN.matWerte[i,mfPX,false])
- - ln.matWerte[i,mfN,false] * ln.matWerte[i,mfIGamma,false] * (-pDNMax - rN.matWerte[i,mfPX,false])
- - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * 2 * pDNMax) * iDX/2;
+ ( rN.matWerte[i,mfN,false] * rN.matWerte[i,mfIGamma,false] * (pDNMax - rN.matWerte[i,mfPX,false])
+ + lN.matWerte[i,mfN,false] * lN.matWerte[i,mfIGamma,false] * (pDNMax + lN.matWerte[i,mfPX,false])
+ - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * 2 * pDNMax) * iDX/2;
// Der zweite Summand entspricht (in etwa) einer thermischen Verschmierung
// und soll die numerische Stabilität bis zu p * dn/dx / n von 2*pDNMax gewährleisten.
// Es handelt sich um einen Diffusionsterm dn/dt = alpha * Laplace n mit
// alpha = pDNMax * dX
end;
+procedure tWertePunkt.berechneNAbleitung(iDX,pDNMax: extended; rechts: boolean);
+var
+ i: longint;
+begin
+ if rechts then begin
+ for i:=0 to length(matWerte)-1 do
+ // rechter Randterm von dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n
+ matWerte[i,mfN,true]:=
+ ( lN.matWerte[i,mfN,false] * lN.matWerte[i,mfIGamma,false] * (pDNMax + lN.matWerte[i,mfPX,false])
+ - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * (pDNMax - matWerte[i,mfPX,false])) * iDX/2;
+ end
+ else begin
+ for i:=0 to length(matWerte)-1 do
+ // linker Randterm von dn/dt = -d(n/Gamma p_x)/dx + (p_max*dx) * Laplace n
+ matWerte[i,mfN,true]:=
+ ( rN.matWerte[i,mfN,false] * rN.matWerte[i,mfIGamma,false] * (pDNMax - rN.matWerte[i,mfPX,false])
+ - matWerte[i,mfN,false] * matWerte[i,mfIGamma,false] * (pDNMax + matWerte[i,mfPX,false])) * iDX/2;
+ end
+end;
+
procedure tWertePunkt.berechneAbleitungen(dX,iDX: extended); // Zeitableitungen berechnen
var
i: longint;
@@ -558,6 +596,68 @@ begin
matWerte[i,maF,false]:= in1.matWerte[i,maF,false] + fak2 * in2.matWerte[i,maF,true];
end;
+procedure tWertePunkt.liKo(in1,in2,in3: tWertePunkt; fak2,fak3: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+var
+ emF: tEMFeldInhalt;
+ maF: tMaterieFeldInhalt;
+ i: longint;
+begin
+(* tEMFeldInhalt = (
+ efA,efAX,efAY,efAZ,
+ efDAXDT,efDAYDT,efDAZDT,
+ efDPhiDX
+ ); *)
+ for emF:=efAX to efDPhiDX do // alles außer efA, welchen Ableitung ja nicht berechnet wurde
+ emWerte[emF,false]:=
+ in1.emWerte[emF,false]
+ + fak2 * in2.emWerte[emF,true]
+ + fak3 * in3.emWerte[emF,true];
+
+(* tMaterieFeldInhalt = (
+ mfN,mfDPsiDX,
+ mfP,mfPX,mfPY,mfPZ,
+ mfGamma,mfIGamma
+ ); *)
+ for i:=0 to length(matWerte)-1 do // siehe oben
+ for maF:=mfN to mfDPsiDX do
+ matWerte[i,maF,false]:=
+ in1.matWerte[i,maF,false]
+ + fak2 * in2.matWerte[i,maF,true]
+ + fak3 * in3.matWerte[i,maF,true];
+end;
+
+procedure tWertePunkt.liKo(in1,in2,in3,in4: tWertePunkt; fak2,fak3,fak4: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+var
+ emF: tEMFeldInhalt;
+ maF: tMaterieFeldInhalt;
+ i: longint;
+begin
+(* tEMFeldInhalt = (
+ efA,efAX,efAY,efAZ,
+ efDAXDT,efDAYDT,efDAZDT,
+ efDPhiDX
+ ); *)
+ for emF:=efAX to efDPhiDX do // alles außer efA, welchen Ableitung ja nicht berechnet wurde
+ emWerte[emF,false]:=
+ in1.emWerte[emF,false]
+ + fak2 * in2.emWerte[emF,true]
+ + fak3 * in3.emWerte[emF,true]
+ + fak4 * in4.emWerte[emF,true];
+
+(* tMaterieFeldInhalt = (
+ mfN,mfDPsiDX,
+ mfP,mfPX,mfPY,mfPZ,
+ mfGamma,mfIGamma
+ ); *)
+ for i:=0 to length(matWerte)-1 do // siehe oben
+ for maF:=mfN to mfDPsiDX do
+ matWerte[i,maF,false]:=
+ in1.matWerte[i,maF,false]
+ + fak2 * in2.matWerte[i,maF,true]
+ + fak3 * in3.matWerte[i,maF,true]
+ + fak4 * in4.matWerte[i,maF,true];
+end;
+
procedure tWertePunkt.liKo(in1,in2,in3,in4,in5: tWertePunkt; fak2,fak3,fak4,fak5: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
var
emF: tEMFeldInhalt;
@@ -596,11 +696,11 @@ function tWertePunkt.maxDT: extended;
var
i: longint;
begin
- result:=1;
+ result:=infinity;
for i:=0 to length(matWerte)-1 do
if (matWerte[i,mfN,true]<0) and
(matWerte[i,mfN,false]>0) then
- result:=min(result,matWerte[i,mfN,false]/matWerte[i,mfN,true]);
+ result:=min(result,-matWerte[i,mfN,false]/matWerte[i,mfN,true]);
end;
{ tFelder }
@@ -613,6 +713,10 @@ begin
inherited create;
par:=parent;
matAnz:=length(teilchen);
+ fillchar(spezLadungen,sizeof(spezLadungen),#0);
+ setlength(spezLadungen,matAnz);
+ for i:=0 to length(spezLadungen)-1 do
+ spezLadungen[i]:=teilchen[i].spezifischeLadung;
fillchar(inhalt,sizeof(inhalt),#0);
setlength(inhalt,groesse+4); // zwei Felder links und rechts extra für Randbedingungen
inhalt[0]:=tWertePunkt.create(nil,matAnz);
@@ -634,7 +738,13 @@ begin
end;
destructor tFelder.destroy;
+var
+ i: longint;
begin
+ setlength(spezLadungen,0);
+ for i:=0 to length(inhalt)-1 do
+ inhalt[i].free;
+ setlength(inhalt,0);
lichters[false].free;
lichters[true].free;
inherited destroy;
@@ -645,7 +755,7 @@ var
emF: tEMFeldInhalt;
rechts: boolean;
i: longint;
- oVal: extended;
+ nVal: extended;
begin
for emF:=efAX to efAZ do begin // Vakuumrandbedingungen für das A-Feld
inhalt[0].emWerte[emF,true]:=
@@ -658,12 +768,12 @@ begin
for rechts:=false to true do
for i:=0 to lichters[rechts].count-1 do begin
- par.kvs.add('t',par.t-dT);
- oVal:=exprToFloat(false,lichters[rechts][i],par.kvs,nil);
+ par.kvs.add('t',par.t+dT);
+ nVal:=exprToFloat(false,lichters[rechts][i],par.kvs,nil);
par.kvs.add('t',par.t);
inhalt[(length(inhalt)-1)*byte(rechts)].emWerte[efAY,true]:=
inhalt[(length(inhalt)-1)*byte(rechts)].emWerte[efAY,true] +
- (exprToFloat(false,lichters[rechts][i],par.kvs,nil) - oVal)*iDT;
+ (nVal - exprToFloat(false,lichters[rechts][i],par.kvs,nil))*iDT;
end;
end;
@@ -680,11 +790,15 @@ begin
-abs(inhalt[length(inhalt)-2].matWerte[i,mfPX,false]);
end;
+ for i:=2 to length(inhalt)-3 do
+ inhalt[i].berechneNAbleitung(iDX,pDNMax);
+
+ inhalt[1].berechneNAbleitung(iDX,pDNMax,false);
+ inhalt[length(inhalt)-2].berechneNAbleitung(iDX,pDNMax,true);
+
setzeRaender(dT,iDX,iDT);
for i:=1 to length(inhalt)-2 do
- inhalt[i].berechneNAbleitung(iDX,pDNMax);
- for i:=1 to length(inhalt)-2 do
inhalt[i].berechneAbleitungen(dX,iDX);
end;
@@ -696,6 +810,22 @@ begin
inhalt[i].liKo(in1.inhalt[i],in2.inhalt[i],fak2);
end;
+procedure tFelder.liKo(in1,in2,in3: tFelder; fak2,fak3: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+var
+ i: longint;
+begin
+ for i:=0 to length(inhalt)-1 do
+ inhalt[i].liKo(in1.inhalt[i],in2.inhalt[i],in3.inhalt[i],fak2,fak3);
+end;
+
+procedure tFelder.liKo(in1,in2,in3,in4: tFelder; fak2,fak3,fak4: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
+var
+ i: longint;
+begin
+ for i:=0 to length(inhalt)-1 do
+ inhalt[i].liKo(in1.inhalt[i],in2.inhalt[i],in3.inhalt[i],in4.inhalt[i],fak2,fak3,fak4);
+end;
+
procedure tFelder.liKo(in1,in2,in3,in4,in5: tFelder; fak2,fak3,fak4,fak5: extended); // Werte werden auf (in1 + \sum_i faki * ini') gesetzt
var
i: longint;
@@ -732,8 +862,10 @@ begin
end
else pDNMax:=pDNMa;
case Zeitverfahren of
- zfEulerVorwaerts: Setlength(Felders,2);
- zfRungeKuttaVier: Setlength(Felders,5);
+ zfEulerVorwaerts:
+ Setlength(Felders,2);
+ zfRungeKuttaDreiAchtel,zfRungeKuttaVier:
+ Setlength(Felders,5);
end{of Case};
xl:=dX/2;
@@ -755,102 +887,178 @@ begin
inherited destroy;
end;
+procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant);
+var
+ i,j: longint;
+begin
+ for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do
+ for j:=0 to felders[aktuelleFelder].matAnz-1 do
+ if abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPx,false]*
+ (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
+ max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10))
+ > pDNMax then begin
+ if pDNMaxDynamisch then begin
+ pDNMax:=
+ 2*
+ abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPX,false]*
+ (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
+ Max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10));
+ pro.schreibe('Neues maximales px * dn/dx / n: '+floattostr(pDNMax)+' (t='+floattostr(t)+')');
+ end
+ else
+ if pDNMax>0 then begin
+ pro.schreibe('Warnung: maximaler Impuls * dn/dx / n von '+floattostr(pDNMax)+' wurde überschritten (t='+floattostr(t)+
+ 'T), die numerische Stabilität ist nicht mehr gewährleistet!',true);
+ pro.schreibe(' Lösung: größeren Diffusionsterm wählen (-D)',true);
+ pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern (-x)',true);
+ pDNMax:=-1;
+ end;
+ end;
+end;
+
+function tGitter.pruefeMaxDT(aF: longint; var mDT,dT: extended): boolean; // wurde verkleinert?
+begin
+ mDT:=min(mDT,felders[aF].maxDT);
+ // das maximale dT, welches bei den Momentanen Ableitungen nicht zu
+ // unphysikalischen Effekten (z.B. negativen Dichten) führt
+
+ result:=mDT<dT*2;
+ if result then begin // beabsichtigter Zeitschritt ist zu groß,
+ dT:=mDT/4; // dann machen wir ihn doch kleiner!
+ {$IFDEF Zeitschrittueberwachung}
+ prot.schreibe('pruefeMaxDT: dT -> '+floattostr(dT));
+ if dT<1E-30 then begin
+ prot.schreibe('pruefeMaxDT: Zeitschritt geht gegen Null (ist bereits '+floattostr(dT)+' < 10^-30) - irgendwas scheint grundsätzlich kaputt zu sein!',true);
+ prot.destroyall;
+ halt(1);
+ end;
+ {$ENDIF}
+ exit;
+ end;
+end;
+
procedure tGitter.iteriereSchritt(var dT: extended);
var i: longint;
{$IFDEF Zeitschrittueberwachung}
- pro: tProtokollant;
j: longint;
{$ENDIF}
+ {$IFDEF Dichteueberwachung}
+ pro: tProtokollant;
+ {$ENDIF}
iDT,mDT: extended;
begin
- kvs.add('t',t+dT);
+ kvs.add('t',t);
diffusionsTermAnpassen(prot);
+ mDT:=infinity;
+
repeat
iDT:=1/dT;
felders[aktuelleFelder].berechneAbleitungen(dT,dX,iDT,iDX,pDNMax); // y' = y'(t,y(t))
- mDT:=felders[aktuelleFelder].maxDT;
- // das maximale dT, welches bei den Momentanen Ableitungen nicht zu
- // unphysikalischen Effekten (z.B. negativen Dichten) führt
-
- if mDT<dT*2 then begin // beabsichtigter Zeitschritt ist zu groß
- dT:=mDT/4; // machen wir ihn kleiner!
- {$IFDEF Zeitschrittueberwachung}
- pro:=tProtokollant.create(prot,'iteriereSchritt');
- pro.schreibe('dT -> '+floattostr(dT));
- if dT<1E-30 then begin
- pro.schreibe('Zeitschritt geht gegen Null (ist bereits '+floattostr(dT)+' < 10^-30) - irgendwas scheint grundsätzlich kaputt zu sein!',true);
- pro.destroyall;
- halt(1);
- end;
- pro.free;
- {$ENDIF}
+ if pruefeMaxDT(aktuelleFelder,mDT,dT) then
continue;
- end;
+
+ case zeitverfahren of
+ zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt
+ felders[1-aktuelleFelder].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT);
+ zfRungeKuttaDreiAchtel:
+ begin
+ felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT/3); // ya = y(t) + y' dt/3
+ felders[2].berechneAbleitungen(dT/3,dX,iDT,iDX,pDNMax); // ya' = y'(t+dt/3,ya)
+
+ if pruefeMaxDT(2,mDT,dT) then
+ continue;
+
+ felders[3].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[2],-dT/3,dT); // yb = y(t) - y' dt/3 + ya' dt
+ felders[3].berechneAbleitungen(dT/3,dX,iDT,iDX,pDNMax); // yb' = y'(t+2dt/3,yb)
+
+ if pruefeMaxDT(3,mDT,dT) then
+ continue;
+
+ felders[4].liKo(felders[aktuelleFelder],felders[aktuelleFelder],felders[2],felders[3],dT,-dT,dT); // yc = y(t) + a' dt - ya' dt + yb' dt
+ felders[4].berechneAbleitungen(dT/3,dX,iDT,iDX,pDNMax); // yc' = y'(t+dt,yc)
+
+ if pruefeMaxDT(4,mDT,dT) then
+ continue;
+
+ felders[1-aktuelleFelder].liKo(
+ felders[aktuelleFelder],
+ felders[aktuelleFelder],
+ felders[2],
+ felders[3],
+ felders[4],
+ dT/8,
+ dT*3/8,
+ dT*3/8,
+ dT/8
+ ); // y(t+dt) = y(t) + (y' + 3(ya' + yb') + yc') dt/8
+ end;
+ zfRungeKuttaVier:
+ begin
+ felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT/2); // ya = y(t) + y' dt/2
+ felders[2].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // ya' = y'(t+dt/2,ya)
+
+ if pruefeMaxDT(2,mDT,dT) then
+ continue;
+
+ felders[3].liKo(felders[aktuelleFelder],felders[2],dT/2); // yb = y(t) + ya' dt/2
+ felders[3].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // yb' = y'(t+dt/2,yb)
+
+ if pruefeMaxDT(3,mDT,dT) then
+ continue;
+
+ felders[4].liKo(felders[aktuelleFelder],felders[3],dT); // yc = y(t) + yb' dt
+ felders[4].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // yc' = y'(t+dt,yc)
+
+ if pruefeMaxDT(4,mDT,dT) then
+ continue;
+
+ felders[1-aktuelleFelder].liKo(
+ felders[aktuelleFelder],
+ felders[aktuelleFelder],
+ felders[2],
+ felders[3],
+ felders[4],
+ dT/6,
+ dT/3,
+ dT/3,
+ dT/6
+ ); // y(t+dt) = y(t) + (y' + 2(ya' + yb') + yc') dt/6
+ end;
+ end{of case};
break;
until false;
- if dT<mDT/1000 then
- dT:=min(1,mDT)/1000;
-
- case zeitverfahren of
- zfEulerVorwaerts: // y(t+dt) = y(t) + y' dt
- felders[1-aktuelleFelder].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT);
- zfRungeKuttaVier:
- begin
- felders[2].liKo(felders[aktuelleFelder],felders[aktuelleFelder],dT/2); // ya = y(t) + y' dt/2
- felders[2].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // ya' = y'(t+dt/2,ya)
-
- felders[3].liKo(felders[aktuelleFelder],felders[2],dT/2); // yb = y(t) + ya' dt/2
- felders[3].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // yb' = y'(t+dt/2,yb)
-
- felders[4].liKo(felders[aktuelleFelder],felders[3],dT); // yc = y(t) + yb' dt
- felders[4].berechneAbleitungen(dT/2,dX,iDT,iDX,pDNMax); // yc' = y'(t+dt,yc)
-
- felders[1-aktuelleFelder].liKo(
- felders[aktuelleFelder],
- felders[aktuelleFelder],
- felders[2],
- felders[3],
- felders[4],
- dT/6,
- dT/3,
- dT/3,
- dT/6
- ); // y(t+dt) = y(t) + (y' + 2(ya' + yb') + yc') dt/6
- end;
- end{of case};
-
aktuelleFelder:=1-aktuelleFelder;
{$IFDEF Dichteueberwachung}
- with felders[aktuelleFelder] do
- for i:=0 to length(inhalt)-1 do
- for j:=0 to matAnz-1 do
- if inhalt[i].matWerte[j,mfN,false]<0 then begin
- pro:=tProtokollant.create(prot,'iteriereSchritt');
- pro.schreibe(
- 'n<=0 bei '+
- floattostr(t)+' '+
- inttostr(i)+' '+
- floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]*
- inhalt[i-1].matWerte[j,mfIGamma,false])+' '+
- floattostr(inhalt[i-1].matWerte[j,mfN,false])+' '+
- floattostr(inhalt[i-1].matWerte[j,mfN,true])+' '+
- floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]*
- inhalt[i].matWerte[j,mfIGamma,false])+' '+
- floattostr(inhalt[i].matWerte[j,mfN,false])+' '+
- floattostr(inhalt[i].matWerte[j,mfN,true])+' '+
- floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]*
- inhalt[i+1].matWerte[j,mfIGamma,false])+' '+
- floattostr(inhalt[i+1].matWerte[j,mfN,false])+' '+
- floattostr(inhalt[i+1].matWerte[j,mfN,true]),true);
- pro.destroyall;
- halt(1);
- end;
+ for i:=0 to length(felders[aktuelleFelder].inhalt)-1 do
+ for j:=0 to felders[aktuelleFelder].matAnz-1 do
+ if felders[aktuelleFelder].inhalt[i].matWerte[j,mfN,false]<0 then
+ with felders[1-aktuelleFelder] do begin
+ pro:=tProtokollant.create(prot,'iteriereSchritt');
+ pro.schreibe('n<0 bei:',true);
+ pro.schreibe(' t = ' +floattostr(t)+',',true);
+ pro.schreibe(' i = ' +inttostr(i)+' (x = '+floattostr(xl+i*dX)+'),',true);
+ pro.schreibe(' v_l = ' +floattostr(inhalt[i-1].matWerte[j,mfDPsiDX,false]*
+ inhalt[i-1].matWerte[j,mfIGamma,false])+',',true);
+ pro.schreibe(' n_l = ' +floattostr(inhalt[i-1].matWerte[j,mfN,false])+',',true);
+ pro.schreibe(' n_l'' = '+floattostr(inhalt[i-1].matWerte[j,mfN,true])+',',true);
+ pro.schreibe(' v = ' +floattostr(inhalt[i].matWerte[j,mfDPsiDX,false]*
+ inhalt[i].matWerte[j,mfIGamma,false])+',',true);
+ pro.schreibe(' n = ' +floattostr(inhalt[i].matWerte[j,mfN,false])+',',true);
+ pro.schreibe(' n'' = '+floattostr(inhalt[i].matWerte[j,mfN,true])+',',true);
+ pro.schreibe(' v_r = ' +floattostr(inhalt[i+1].matWerte[j,mfDPsiDX,false]*
+ inhalt[i+1].matWerte[j,mfIGamma,false])+',',true);
+ pro.schreibe(' n_r = ' +floattostr(inhalt[i+1].matWerte[j,mfN,false])+',',true);
+ pro.schreibe(' n_r'' = '+floattostr(inhalt[i+1].matWerte[j,mfN,true]),true);
+ pro.destroyall;
+ halt(1);
+ end;
{$ENDIF}
t:=t+dT;
+ dT:=max(dT,min(1,mDT)/100);
end;
function tGitter.dumpErhaltungsgroessen: boolean;
@@ -885,35 +1093,6 @@ begin
pro.free;
end;
-procedure tGitter.diffusionsTermAnpassen(pro: tProtokollant);
-var
- i,j: longint;
-begin
- for i:=1 to length(felders[aktuelleFelder].inhalt)-2 do
- for j:=0 to felders[aktuelleFelder].matAnz-1 do
- if abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPx,false]*
- (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
- max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10))
- > pDNMax then begin
- if pDNMaxDynamisch then begin
- pDNMax:=
- 2*
- abs(felders[aktuelleFelder].inhalt[i].matWerte[j,mfPX,false]*
- (felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] - felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false])/
- Max(felders[aktuelleFelder].inhalt[i+1].matWerte[j,mfN,false] + felders[aktuelleFelder].inhalt[i-1].matWerte[j,mfN,false],1e-10));
- pro.schreibe('Neues maximales px * dn/dx / n: '+floattostr(pDNMax)+' (t='+floattostr(t)+')');
- end
- else
- if pDNMax>0 then begin
- pro.schreibe('Warnung: maximaler Impuls * dn/dx / n von '+floattostr(pDNMax)+' wurde überschritten (t='+floattostr(t)+
- 'T), die numerische Stabilität ist nicht mehr gewährleistet!',true);
- pro.schreibe(' Lösung: größeren Diffusionsterm wählen (-D)',true);
- pro.schreibe(' außerdem empfohlen: Ortsauflösung in gleichem Maße verbessern (-x)',true);
- pDNMax:=-1;
- end;
- end;
-end;
-
procedure tGitter.berechne(was: char; teilchen: longint);
var
i: longint;
@@ -957,6 +1136,7 @@ var
teilchen: array of tTeilchenSpezies;
lichter: tMyStringlist;
pro: tProtokollant;
+ na: pSigActionRec;
begin
inherited create;
prot:=tProtokollant.create(protokollant,name);
@@ -986,6 +1166,7 @@ begin
breite:=10.0;
pDNMax:=-1;
fortschrittsAnzeige:=false;
+ gotSighup:=false;
// Standardeinstellungen Bereich 'ausgaben'
aPrefix:=extractfilepath(paramstr(0));
@@ -1016,6 +1197,10 @@ begin
halt(1);
end;
if s='allgemeinEnde' then break;
+ if s='runge-Kutta-3/8' then begin
+ Zeitverfahren:=zfRungeKuttaDreiAchtel;
+ continue;
+ end;
if s='runge-Kutta-4' then begin
Zeitverfahren:=zfRungeKuttaVier;
continue;
@@ -1140,6 +1325,11 @@ begin
kvs.add('nmax'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].nMax);
continue;
end;
+ if startetMit('minimaldichte ',s) then begin
+ teilchen[length(teilchen)-1].nMin:=exprToFloat(false,s,kvs,nil);
+ kvs.add('nmin'+inttostr(length(teilchen)),teilchen[length(teilchen)-1].nMin);
+ continue;
+ end;
if startetMit('verteilung ',s) then begin
teilchen[length(teilchen)-1].liesDichteFunktion(s,ifile,kvs,teilchen,prot);
continue;
@@ -1178,6 +1368,8 @@ begin
case zeitverfahren of
zfEulerVorwaerts:
pro.schreibe('Iteration mittels Euler-Vorwärts');
+ zfRungeKuttaDreiAchtel:
+ pro.schreibe('Iteration mittels Runge-Kutta-3/8');
zfRungeKuttaVier:
pro.schreibe('Iteration mittels Runge-Kutta-4');
else
@@ -1208,6 +1400,23 @@ begin
end;
pro.free;
+ if length(sighupSimulationen)=0 then begin
+ new(na);
+ na^.sa_Handler := sigActionHandler(@signalCapture);
+ fillchar(na^.sa_Mask,sizeof(na^.sa_mask),#0);
+ na^.sa_Flags := 0;
+ {$ifdef Linux} // Linux specific
+ na^.sa_Restorer := Nil;
+ {$endif}
+ if (fPSigaction(SIGHUP, na, nil) <> 0) then begin
+ writeln(stderr, 'Error: ', fpgeterrno, '.');
+ halt(1);
+ end;
+ dispose(na);
+ end;
+ setlength(sighupSimulationen,length(sighupSimulationen)+1);
+ sighupSimulationen[length(sighupSimulationen)-1]:=self;
+
gitter:=tGitter.create(round(breite/deltaX)+1,dT,deltaX,pDNMax,kvs,teilchen,lichter,zeitverfahren,prot,'gitter');
for i:=0 to length(teilchen)-1 do
teilchen[i].free;
@@ -1216,9 +1425,17 @@ begin
end;
destructor tSimulation.destroy;
+var
+ i,j: longint;
begin
gitter.free;
prot.free;
+ for i:=length(sighupSimulationen)-1 downto 0 do
+ if sighupSimulationen[i]=self then begin
+ for j:=i+1 to length(sighupSimulationen)-1 do
+ sighupSimulationen[j-1]:=sighupSimulationen[j];
+ setlength(sighupSimulationen,length(sighupSimulationen)-1);
+ end;
inherited destroy;
end;
@@ -1240,22 +1457,37 @@ begin
end;
zeitDatei:=zeitDatei+now;
- if fortschrittsAnzeige then begin
- if (floor(100*Gitter.t/tEnde) < floor(100*(Gitter.t+dT)/tEnde)) then begin
- prot.schreibe(inttostr(round(100*Gitter.t/tEnde))+'% (t='+floattostr(Gitter.t)+'T)',true);
- prot.schreibe(timetostr(now-start)+' ('+floattostr(zeitPhysik/max(1e-11,zeitPhysik+zeitDatei))+')',true);
- prot.schreibe('ETA: '+timetostr((now-start)*(tEnde-Gitter.t)/max(Gitter.t,dT)),true);
- prot.schreibe('aktueller Zeitschritt: '+floattostr(dT)+'T',true);
- if errorCode<2 then
- if not gitter.dumpErhaltungsgroessen then begin
- errorCode:=2;
- exit;
- end;
- end;
+ if gotSighup or
+ (fortschrittsAnzeige and (floor(100*Gitter.t/tEnde) < floor(100*(Gitter.t+dT)/tEnde))) then begin
+ gotSighup:=false;
+ prot.schreibe(inttostr(round(100*Gitter.t/tEnde))+'% (t='+floattostr(Gitter.t)+'T)',true);
+ prot.schreibe(timetostr(now-start)+' ('+floattostr(zeitPhysik/max(1e-11,zeitPhysik+zeitDatei))+')',true);
+ prot.schreibe('ETA: '+timetostr((now-start)*(tEnde-Gitter.t)/max(Gitter.t,dT)),true);
+ prot.schreibe('aktueller Zeitschritt: '+floattostr(dT)+'T',true);
+ if errorCode<2 then
+ if not gitter.dumpErhaltungsgroessen then begin
+ errorCode:=2;
+ exit;
+ end;
end;
result:=gitter.t<tEnde;
end;
+// allgemeine Funktionen *******************************************************
+
+procedure SignalCapture(signal : longint); cdecl;
+var
+ i: longint;
+begin
+ if signal=SIGHUP then
+ for i:=0 to length(sighupSimulationen)-1 do
+ sighupSimulationen[i].gotSighup:=true;
+end;
+
+initialization
+ setlength(sighupSimulationen,0);
+finalization
+ setlength(sighupSimulationen,0);
end.
diff --git a/Plasmapropagation b/Plasmapropagation
deleted file mode 100755
index 561167f..0000000
--- a/Plasmapropagation
+++ /dev/null
Binary files differ
diff --git a/Plasmapropagation.lps b/Plasmapropagation.lps
index b5afdb3..e1cfd0c 100644
--- a/Plasmapropagation.lps
+++ b/Plasmapropagation.lps
@@ -7,67 +7,64 @@
<Unit0>
<Filename Value="Plasmapropagation.lpr"/>
<IsPartOfProject Value="True"/>
- <TopLine Value="25"/>
+ <TopLine Value="11"/>
<CursorPos X="34" Y="11"/>
- <UsageCount Value="153"/>
+ <UsageCount Value="164"/>
<Loaded Value="True"/>
</Unit0>
<Unit1>
<Filename Value="Physikunit.pas"/>
<IsPartOfProject Value="True"/>
- <IsVisibleTab Value="True"/>
<EditorIndex Value="1"/>
- <TopLine Value="760"/>
- <CursorPos X="45" Y="781"/>
- <FoldState Value=" T3kO0-4 pigkU095 pjDja083]90jG0M[I4Bk40J513[c4nL0~E2%"/>
- <UsageCount Value="94"/>
+ <TopLine Value="76"/>
+ <CursorPos X="60" Y="110"/>
+ <FoldState Value=" T3kf0-4 pigkU0A5 pjDjb084]9Ija0M3]94jY09[I43jO067[I5S0J1]9akG0U2 poOpJ0D2 T0\1Pc071'"/>
+ <UsageCount Value="105"/>
<Bookmarks Count="1">
- <Item0 Y="790"/>
+ <Item0 Y="1031"/>
</Bookmarks>
<Loaded Value="True"/>
</Unit1>
<Unit2>
<Filename Value="../units/protokollunit.pas"/>
<IsPartOfProject Value="True"/>
- <EditorIndex Value="2"/>
+ <EditorIndex Value="-1"/>
<TopLine Value="20"/>
<CursorPos X="22" Y="31"/>
- <UsageCount Value="56"/>
- <Loaded Value="True"/>
+ <UsageCount Value="67"/>
</Unit2>
<Unit3>
<Filename Value="input.plap"/>
<IsPartOfProject Value="True"/>
- <EditorIndex Value="4"/>
- <TopLine Value="9"/>
- <CursorPos X="30" Y="43"/>
- <UsageCount Value="55"/>
+ <IsVisibleTab Value="True"/>
+ <EditorIndex Value="2"/>
+ <CursorPos Y="5"/>
+ <UsageCount Value="66"/>
<Loaded Value="True"/>
<DefaultSyntaxHighlighter Value="None"/>
</Unit3>
<Unit4>
<Filename Value="input.epost"/>
- <EditorIndex Value="5"/>
+ <EditorIndex Value="3"/>
<CursorPos X="32" Y="21"/>
- <UsageCount Value="60"/>
+ <UsageCount Value="66"/>
<Loaded Value="True"/>
<DefaultSyntaxHighlighter Value="None"/>
</Unit4>
<Unit5>
<Filename Value="../units/matheunit.pas"/>
- <EditorIndex Value="3"/>
- <TopLine Value="83"/>
+ <EditorIndex Value="-1"/>
+ <TopLine Value="53"/>
<CursorPos Y="53"/>
<FoldState Value=" T3i905B pj0jV034 piaj60U2-"/>
- <UsageCount Value="21"/>
- <Loaded Value="True"/>
+ <UsageCount Value="20"/>
</Unit5>
<Unit6>
<Filename Value="../units/lowlevelunit.pas"/>
<EditorIndex Value="-1"/>
- <TopLine Value="540"/>
- <CursorPos X="20" Y="573"/>
- <UsageCount Value="11"/>
+ <TopLine Value="4"/>
+ <CursorPos X="86" Y="23"/>
+ <UsageCount Value="10"/>
</Unit6>
<Unit7>
<Filename Value="../units/mystringlistunit.pas"/>
@@ -75,155 +72,155 @@
<TopLine Value="367"/>
<CursorPos X="17" Y="390"/>
<FoldState Value=" T3i3075 piZjD0WQ"/>
- <UsageCount Value="12"/>
+ <UsageCount Value="11"/>
</Unit7>
<Unit8>
<Filename Value="../epost/werteunit.pas"/>
<EditorIndex Value="-1"/>
<TopLine Value="950"/>
<CursorPos X="30" Y="1054"/>
- <UsageCount Value="9"/>
+ <UsageCount Value="8"/>
</Unit8>
<Unit9>
<Filename Value="../epost/typenunit.pas"/>
<EditorIndex Value="-1"/>
<TopLine Value="347"/>
<CursorPos X="62" Y="358"/>
- <UsageCount Value="9"/>
+ <UsageCount Value="8"/>
</Unit9>
<Unit10>
<Filename Value="../units/systemunit.pas"/>
<EditorIndex Value="-1"/>
<CursorPos X="3" Y="79"/>
- <UsageCount Value="9"/>
+ <UsageCount Value="8"/>
</Unit10>
<Unit11>
<Filename Value="/usr/lib/fpc/src/rtl/inc/objpash.inc"/>
<EditorIndex Value="-1"/>
<TopLine Value="232"/>
<CursorPos X="23" Y="192"/>
- <UsageCount Value="9"/>
+ <UsageCount Value="8"/>
</Unit11>
</Units>
<JumpHistory Count="30" HistoryIndex="29">
<Position1>
<Filename Value="Physikunit.pas"/>
- <Caret Line="112" Column="37" TopLine="92"/>
+ <Caret Line="121" Column="45" TopLine="98"/>
</Position1>
<Position2>
<Filename Value="Physikunit.pas"/>
- <Caret Line="736" TopLine="686"/>
+ <Caret Line="837" Column="13" TopLine="806"/>
</Position2>
<Position3>
<Filename Value="Physikunit.pas"/>
- <Caret Line="1195" Column="41" TopLine="891"/>
+ <Caret Line="843" Column="10" TopLine="806"/>
</Position3>
<Position4>
<Filename Value="Physikunit.pas"/>
- <Caret Line="121" Column="35" TopLine="101"/>
+ <Caret Line="835" Column="63" TopLine="806"/>
</Position4>
<Position5>
<Filename Value="Physikunit.pas"/>
- <Caret Line="806" Column="58" TopLine="678"/>
+ <Caret Line="121" Column="79" TopLine="112"/>
</Position5>
<Position6>
<Filename Value="Physikunit.pas"/>
- <Caret Line="807" Column="39" TopLine="678"/>
+ <Caret Line="861" TopLine="855"/>
</Position6>
<Position7>
<Filename Value="Physikunit.pas"/>
- <Caret Line="832" Column="36" TopLine="678"/>
+ <Caret Line="120" TopLine="120"/>
</Position7>
<Position8>
<Filename Value="Physikunit.pas"/>
- <Caret Line="1195" Column="31" TopLine="862"/>
+ <Caret Line="121" Column="59" TopLine="98"/>
</Position8>
<Position9>
<Filename Value="Physikunit.pas"/>
- <Caret Line="121" Column="35" TopLine="101"/>
+ <Caret Line="851" TopLine="806"/>
</Position9>
<Position10>
<Filename Value="Physikunit.pas"/>
- <Caret Line="813" Column="76" TopLine="678"/>
+ <Caret Line="911" Column="64" TopLine="879"/>
</Position10>
<Position11>
<Filename Value="Physikunit.pas"/>
- <Caret Line="121" Column="37" TopLine="98"/>
+ <Caret Line="913" Column="16" TopLine="897"/>
</Position11>
<Position12>
<Filename Value="Physikunit.pas"/>
- <Caret Line="811" Column="17" TopLine="736"/>
+ <Caret Line="18" Column="77"/>
</Position12>
<Position13>
<Filename Value="Physikunit.pas"/>
- <Caret Line="827" Column="144" TopLine="736"/>
+ <Caret Line="785" Column="44" TopLine="674"/>
</Position13>
<Position14>
<Filename Value="Physikunit.pas"/>
- <Caret Line="1198" Column="20" TopLine="894"/>
+ <Caret Line="915" Column="23" TopLine="883"/>
</Position14>
<Position15>
<Filename Value="Physikunit.pas"/>
+ <Caret Line="1072" Column="34" TopLine="981"/>
</Position15>
<Position16>
<Filename Value="Physikunit.pas"/>
- <Caret Line="827" Column="40" TopLine="735"/>
+ <Caret Line="1119" Column="48" TopLine="1087"/>
</Position16>
<Position17>
<Filename Value="Physikunit.pas"/>
- <Caret Line="1200" Column="5" TopLine="894"/>
+ <Caret Line="1123" Column="42" TopLine="1091"/>
</Position17>
<Position18>
<Filename Value="Physikunit.pas"/>
- <Caret Line="820" Column="52" TopLine="686"/>
+ <Caret Line="1291" Column="6" TopLine="1257"/>
</Position18>
<Position19>
<Filename Value="Physikunit.pas"/>
- <Caret Line="121" Column="45" TopLine="103"/>
+ <Caret Line="18" Column="77"/>
</Position19>
<Position20>
<Filename Value="Physikunit.pas"/>
- <Caret Line="828" Column="13" TopLine="736"/>
</Position20>
<Position21>
<Filename Value="Physikunit.pas"/>
- <Caret Line="807" TopLine="736"/>
+ <Caret Line="18" Column="77"/>
</Position21>
<Position22>
<Filename Value="Physikunit.pas"/>
- <Caret Line="1175" TopLine="782"/>
+ <Caret Line="785" Column="44" TopLine="674"/>
</Position22>
<Position23>
<Filename Value="Physikunit.pas"/>
- <Caret Line="120" Column="33" TopLine="113"/>
+ <Caret Line="915" Column="23" TopLine="883"/>
</Position23>
<Position24>
<Filename Value="Physikunit.pas"/>
- <Caret Line="736" Column="35" TopLine="695"/>
+ <Caret Line="1072" Column="34" TopLine="981"/>
</Position24>
<Position25>
<Filename Value="Physikunit.pas"/>
- <Caret Line="120" Column="35" TopLine="87"/>
+ <Caret Line="1123" Column="42" TopLine="1091"/>
</Position25>
<Position26>
<Filename Value="Physikunit.pas"/>
- <Caret Line="802" Column="51" TopLine="790"/>
+ <Caret Line="1291" Column="21" TopLine="1259"/>
</Position26>
<Position27>
<Filename Value="Physikunit.pas"/>
- <Caret Line="6"/>
+ <Caret Line="891" Column="95" TopLine="872"/>
</Position27>
<Position28>
<Filename Value="Physikunit.pas"/>
- <Caret Line="760" Column="121" TopLine="726"/>
+ <Caret Line="88" Column="64" TopLine="64"/>
</Position28>
<Position29>
<Filename Value="Physikunit.pas"/>
- <Caret Line="736" TopLine="558"/>
+ <Caret Line="597" TopLine="526"/>
</Position29>
<Position30>
<Filename Value="Physikunit.pas"/>
- <Caret Line="701" Column="41" TopLine="465"/>
+ <Caret Line="811" Column="19" TopLine="778"/>
</Position30>
</JumpHistory>
</ProjectSession>
diff --git a/dAYDTtest b/dAYDTtest
deleted file mode 100644
index 5c8f0b7..0000000
--- a/dAYDTtest
+++ /dev/null
Binary files differ
diff --git a/dPhiDXtest b/dPhiDXtest
deleted file mode 100644
index a37f979..0000000
--- a/dPhiDXtest
+++ /dev/null
Binary files differ
diff --git a/dPsiDXtest b/dPsiDXtest
deleted file mode 100644
index 7ba5d78..0000000
--- a/dPsiDXtest
+++ /dev/null
Binary files differ
diff --git a/input.plap b/input.plap
index 5e98ccc..b061987 100644
--- a/input.plap
+++ b/input.plap
@@ -1,7 +1,9 @@
# Parameterdatei für Plasmapropagation
allgemein
- runge-Kutta-4
+ runge-Kutta-3/8
+# runge-Kutta-4
+# euler-Vorwärts
ortsschritt 10^-2 * λ
zeitschritt 10^-2 * T
zeit 20 * T
@@ -24,6 +26,7 @@ licht von links 0.5 * 2^(-2*((t-$tMitte)/$tFwhm)^2) * sin(ω*t)
teilchen1
spezifische Ladung -q/me
maximaldichte 10
+# minimaldichte 1
!setze $profilbreite: (4 * λ)
!setze $randbreite: (0.1 * λ)
verteilung stückweise
@@ -42,6 +45,7 @@ teilchen1Ende
teilchen2
spezifische Ladung (q/me)/2
maximaldichte nmax1*2
+# minimaldichte nmin1*2
verteilung wie teilchen1
teilchen2Ende