diff options
-rw-r--r-- | .gitignore | 1 | ||||
-rw-r--r-- | AXtest | bin | 15079072 -> 0 bytes | |||
-rw-r--r-- | AYtest | bin | 15079072 -> 0 bytes | |||
-rw-r--r-- | Gammatest | bin | 15079072 -> 0 bytes | |||
-rw-r--r-- | Ntest | bin | 15079072 -> 0 bytes | |||
-rw-r--r-- | Physikunit.pas | 520 | ||||
-rwxr-xr-x | Plasmapropagation | bin | 3316560 -> 0 bytes | |||
-rw-r--r-- | Plasmapropagation.lps | 113 | ||||
-rw-r--r-- | dAYDTtest | bin | 15079072 -> 0 bytes | |||
-rw-r--r-- | dPhiDXtest | bin | 15079072 -> 0 bytes | |||
-rw-r--r-- | dPsiDXtest | bin | 15079072 -> 0 bytes | |||
-rw-r--r-- | input.plap | 6 |
12 files changed, 437 insertions, 203 deletions
@@ -7,6 +7,7 @@ error *.o *.zip *.tar.gz +Plasmapropagation Ergebnisse lib epost Binary files differBinary files differdiff --git a/Gammatest b/Gammatest Binary files differBinary files differdeleted file mode 100644 index eb6cdf5..0000000 --- a/Gammatest +++ /dev/null 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 Binary files differdeleted file mode 100755 index 561167f..0000000 --- a/Plasmapropagation +++ /dev/null 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 Binary files differdeleted file mode 100644 index 5c8f0b7..0000000 --- a/dAYDTtest +++ /dev/null diff --git a/dPhiDXtest b/dPhiDXtest Binary files differdeleted file mode 100644 index a37f979..0000000 --- a/dPhiDXtest +++ /dev/null diff --git a/dPsiDXtest b/dPsiDXtest Binary files differdeleted file mode 100644 index 7ba5d78..0000000 --- a/dPsiDXtest +++ /dev/null @@ -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 |